9

I have a set of polygons representing large areas, say city neighborhoods. I want to identify the large overlapping areas between them.

But there's a problem: sometimes these polygons will overlap along their perimeters (because they were drawn with little precision). This will generate long and narrow overlaps that I do not care about.

But other times there will be big overlaps of robust polygons, meaning large areas where a neighborhood's polygon overlaps another. I want to select only these.

See the picture below of just the overlaps. Imagine I wanted to select only the blue polygon in the lower left corner.

overlap

I could look at areas, but sometimes the narrow ones are so long they end up having areas as large as the blue polygon. I've tried to do a ratio of area / perimeter, but that has also yielded mixed results.

I've even tried using ST_MinimumClearance, but sometimes the large areas will have a narrow part attached to it, or two very close vertices.

Any ideas of other approaches?


In the end what worked best for me was using a negative buffer, as suggested by @Cyril and @FGreg below.

I used something like:

ST_Area(ST_Buffer(geom, -10)) as neg_buffer_area

In my case, units were meters, so 10 m negative buffer.

For narrow polygons, this area returned zero (also, the geometry would be empty). Then I used this column to filter out the narrow polygons.

5

I would try to create a negative buffer, if it eats thin polygons, then it’s good, if it doesn’t eat the polygon, then it’s mine ... :-)

run this script, having previously set 2/3 of the width of the linear polygons ...

create table name_table as
SELECT ST_Buffer(
(ST_Dump(
(ST_Union(
ST_Buffer(
(geom),-0.0001))))).geom,
0.0001)) as geom from source_table

OS :-)...

  • in the end your suggestion is what worked best for me. I ended using something like ST_Area(ST_Buffer(geom, -10)), the -10 being -10 meters in my case. If anything returned 0 from that expression then I could filter it out. – bplmp Mar 21 at 14:11
  • Thank you, I am glad that it helped you, use this approach in such cases, with respect ... – Cyril Mar 21 at 14:12
9

Instead of area/perimeter, it is better to use the area divided by the square of the perimeter (or its inverse).

This is also called "shape index". The square of the perimeter divided by the area has a minimum value of 4*Pi() (in the case of a disk, which is the most compact 2D geometry), so it can be normalized by 4*Pi() for an easy interpretation (normalized values close to 1 then mean that you have very compact objects and squares have a values of approximately 1.27).

EDIT: A threshold on the area would be usefull to remove the very small artefacts, which could be compact. Then the shape index would show better contrast. EDIT: in addition to this answer, the use of ST_Snap could help you solve the problem before it occurs.

  • Thanks! But I'm unsure how ST_Snap could help in this case... If I got it right, you're suggesting something like (o.overlap_perimeter^2 / o.overlap_area) / (4 * Pi()) as overlap_ratio? This is having worse results for me than just area / perimeter. – bplmp Mar 20 at 15:55
  • Now using o.overlap_perimeter / (4 * sqrt(o.overlap_area)) as overlap_ratio according to this paper, but still worse results (although that's hard to quantify what I mean by worse) isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/135/…, page 183. – bplmp Mar 20 at 16:09
  • 2
    Thank you for this, I had never heard of the "shape index". I had always thought that using a minimum bounding rectangle was the best way to answer this sort of question. I found this, repository.asu.edu/attachments/111230/content/…, which is interesting. – John Powell Mar 20 at 16:50
  • @JohnPowell intersting paper, thanks. I see that what I know as a shape index is called circularity index in the paper. My problem with minimum bounding rectangles is that it doesn't work with very concave objects (e.g. U-shaped) – radouxju Mar 21 at 7:15
  • @bplmp ST_Snap would help you snap the vertices of "nearly" adjacent polygons so that they do no overlap anymore. There is no scale on your figures, but your artefact look like lines, so I guess that you can use a tolerance value theat is enough to avoid artefacts but does not affect the large polygons. – radouxju Mar 21 at 7:20
5

One option would be to use the ratio of the area of the polygon to the longest line that can be drawn using its extremities. Identifying long narrow polygons.

select * from polygons where ST_Length(ST_LongestLine(geom, geom)) < ST_Area(geom) * 4

This works pretty well for sliver polygons. You can adjust what the ratio (what you multiply the area with) to suit your needs and projection.

1

It sounds like this might match your use case: Eliminate selected polygons

Combines selected polygons of the input layer with certain adjacent polygons by erasing their common boundary. The adjacent polygon can be either the one with the largest or smallest area or the one sharing the largest common boundary with the polygon to be eliminated.

Eliminate is normally used to get rid of sliver polygons, i.e. tiny polygons that are a result of polygon intersection processes where boundaries of the inputs are similar but not identical.

It sounds like you'd want to try the "Largest Common Boundary" option.

  • I realize now you were asking for postgis solutions not qgis solutions. My apologies, I don't think postgis has an equivalent function but I'll leave this up for posterity. – FGreg Mar 20 at 22:49
0

This looks to me like a perfect use case for PostGIS topology extension. The topology's tolerance parameter will determine how far you allow vertices to snap to other existing polygons, to cope with the low precision of the source data and to clean it.

In short, the strategy is:

1. Enable the topology extension

CREATE EXTENSION postgis_topology;

2. Create a new empty topology

SELECT topology.CreateTopology('neighborhoods_topo', 4326, 1e-7);

The third parameter is the tolerance, in the units of the CRS; choose it wisely. Ideally, you want a CRS where unit is meters. If the CRS unit is not meters, as with WGS 84 aka 4326, use ST_Transform to reproject your polygons.

3. Add a TopoGeometry column to the polygons table

SELECT topology.AddTopoGeometryColumn('neighborhoods_topo', 'public', 'neighborhoods', 'topogeom', 'POLYGON');

This returns a new layer_id. Save it, it will be needed later. It will be layer 1 if your start from scratch, and incremented at every new call.

4. Add all polygons into the topology

UPDATE public.neighborhoods
SET topogeom = topology.toTopoGeom(geom, 'neighborhoods_topo', 1, 1e-7);

This can take several hours for a large dataset, be patient. 1 is the layer_id returned earlier.

5. Find faces appearing in several neighbourhoods

Find all faces from the topology that are present in 2 or more topogeometries. I will leave the query as an exercise. Easiest is probably with the GetTopoGeomElements function, then group by face id, and look at the ones with a count of 2 or more. Alternatively, you could create a new table with the cleaned geometry from the topogeom column, just cast it to standard geometry topogeom::geometry, and repeat what you already have now, but now with a clean dataset without the sliver overlaps.

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.