Geographic Even Grid using ECEF

ecef_grid.sql

CREATE TABLE global_geodetic_grid_rotation (
    id BIGSERIAL PRIMARY KEY,
    geom GEOMETRY(Point, 4326) NOT NULL
);

INSERT INTO global_geodetic_grid_rotation (geom)
WITH constants AS (
    SELECT
        9000.0::float AS spacing_meters,       -- Target grid spacing
        11.5::float AS rotation_degrees,       -- Desired rotation angle
        6371000.0::float AS earth_radius       -- Average Earth radius
), 
grid_params AS (
    SELECT
        spacing_meters,
        rotation_degrees,
        -- Convert spacing to degrees at equator
        degrees(spacing_meters / earth_radius) AS base_lat_step
    FROM constants
),
-- Create unrotated base grid
base_grid AS (
    SELECT
        -- Calculate unrotated lat/lon for grid point
        -85.0 + (lat_idx * gp.base_lat_step) AS base_lat,
        -180.0 + (lon_idx * (gp.base_lat_step / GREATEST(cos(radians(-85.0 + (lat_idx * gp.base_lat_step))), 0.01))) AS base_lon
    FROM 
        grid_params gp,
        generate_series(0, floor(170.0/gp.base_lat_step)::int) AS lat_series(lat_idx),
        generate_series(0, 
            floor(
                360.0 / (gp.base_lat_step / GREATEST(cos(radians(-85.0 + (lat_idx * gp.base_lat_step))), 0.01))
            )::int
        ) AS lon_series(lon_idx)
    -- WHERE 
    --     lat_idx % 10 = 0 AND lon_idx % 10 = 0  -- COMMENTED OUT: Generate full density grid
),
-- Convert to 3D Cartesian coordinates (ECEF)
cartesian_points AS (
    SELECT
        -- X coordinate (ECEF)
        cos(radians(base_lat)) * cos(radians(base_lon)) AS x,
        -- Y coordinate (ECEF)
        cos(radians(base_lat)) * sin(radians(base_lon)) AS y,
        -- Z coordinate (ECEF)
        sin(radians(base_lat)) AS z
    FROM base_grid
),
-- Apply 3D rotation matrix around Z axis (simplest global rotation)
rotated_cartesian AS (
    SELECT
        -- Apply rotation matrix [cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1]
        x * cos(radians(c.rotation_degrees)) - y * sin(radians(c.rotation_degrees)) AS x_rot,
        x * sin(radians(c.rotation_degrees)) + y * cos(radians(c.rotation_degrees)) AS y_rot,
        z AS z_rot  -- Z remains unchanged in this rotation
    FROM cartesian_points, constants c
),
-- Convert back to geographic coordinates
rotated_geographic AS (
    SELECT
        -- Convert ECEF back to lat/lon
        degrees(atan2(z_rot, sqrt(x_rot*x_rot + y_rot*y_rot))) AS rotated_lat,
        degrees(atan2(y_rot, x_rot)) AS rotated_lon
    FROM rotated_cartesian
)
-- Final query with geometry creation
SELECT
    ST_SetSRID(ST_MakePoint(rotated_lon, rotated_lat), 4326) AS geom
FROM rotated_geographic;

-- Create spatial index
CREATE INDEX global_geodetic_grid_rotation_gist 
ON global_geodetic_grid_rotation USING GIST (geom);