4

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.

  • 1
    Calculating the area/perimeter ratio might actually be more expensive than adding a char(1) column. – Vince Jul 22 at 22:38
  • Buffers are not perfect circles but only a more or less rough approximation. Also beware of the projection you use (a circle of 1 degree is not a circle if measured in meters.) That being said, the approach from @OClark could be complemented by a filtering on the shape bounding box X to Y ratio (should be close to a square) – JGH Jul 23 at 1:42
5

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/

  • As @Vince points out, incorporating this alongside an isCircle field would save a lot of processing. I would look at putting this function into a trigger to automatically update the isCircle field after an individual row is inserted or updated. – Nate Wanner Jul 23 at 14:03
  • Keep in mind that ST_Buffer does not create true circles. Most "circles" in GIS will actually be a many-sided polygon, and ST_Buffer defaults to creating a 32-sided regular polygon. (Though PostGIS does support circular arcs, using it is still pretty niche.) The Polsby-Popper score will be 0.996785… – Lee Hachadoorian Aug 6 at 20:43
4

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

  • 2
    Turning the equation into a function and creating an index on it CREATE INDEX example1_ix2 ON example1(compactness(geomcol)); would get the same benefit, while allowing playing with the threshold value (SELECT count(*) FROM example1 WHERE compactness(geomcol) > 0.99;) and not having to handle an extra column (that would eventually be updated to a wrong value by someone, one day) – JGH Jul 23 at 21:55
  • Not exactly the same benefit, since the Y/null index is going to be faster than the double precision covering index (30ms v. 140ms), and you have to pay that function execution cost build the index at some point(s), even if it isn't on every query. – Vince Jul 23 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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