2

I've got a polygon table with a UTM 10x10 km grid in a PostgreSQL/PostGIS database and I need to identify those cells that are not a regular square (i.e. changes in UTM numbers):

enter image description here

Following the comments below, I've built this code to calculate the difference in lenght of two consecutive sides of each cell:

with tabla as (
     select utmcode, 
                ST_Distance(ST_PointN(ST_ExteriorRing(geom), 1),
                ST_PointN(ST_ExteriorRing(geom), 2)) AS d1, 
                ST_Distance(ST_PointN(ST_ExteriorRing(geom), 2),
                ST_PointN(ST_ExteriorRing(geom), 3)) AS d2
     from malla_p)
select utmcode, abs(round(d1-d2)) AS difference         
from tabla          
                ;

But it returns, for example, difference = 0 for this cell:

enter image description here

or difference = 498 for this other:

enter image description here

How can it be solved?

7
  • 1
    Check out this answer: gis.stackexchange.com/a/366837/120129 Dec 11 '21 at 13:49
  • 1
    Do you know where to download your grid? Or did you create it yourself?
    – BERA
    Dec 11 '21 at 18:16
  • 1
    why not just selecting geometries whose area is not 100km2?
    – JGH
    Dec 12 '21 at 13:29
  • 2
    the issue with computing distances between points is that you must be sure that the polygon has only 4 points (+1 for the last point)
    – JGH
    Dec 12 '21 at 13:30
  • 1
    Possible workflow - extract vertices > split exterior rings with vertices > count number of segments/flag polygons with more than 4 > flag polygons with less than 2 matching segment lengths.
    – jbalk
    Dec 13 '21 at 21:53
4

Modified from https://stackoverflow.com/questions/62467829/python-check-if-shapely-polygon-is-a-rectangle.

Select polygons having an area that does not match the area of the oriented envelope https://postgis.net/docs/ST_OrientedEnvelope.html of the polygon.

SELECT * from the_table 
where ST_Area(geometry)!=
ST_Area(ST_OrientedEnvelope(geometry));
2
  • This also selects cells like this one cloud.mail.ru/public/FF7M/xHBBS2LZD that, although being a regular square, are not in a "vertical" position Dec 14 '21 at 8:19
  • You should probably add some tolerance depending on the application - probably don't need full floating point accuracy - where round(st_area(geom)::numeric,1) != round(st_area(st_orientedenvelope(geom))::numeric,1) - adjust the rounding to suit your level of tolerance.
    – jbalk
    Dec 16 '21 at 20:50

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.