6

I'm trying to create circles of equal area in different places, but depending on the location of the centroid different areas are calculated. For example

SELECT ST_Area(ST_Buffer('POINT(0 0)'::geography, 100000))/1000000

results in an area of 31,291 km², whereas

SELECT ST_Area(ST_Buffer('POINT(30 30)'::geography, 100000))/1000000

results in an area of 31,233 km².

Is this the correct approach to take for this problem, and if so, is there a method that would allow the creation of circles of equal area?

3
  • did you mean to use to coordinate examples? Commented Mar 10, 2016 at 20:06
  • You seem to have used the same centroid in both examples?
    – Andy
    Commented Mar 10, 2016 at 21:12
  • Sorry, I meant to use different centroids. Updated to reflect that.
    – jczaplew
    Commented Mar 10, 2016 at 21:54

2 Answers 2

11

Internally, ST_Buffer(geography, ...) uses a fixed projection guess with _ST_BestSRID, which are typically UTM zones or whatever makes sense to the algorithm. This is why you see the differences, because they are different projections that are not optimized for the location of the points.


For simple point buffers, you could use a custom azimuthal equidistant projection, which are better suited for unique points around the globe.

For PostGIS 2.2 and earlier, get a new ST_Transform function to use custom projections:

CREATE FUNCTION ST_Transform(geom geometry, from_proj text, to_srid integer)
  RETURNS geometry AS
'SELECT postgis_transform_geometry(ST_SetSRID($1, 99), $2, proj4text, $3)
FROM spatial_ref_sys WHERE srid=$3;'
  LANGUAGE sql IMMUTABLE STRICT;

And use it to determine buffers of various places:

SELECT ST_Area(
    ST_Transform(ST_Buffer(ST_MakePoint(0, 0), 100000.0),
        concat('+proj=aeqd +lat_0=', ST_Y(ST_Centroid(geom)),
               ' +lon_0=', ST_X(ST_Centroid(geom)),
               ' +x_0=0 +y_0=0 +ellps=WGS84 +no_defs '),
        4326)::geography) / 1000000 AS area_km2
FROM (
     SELECT 'SRID=4326;POINT(0 0)'::geometry geom
     UNION SELECT 'SRID=4326;POINT(30 30)'::geometry
 ) s;
      area_km2
------------------
 31213.8449025783
 31213.8469309491
(2 rows)

I hope a 0.002 km² difference is acceptable. You can also experiment with Lambert azimuthal equal-area projections by using +proj=laea for the PROJ text; I get differences as low as 0.00007 km².

2
  • The above works perfectly on PostGIS 2.2, but on 2.1 I get the following areas: 31353.9839957326 and 31213.8099507586. Why is that?
    – jczaplew
    Commented Mar 10, 2016 at 22:50
  • Yes, this was expected. PostGIS 2.1 and before use a Vincenty-based technique, while 2.2 and after use a more precise approach from Karney (2013) with GeographicLib
    – Mike T
    Commented Mar 10, 2016 at 23:01
1

If your data is in a geographic coordinate system, rather than a cartesian system, your buffers could vary in area due to the shape of the Earth.

58 square kilometers in area difference may be possible if say one buffer is at the equator, and the other is near to a pole.

In WGS84, the Earth flattens approximately 1 meter for every 298 meters as you travel from the equator to a pole. This would make your buffer elliptical if you are calculating the buffer on the geodesic.

I would make certain that the buffer radii are the same first though.

Then, I would save the database to something that is projected rather than geographic, and try to see if the areas are closer to one another.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.