Ggtkuk

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.

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.

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

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.

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.

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.

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.