Different ST_AZIMUTH results for ST_PROJECT points along the same original bearing

Question

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

Answer 1

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).)