2

In PostGIS 2.1.1, I've been plotting elevations from a raster over distances, but noticed that beyond 1km, the earth's curvature seems to be greatly exaggerated. In this example over a 40km geodesic entirely at 0 altitude (for simplification), the earth appears to drop almost 130m below the starting point. From estimating on the great circle, I would expect this drop to be around only 3m.

130m drop over 40km geodesic

This problem is simplified so as not to use raster data until I can figure out this curvature problem. What I've done is taken two lat/lon points on the earth's surface and created a geodesic line connecting them. Then I interpolated 100 points along that line and assigned a 0 altitude to all of those points. I also added an additional preceding point directly above the first point by 3m, then did the same following the last point, raised by 2m. These points represent antennas, which should be projecting directly above the start/end points.

With all of this in lat/long/altitude (SRID 4326), I transformed those points to cartesian (SRID 94978) per PostGIS - How do I create a line that directly connects two earth points without travelling over a geodesic?

Then I translated these points to work around (0,0,0) using the first ground (0 altitude) level point as the center. To normalize the image for my plot, I then used the 3m starting antenna point as a reference to snap the points upright per the Y axis by rotating along the XY (ST_RotateZ) and ZY (ST_RotateX) planes.

Finally, I had to rotate along the XZ plane (ST_RotateY) in order for the whole geodesic line to project flat against the X axis, except it actually shows a bulge approaching 30m towards the middle of the geodesic, which I'm hoping someone can explain to me (see Z column from output of the SQL below).

Plot X as distance and Y as elevation from the below SQL in order to duplicate my image, above.

WITH
ll AS(
SELECT
    POINT(-74.35, 40.5) AS startpt 
    ,POINT(-74.7, 40.25) AS endpt 
    ,3 AS start_antenna_height 
    ,2 AS end_antenna_height
)
,geog_ll AS(
SELECT
    ST_SETSRID(ST_MAKEPOINT(startpt[0], startpt[1]), 4326)::geography AS startloc
    ,ST_SETSRID(ST_MAKEPOINT(endpt[0], endpt[1]), 4326)::geography AS endloc
FROM
    ll
)
,geodesic AS(
SELECT
    ST_MAKELINE(startloc::geometry, endloc::geometry) AS geo_line
FROM
    geog_ll
)
,geodesic_ll AS(
SELECT
    ST_LINEINTERPOLATEPOINT(geo_line, i/100::float) AS geo_ll
FROM
    geodesic
    ,generate_series(0,100) AS i
)
,geodesic_lla AS(
SELECT
    ST_DISTANCE(startloc, geo_ll::geography, true) AS geodistance
    ,ST_SETSRID(ST_MAKEPOINT(ST_X(geo_ll), ST_Y(geo_ll), 0), 4326) AS geo_lla
FROM
    geodesic_ll
    ,geog_ll
)
,antennas_and_geodesic_lla AS(
SELECT
    0 AS draw_order
    ,ST_SETSRID(ST_MAKEPOINT(ST_X(startloc::geometry), ST_Y(startloc::geometry), start_antenna_height), 4326) AS geo_lla
FROM
    geog_ll
    ,ll
UNION ALL
SELECT
    1 AS draw_order
    ,geo_lla
FROM
    geodesic_lla
UNION ALL
SELECT
    2 AS draw_order
    ,ST_SETSRID(ST_MAKEPOINT(ST_X(endloc::geometry), ST_Y(endloc::geometry), end_antenna_height), 4326) AS geo_lla
FROM
    geog_ll
    ,ll
)
,antennas_and_geodesic_lla_and_cart AS(
SELECT
    ST_DISTANCE(startloc, geo_lla::geography, true) AS geodistance
    ,draw_order
    ,geo_lla
    ,ST_TRANSFORM(geo_lla, 94978) AS geo_cart
FROM
    antennas_and_geodesic_lla
    ,geog_ll
)
,start_ground_point AS(
SELECT 
    geo_lla AS start_ground_lla
    ,geo_cart AS start_ground_cart
FROM 
    antennas_and_geodesic_lla_and_cart 
WHERE 
    geodistance = 0
    AND draw_order = 1
)
,translated_geodesic_cartesian AS(
SELECT  
    geodistance
    ,draw_order
    ,ST_TRANSLATE(geo_cart, -ST_X(start_ground_cart), -ST_Y(start_ground_cart), 
FROM
    antennas_and_geodesic_lla_and_cart
    ,start_ground_point
)
,translated_start_antenna_point AS(
SELECT
    trans_geo_cart AS t_start_antenna_cart
FROM
    translated_geodesic_cartesian
WHERE
    draw_order = 0
)
,z_rotated_geodesic_cartesian AS(
SELECT
    geodistance
    ,draw_order
    ,ST_ROTATEZ(trans_geo_cart, CASE WHEN ST_Y(t_start_antenna_cart) > 0 THEN xy_rads_from_y ELSE pi() + xy_rads_from_y END) AS zr_geo_cart
FROM
    translated_geodesic_cartesian
    ,translated_start_antenna_point
    ,atan(ST_X(t_start_antenna_cart) / ST_Y(t_start_antenna_cart)) AS xy_rads_from_y
)
,z_rotated_start_antenna_point AS(
SELECT
    zr_geo_cart AS zr_start_antenna_cart
FROM
    z_rotated_geodesic_cartesian
WHERE
    draw_order = 0
)
,xz_rotated_geodesic_cartesian AS(
SELECT
    geodistance
    ,draw_order
    ,ST_ROTATEX(zr_geo_cart, -zy_rads_from_y) AS xzr_geo_cart
FROM
    z_rotated_geodesic_cartesian
    ,z_rotated_start_antenna_point
    ,atan(ST_Z(zr_start_antenna_cart) / ST_Y(zr_start_antenna_cart)) AS zy_rads_from_y
)
,xz_rotated_end_ground_point AS(
SELECT
    xzr_geo_cart AS xzr_end_ground_cart
FROM
    xz_rotated_geodesic_cartesian
WHERE
    draw_order = 1
ORDER BY
    geodistance DESC
LIMIT 1
)
,xyz_rotated_geodesic_cartesian AS(
SELECT
    geodistance
    ,draw_order
    ,ST_ROTATEY(xzr_geo_cart, CASE WHEN ST_X(xzr_end_ground_cart) < 0 THEN pi() + zx_rads_from_x ELSE zx_rads_from_x END) AS xyzr_geo_cart
FROM
    xz_rotated_geodesic_cartesian
    ,xz_rotated_end_ground_point
    ,atan(ST_Z(xzr_end_ground_cart) / ST_X(xzr_end_ground_cart)) AS zx_rads_from_x
)
SELECT
    geodistance
    ,draw_order
    ,ST_X(xyzr_geo_cart) AS x
    ,ST_Y(xyzr_geo_cart) AS y
    ,ST_Z(xyzr_geo_cart) AS z
FROM
    xyz_rotated_geodesic_cartesian
ORDER BY
    geodistance
    ,draw_order
1

Apparently my underlying assumption of expecting a 3m drop over 40km was incorrect. I had taken an answer I found stating that the earth dropped off by .078m over 1km, and I attempted to apply it linearly over 40km. But due to the nature of this problem, the drop over 1km does not linearly scale as the geodesic distance increases.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.