2

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

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:

height graph

  • red - original equal
  • green - cartesian equal
  • yellow - original diff
  • blue - cartesian diff
  • Thanks! I guess I'm surprised they didn't have a cartesian coordinate system in the default spatial_ref_sys table. I'm curious though, why add this as 94978 instead of just 4978? – Shawn Sep 22 '14 at 14:14
  • @Shawn That is a very good question but I don't know to be honest. spatialreference.org supplies data in such a way, probably to avoid collision with the default PostGIS data, but I don't know if there is some signifance to 9 as a prefix. – Jakub Kania Sep 22 '14 at 14:27
  • Also, I just want to point out that there's a space in the portion of text before +proj=geocent and the preceding single quote. I noticed that none of the other entries in the table have that space. Is that a typo? Does it matter? – Shawn Sep 22 '14 at 15:30
  • Yeah, that's a messy cut and paste, it works though so it doesn't seem to matter. – Jakub Kania Sep 22 '14 at 17:50
  • So I've gone ahead and implemented this by creating 3D points along a line in SRID 4326 using Z values as elevation data from my raster, then transforming them to 94978 and rotated everything to normalize. At short distances ~1km, things look pretty good and very slightly curved, but I also did a quick test at 140km with a 58m starting elevation and a 355m ending elevation, and my graph of the cartesian representation shows a curved elevation drop of over 1200m along the 140km path. Estimates that I've researched indicate that I shouldn't see more than 20m of drop over 140km. Any ideas? – Shawn Sep 24 '14 at 11:31

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.