# ST_Buffer on Geography — Iteration

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);
```

## 2 thoughts on “ST_Buffer on Geography — Iteration”

1. 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…

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

This site uses Akismet to reduce spam. Learn how your comment data is processed.