How to create a circle in meters in postgis?

Question

I would like to ask how to create a circle with radius=4km. I have tried the ST_Buffer function but it creates a larger circle. (I see the created circle by inserting its polygon into an new kml file.)

This is what i am trying.

INSERT INTO camera(geom_circle) VALUES(geometry(ST_Buffer(georgaphy(ST_GeomFromText('POINT(21.304116745663165 38.68607570952619)')), 4000)))

The center of the circle is a lon lat point but I don't know its SRID because I have imported it from a kml file. Do I need the SRID in order to transform the geometries etc?

Mike T · Accepted Answer · 2019-03-14T20:41:24.760

KML files are always lat/long and use SRID=4326. This SRID is implied if you use geography. Geography is a good way to mix-in the 4 km metric measure on lat/long data ... excellent you tried this!

Try this statement to fix up the casts, and use a parameterized point constructor:

SELECT ST_Buffer(ST_MakePoint(21.304116745663165, 38.68607570952619)::geography, 4000);

And if you need to cast this back to geometry, add a ::geometry cast to the end.

Update on accuracy

The previous answer internally re-projects the geometry (usually) to a UTM zone that the point fits within (see ST_Buffer). This may cause minor distortions if the point is on the edge of two UTM boundaries. Most folks won't care about the size of these errors, but it will often be several meters. However, if you require sub millimeter precision, consider building a dynamic azimuthal equidistant projection. This requires PostGIS 2.3's ST_Transform, and is adapted from another answer:

CREATE OR REPLACE FUNCTION geodesic_buffer(geom geometry, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geometry AS $$
  SELECT ST_Transform(
    ST_Buffer(ST_Point(0, 0), $2, $3),
      ('+proj=aeqd +x_0=0 +y_0=0 +lat_0='
       || ST_Y(ST_Centroid($1))::text || ' +lon_0=' || ST_X(ST_Centroid($1))::text),
      ST_SRID($1))
  $$ LANGUAGE sql IMMUTABLE STRICT COST 100;
CREATE OR REPLACE FUNCTION geodesic_buffer(geom geometry, dist double precision)
  RETURNS geometry AS 'SELECT geodesic_buffer($1, $2, 8)'
  LANGUAGE sql IMMUTABLE STRICT COST 100;
-- Optional warppers for geography type
CREATE OR REPLACE FUNCTION geodesic_buffer(geog geography, dist double precision)
  RETURNS geography AS 'SELECT geodesic_buffer($1::geometry, $2)::geography'
LANGUAGE sql IMMUTABLE STRICT COST 100;
CREATE OR REPLACE FUNCTION geodesic_buffer(geog geography, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geography AS 'SELECT geodesic_buffer($1::geometry, $2, $3)::geography'
  LANGUAGE sql IMMUTABLE STRICT COST 100;

A simple example to run one of the functions is:

SELECT geodesic_buffer(ST_MakePoint(21.304116745663165, 38.68607570952619)::geography, 4000);

And to compare the distances to each of the buffered points, here are the lengths of each geodesic (shortest path on an ellipsoid of revolution, i.e. WGS84). First this function:

SELECT count(*), min(buff_dist), avg(buff_dist), max(buff_dist)
FROM (
  SELECT ST_Distance((ST_DumpPoints(geodesic_buffer(poi, dist)::geometry)).geom, poi) AS buff_dist
  FROM (SELECT ST_MakePoint(21.304116745663165, 38.68607570952619)::geography AS poi, 4000 AS dist) AS f
) AS f;

 count |      min       |       avg       |      max
-------+----------------+-----------------+----------------
    33 | 3999.999999953 | 3999.9999999743 | 4000.000000001

Compare this to ST_Buffer (first part of answer), that shows it's off by about 1.56 m:

SELECT count(*), min(buff_dist), avg(buff_dist), max(buff_dist)
FROM (
  SELECT ST_Distance((ST_DumpPoints(ST_Buffer(poi, dist)::geometry)).geom, poi) AS buff_dist
  FROM (SELECT ST_MakePoint(21.304116745663165, 38.68607570952619)::geography AS poi, 4000 AS dist) AS f
) AS f;

 count |      min       |       avg        |      max
-------+----------------+------------------+----------------
    33 | 4001.560675049 | 4001.56585986067 | 4001.571105793

This is a great answer and worked for me. If you try using only geometry for your calculations, depending on where you are at in the world, you'll end up with an ellipse. A good way to test the information you're storing in your geometry column is to select ST_AsText(geoshape) AS geom and render it on your map. — Ryan Charmley, Oct 07 '13 at 19:39
@RyanCharmley although it may look like an ellipse in Cartesian space (a flat map), these are always circles on a sphere. For instance, view the output in Google Earth and you should see circles, even away from the equator. — Mike T, Oct 08 '13 at 00:44
Hello there, long ago, but the question is still relevant. The approach is not so accurate because what you have is still an inscribed polygon, which is an approximation of a circle. `st_segmentize` on neighboring points, taking a middle point and then `st_distance` with the circle center gives a whopping 940 meters of difference for a 200 km circle, if I call your function without passing `num_seg_quarter_circle`. I haven't found anything which can give me the real geography circle (or the less approximated, but still performant one), so I'm going to mess around with `st_dwithin` instead. — 1valdis, Sep 20 '19 at 12:26
@1valdis Correct, this answer provides a polygon approximation of a circular geometry, thus will have discrepancies. There is probably a way to describe this circular geometry in WKT, however in practice, circle geometries are seldom used with most GIS software today. And yes, `ST_DWithin` is the correct and most efficient approach to find geometries within an exact distance of a radius from another point. — Mike T, Sep 23 '19 at 01:37

How to create a circle in meters in postgis?

1 Answers1

Update on accuracy

Linked