Determine If PostGIS Geometry Is a Circle

Question

I have a PostGIS database table that contains several generic geometries in one of it's columns. Sometimes I have geometries that are intended to be circles which I insert into the database as the result of a center point and a radius point using something like this:

INSERT INTO myTable (drawing) 
VALUES (
    ST_Buffer(
        ST_MakePoint(0, 0)::geography, 
        (SELECT 
            ST_Distance(
                ST_GeographyFromText('POINT(0 0)'), 
                ST_GeographyFromText('POINT(1 0)')
                )
            )
        )::geometry
    )

Now I would like to query the table for all the values inserted in this way, using something like this:

SELECT * FROM myTABLE WHERE FUNCTION_SAYS_GEOMETRY_IS_CIRCLE(geometry)

Does anybody have a recommendation? If this involves inserting my data differently that is okay. I have thought about adding an extra column for "is_circle" or something but that feels like a hack. I would prefer some insight from the geometry itself.

Answer 1

You could try calculating a compactness score for your geometries to see if they are a circle. Something like the Polsby-Popper test will calculate a ratio between 1 and 0, 1 being a perfect circle and any other geometric shape will have a smaller ratio.

4 * pi() * (area/(perimeter^2))

If you are working with perfect circles you can select anything with a ratio of 1.

Here is a blog post on calculating the compactness of congressional districts in PostGIS where they use a few compactness tests one of which is the Polsby-Popper.

https://www.azavea.com/blog/2016/07/11/measuring-district-compactness-postgis/

Answer 2

You can compute a compactness test in a query, but you really don't want to. Here's an example of a small (100k row) table:

DROP TABLE IF EXISTS example1
;

CREATE TABLE example1 (
    idcol       serial      NOT NULL,
    isCircle    char(1)         NULL,
    geomcol     geometry        NULL,
    CONSTRAINT  example1_pk PRIMARY KEY (idcol),
    CONSTRAINT  enforce_srid CHECK (st_srid(geomcol) = 4326)
)
with (
    OIDS=FALSE
);

INSERT INTO example1(isCircle,geomcol)
SELECT isCircle, (CASE WHEN isCircle = 'Y'
        THEN ST_Buffer(point::geography,km*1000)::geometry
        ELSE ST_Envelope(ST_Buffer(point::geography,km*1000)::geometry)::geometry
        END) as geomcol
FROM (
    SELECT  ST_SetSRID(
        ST_MakePoint(
        (random()*360.0) - 180.0,
        (acos(1.0 - 2.0 * random()) * 2.0 - pi()) * 90.0 / pi()),
        4326) as point,
        (CASE WHEN random() > 0.4 THEN 'Y' ELSE NULL END)::char(1) as isCircle,
        floor(random()*50 + 50) as km
    FROM  generate_series(1, 100000) vtab
) vt ;
-- 100000 rows affected, 28.5 secs execution time.

CREATE INDEX example1_spx ON example1 USING GIST (geomcol);
CREATE INDEX example1_ix1 ON example1(isCircle);

-- Query time

SELECT count(*) FROM example1 
WHERE isCircle = 'Y';
-- 60174 rows in 30ms

SELECT count(*) FROM example1 
WHERE 4 * pi() * (ST_Area(geography(geomcol),true) / (ST_Perimeter(geography(geomcol),true) * ST_Perimeter(geography(geomcol),true))) > 0.99;
-- 60174 rows in 18300ms

Computing the isCircle property once per geometry is going to represent a significant savings vice computing it with each query, even if you aren't working in geodetic space. This isn't a hack, just good SQL practice.

---

Addendum: As noted in @JGH's comment below, you can reduce the cost of the geodetic area and perimeter calculation by creating a function and a covering index on that function.

It's also faster to use the pow function to square the perimeter value (at least on my PG 9.5 instance), and using that also drops the index query to effectively identical (140ms down to 32ms)

CREATE FUNCTION compactness( geometry ) RETURNS double precision AS $$
    BEGIN
        RETURN 4.0 * pi() * (ST_Area( $1 ::geography,true) / 
                     pow(ST_Perimeter( $1 ::geography,true),2.0));
    END;
$$ LANGUAGE plpgsql
    IMMUTABLE
    RETURNS NULL ON NULL INPUT;

CREATE INDEX example1_cx1 on example1(compactness(geomcol));
-- 12100ms (due to pow() function)

SELECT count(*) FROM example1 where compactness(geomcol) > 0.99;
--    32ms

Still, if you have the isCircle value up-front, it's not a hack to use it (and if it gets corrupted, you can use the compactness function to recover).