Different ST_AZIMUTH results for ST_PROJECT points along the same original bearing

I'm trying to develop some way of determining line of sight (LOS) between two geographical points about 1km apart or less using ST_3DINTERSECTS or something similar.

BASIC PROBLEM:

Along the way, I've found that if I calculate the ST_AZIMUTH between two points and then use that bearing in a series of ST_PROJECT statements to create additional points along the line connecting the two original points, I get a close but slightly different ST_AZIMUTH when calculating from the original start point out to each of the new points along the line.

I believe this slight difference in bearings between the lines is causing my ST_3DINTERSECTS function call to return false for all scenarios, even when terrain elevation should be blocking LOS between the two original points.

DETAILS:

To begin, I have a start point (SP) and an end point (EP) both represented as geographies. I calculate the original bearing: (OB) = ST_AZIMUTH(SP, EP). Then I use my STEP size (STEP = 30 meters), SP and OB to create points along the line connecting SP and EP in order to pull elevation samples along the way: ST_PROJECT(SP, STEP, OB). In this manner, I now have roughly 30 geography points with elevations along the line connecting SP and EP.

In all cases, elevation data is retrieved using my raster data in conjunction with these geography points casted as geometries e.g.: ELEV = ST_VALUE(myraster, SP::geometry). I also extract the X and Y coordinates from the geography points casted as geometries e.g.: X = ST_X(SP::geometry) and Y = ST_Y(SP::geometry). Using this information, I create a series of 3-dimensional points e.g.: ST_MAKEPOINT(X, Y, ELEV) in order to represent the original SP and EP as well as all sample points along the line connecting SP and EP.

Using these 3D points, I create two lines e.g.: ST_MAKELINE(my3dpoints). One is a simple line connecting just the SP and the EP in 3D. The other is a segmented line comprised of all the sample points taken every 30 meters along the straight line connecting SP and EP. When I use ST_3DINTERSECTS on both of these lines, it always returns a false result, even if I know the terrain should be obstructing LOS.

QUESTION:

Can anyone verify if the slight bearing difference is what's causing my lines not to intersect when they should? Is there a better way to do this?

SETUP:

Example SP = POINT(-74.714919, 41.282384), EP = POINT(-74.7147106468, 41.2913868143398) SP is always 3 meters above the actual terrain elevation. EP is always 2 meters above the actual terrain elevation. These two points should not have LOS due to terrain.

My elevation data source is: http://www.viewfinderpanoramas.org/dem3/K18.zip loaded using raster2pgsql with SRID 4326.

My postgis_full_version() is: POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER

• You are not doing what you think you are with your repeated use of the initial azimuth. You're treating it like a loxodrome when it's really an orthodrome gis.stackexchange.com/a/112261/457 – Paul Ramsey Sep 3 '14 at 15:23
• @PaulRamsey, thanks for the help. I wasn't familiar with that terminology. I just created a large-scale demonstration to help me visualize how these two lines don't intersect. Do you think it would be appropriate to calculate for intersection by converting the direct line into a polygon by projecting corner points a certain distance away from the original SP and EP on a bearing perpendicular to the OB? – Shawn Sep 4 '14 at 15:24

You are creating a specific example of the general observation that introducing vertices into a geometry tends to make it not the same as the original.

SELECT ST_Equals(ST_Segmentize(geom,10),geom);

Why is this so? Because it's possible (heck, it is common) for the "lines" between vertices to traverse real numbers that are not representable using double precision floating point numbers. Doubles have a finite number of bits to hold information; real numbers have infinite precision. When you create a new vertex in a line, odds are that it is almost on the line, but not quite exactly, because the line is traversing some portion of cartesian space that just is not exactly representable using double precision numbers.

What is an alternate approach? Well, you're creating points at fixed intervals from SP to EP. Why not, in addition to grabbing their elevation from the DEM, also calculate the Z value of the interpolated line from SP to EP at that point: if the DEM Z is > than the interpolated Z, you know the LOS is obstructed and you stop; otherwise, move on to the next point.

(Another approach, similar to your original one, would be to substitute 3D intersects with 3d distance, and test for distance being almost zero, which would indicate a crossing point (or near enough).)

• Makes sense. But doesn't Postgres allow for infinite precision in its 'numeric' data type? I wonder why that isn't utilized in these functions. – Shawn Sep 4 '14 at 17:32
• That would involve having all the vertices represented w/ infinite precision, which would make all calculations much slower, and all geometries much larger. A "better" approach would be allowing all predicates to work with tolerances, though that's such a huge lift as to be unlikely to ever happen. – Paul Ramsey Sep 4 '14 at 19:58
• It sure makes for some interesting problem solving. In any case, I'm trying to employ a method of calculating the Z value of the interpolated line, but I'm struggling since we're not just talking about simple triangles here. The bottom of the 'triangle' is an orthodrome. Isn't the interpolated line also curved somehow? I tested this theory on a large-scale model and while Google Earth draws a straight direct line connecting two points, I used ST_LINEINTERPOLATEPOINT to discover that the Postgis direct line curves around the earth in quite a funny way. – Shawn Sep 5 '14 at 11:37
• I feel foolish and elated at the same time for having just discovered the wonders of ST_Z().... – Shawn Sep 5 '14 at 11:58