Tuesday, 8 October 2019

How does PostGIS store circles?


I was reading this answer by MikeT to Creating circle in PostGIS? which mentions that the following


SELECT ST_Area(geom)
FROM ( VALUES
('CURVEPOLYGON(CIRCULARSTRING(2 1, 1 2, 0 1, 1 0, 2 1))'::geometry)

) AS t(geom);

returns 3.14033115695475. I tried it myself and it's true. He says it's only a rough approximation of pi? Why is that?



Answer



In PostGIS in Action, 2nd Edition by Regina O. Obe and Leo S. Hsu. They mention this here,



Why is a PostGIS 10-meter buffer of a point 312 and not 314 sq m? If you do the calculation, a perfect 10-meter buffer will give you an area of 10*10*pi(), which is around 314 square meters. The default buffer in PostGIS is a 32-sided polygon (eight points approximate a quarter segment of a circle). You can make this more accurate by using the overloaded version of the ST_Buffer function that allows you to pass in the number of points used to approximate a quarter segment.



It seems that it's not just the ST_Buffer that makes this assumption. ST_Area() calls lwgeom_area_polygon which calls lwgeom_area which eventually calls lwcurvepoly_area. And there again we see it with it's 32-sides.


/**

* This should be rewritten to make use of the curve itself.
*/
double
lwcurvepoly_area(const LWCURVEPOLY *curvepoly)
{
double area = 0.0;
LWPOLY *poly;
if( lwgeom_is_empty((LWGEOM*)curvepoly) )
return 0.0;
poly = lwcurvepoly_stroke(curvepoly, 32);

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...