Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

ST_Buffer on Geography — Iteration

Posted by smathermather on January 12, 2016

In my previous post, Imagico said:

This is nice but still fails in one important aspect – it will create significant inaccuracies once your buffering distance reaches values where buffering in a single local coordinate system is getting problematic which, depending on your precision requirements, will likely happen at distances of a few degrees. A solution to this would require either an iterative split and merge technique of a completely different approach on a very basic level (i.e. attempting truly spherical buffering).

Ok, ok. I kinda forgot about this aspect. The proper way to do this is, as always is to do it at a lower level thank in an SQL function, preferably calculating the buffer on the spheroid, ellipsoid, or geoid. But, ouch. I think we can munge something good enough for now.

We’ll still do the buffer:

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
...

But now that we’ve performed the buffer, we need to determine that the buffered geometry is not too large for our local coordinate system. To do that, we’ll compare the best SRID for our buffered geometry to the original best SRID that we calculated based on the original input geometry.

SELECT
CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography

If the two SRIDs aren’t equal, then we need to use the new best SRID based on buffered geometry and buffer using that new best SRID.

ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography

Now let’s put it all together:

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
buffer_style_parameters text)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, buffer_style_parameters), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;

So, what’s the affect? Pretty significant really.

SELECT ST_Buffer(
	 ST_GeomFromText(
	  'POINT(0 0)'
	 )::geography, 5000000, 50);
Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

2 Responses to “ST_Buffer on Geography — Iteration”

  1. imagico said

    Great example – since, as you probably have not realized, the blue and the red line are equally wrong (not necessarily exactly equally but similarly). One of them is stretched vertically, one of them horizontally: blue reaches north to northern Spain, less than 42°N, red up to southern France, more than 44°N and red reaches east to about the eastern end of Kenya (less than 42°E) while blue goes way into Somalia, more than 44°E. This is because one uses UTM and the other one uses World Mercator which is more or less the same projection, just rotated 90 degrees.

    _ST_BestSRID() suggest UTM if the geometry is small enough and not at high latitudes. At high latitudes it suggests Lambert Azimuthal Equal Area and in all other cases it suggests World Mercator. This is only better than UTM as far as it traditionally uses a more accurate formula that works also at larger distances from the equator (which is why both circles could likely be wrong in different ways – the error in UTM is a combination of geometric and numerical problems) but it has the same variation of scale that leads to incorrect buffering of larger objects. I wrote some more about this some time ago in http://blog.imagico.de/projection-issues-and-tools-to-deal-with-them/.

    As usual: there is no such thing as a free lunch…

    • Re: free lunches, of course there aren’t. We should be calculating all buffers on a dynamic gravimetric geoid and should require the user to input a date time code… :).

      I’ll work on a version with ST_SubDivide for further optimization but golly this is a deep hole.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: