[postgis-users] Geog Buffering, Distance, and Line contruction

[Lopez, Anthony]

All,

I have been seeing some interesting results using some of the geog functions. First I'll layout what I am attempting to do, my attempted workaround, and the resulting problems.

The goal:
On a global offshore wind dataset of 41 million points, I am attempting to calculate the distance to the closest shoreline (done), and build a line to that location (trouble).
The beginning problem:
There are no geog functions for constructing a line in this manner e.g. No ST_ShortestLine(geog, geog) nor are there geog functions to return the closest point e.g. No ST_ClosestPoint.
My workaround (I thought…):
My workaround was to be quite simple, use ST_Distance(a.geog, b.geog) to get the distance between a point and the closest shoreline, then use the resulting distance in ST_Buffer(geog, returned_dist) to construct a buffered polygon. Lastly, take the buffer and intersect with the the shoreline and return the point of intersection. One limitation with this method is that you cannot specify the number of quad_segs in ST_Buffer(geog), which depending on the buffer/shoreline relationship could potentially produce a near miss…But, I could live with this limitations and could come up with a percentage increase to always produce an intersection.
The larger problem:
My method failed, and when investigating the failure I found these results.

 1.  In comparing ST_Distance(geog) and ST_Distance(geom) the resulting distance were different, staggeringly different depending on latitude and orientation. This WAS expected, even if the two functions measure distance to the same point they would return different measurements due to the nature of projections.
 2.  I tried to determine if constructing a line using ST_ShortestLine(geom) then casting to geog then calculating length would return a similar distance as ST_Distance(geog). It did not, and proves that the geog, geom functions are not measuring distance to the same point (on the shoreline). I'll mention that I expected this as well.
 3.  When buffering by the distance ST_Distance(geog) returns, using ST_Buffer(geog), you would expect a point of intersection between the shoreline and the buffer. I understand that since ST_Buffer(geog) does not have the quad_segs option it sets up the possibility of a near miss (hit). It seems that ST_Buffer(geog) comes up short at all latitudes and orientations…

Please see my testing SQL below. I am completely up for other workarounds and suggestions.

Again, the goal is to be able to construct lines that represent the same distance/point of intersect that ST_Distance(geog) uses.

Thanks in advance,
Anthony

INSERT INTO global_wind.buff_tests
SELECT
buff.gid,
buff.country,
buff.lat,
buff.geom_dist / 1000 as geom_dist,
buff.geog_dist / 1000 as geog_dist,
buff.geom_line,
Geometry(buff.geog_buff) as buffer_geom,
Geometry(ST_Intersection(buff.geog_buff, buff.cntry_geog)) as buff_cntry_inter
ST_Intersects(buff.geog_buff, buff.cntry_geog)

FROM (

SELECT
r.gid,
r.country,
r.lat,
r.geom_dist,
r.geog_dist,
ST_Buffer(r.pt_geog, r.geog_dist) as geog_buff,
r.geom_line,
r.cntry_geog

FROM (

SELECT
a.gid,
a.country,
ST_Y(a.the_geom) as lat,
ST_Length(Geography(ST_ShortestLine(a.the_geom, b.the_geom))) as geom_dist,
ST_Distance(a.the_geog, b.geog) as geog_dist,
ST_ShortestLine(a.the_geom, b.the_geom) as geom_line,
a.the_geog as pt_geog,
b.geog as cntry_geog
FROM
global_wind.offshore_wind a, global_wind.coasts_simp b
WHERE
a.gid = 22000000 AND a.eez_id = b.zones ) r ) buff