PostGIS - How do I create a line that directly connects two earth points without travelling over a geodesic?

Question

In PostGIS 2.1.1, I've noticed that if I create two distant earth points with equal altitude and then create a line connecting them, as I gather interpolated Z coordinate data along the line my Z (height) results remain equal to my original altitude_in_meters.

I've confirmed that it doesn't simply ignore the Z value because I created a line between an sp and ep with different heights, and I get different Z values along the line, indicating slope from sp to ep. I've also confirmed that these lines are also still curved such that they will always lie above the geoid and never cut through it, which should happen at very great distances with an actual straight line.

This tells me that the line being projected between sp and ep is actually a geodesic. How can I create an actual straight line connecting these points in 3D space? My desired result is to see the relative decrease in altitude along the line connecting two distant earth points, which would be most prevalent at its midpoint.

Here's some SQL to show how ST_MakeLine() is creating a geodesic:

WITH
points AS(
 SELECT
   ST_SETSRID(ST_MAKEPOINT(-74.5, 40.5, 50), 4326) AS sp_equal_height
   ,ST_SETSRID(ST_MAKEPOINT(-74.45, 40.45, 50), 4326) AS ep_equal_height
   ,ST_SETSRID(ST_MAKEPOINT(-74.5, 40.5, 100), 4326) AS sp_diff_height
   ,ST_SETSRID(ST_MAKEPOINT(-74.45, 40.45, 50), 4326) AS ep_diff_height
 )
,line AS(
 SELECT
   ST_MAKELINE(sp_equal_height, ep_equal_height) AS line_equal_height
   ,ST_MAKELINE(sp_diff_height, ep_diff_height) AS line_diff_height
 FROM
   points
 )
SELECT
  ST_Z(ST_LineInterpolatePoint(line_equal_height, 0.5)) AS mid_equal_height_z
  ,ST_Z(ST_LineInterpolatePoint(line_diff_height, 0.5)) AS mid_diff_height_z
FROM
  points
  ,line

Results: mid_equal_height_z = 50 and mid_diff_height_z = 75

Answer 1

What you're trying to do can't be done in normal projection, we need a cartesian coordinates for that. With a cartesian coordinate system we can cut straight thru the Earth instead of working on it's surface Projection 4978 is just what we need [​1], however the code from the site needs to be fixed [​2].

This code adds the projection we need:

DELETE FROM spatial_ref_sys WHERE srid = 94978;
INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 94978, 'epsg', 4978, '+proj=geocent +datum=WGS84 +units=m +no_defs
', 
'GEOCCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["m",1.0],AXIS["Geocentric X",OTHER],AXIS["Geocentric Y",EAST],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4978"]]');

So what we need to do is take your points with 4326 projection, make lines, interpolate needed points and transform then back (we could transform the lines instead of initial points but still we need to do interpolation before transformation).

WITH
points AS(
 SELECT
   ST_Transform(ST_SETSRID(ST_MAKEPOINT(-74.5, 40.5, 50), 4326),94978) AS sp_equal_height
   ,ST_Transform(ST_SETSRID(ST_MAKEPOINT(-74.45, 40.45, 50), 4326),94978) AS ep_equal_height
   ,ST_Transform(ST_SETSRID(ST_MAKEPOINT(-74.5, 40.5, 100), 4326),94978) AS sp_diff_height
   ,ST_Transform(ST_SETSRID(ST_MAKEPOINT(-74.45, 40.45, 50), 4326),94978) AS ep_diff_height
 )
,line AS(
 SELECT
   ST_MakeLine(sp_equal_height, ep_equal_height) AS line_equal_height
   ,ST_MakeLine(sp_diff_height, ep_diff_height) AS line_diff_height
 FROM
   points
 )
SELECT
  ST_Transform(ST_LineInterpolatePoint(line_equal_height, i/100::float),4326) AS mid_equal_height_z
  ,ST_Transform(ST_LineInterpolatePoint(line_diff_height, i/100::float),4326) AS mid_diff_height_z
FROM
  line
  ,generate_series(0,100) as i

I also add more points so we could see it on the graph:

• red - original equal
• green - cartesian equal
• yellow - original diff
• blue - cartesian diff