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