Cleaning geometries in PostGIS?

Question

I am trying to do some processing on some very large polygon layers. However I am running into various geometry errors such as:

NOTICE:  Ring Self-intersection at or near point 470396.52017068537 141300.52235257279
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 504154.61769969884 140782.04115761846
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 505255.50242871145 140803.34860398644
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 510312.46970004693 141215.29256710084
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 510312.46970004693 141215.29256710084
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 511839.50335641927 141115.85781738357
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 515064.03024010791 140895.68087158105
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 519233.18724611058 140881.47590733573
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 521072.73011588014 141044.83299615697
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 523331.31943088671 141144.26774587421
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 523331.31943088671 141144.26774587424
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 523395.24176999065 140725.22130063715
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1
NOTICE:  Ring Self-intersection at or near point 524531.63890961662 140810.45108610913
CONTEXT:  PL/pgSQL function st_intersection(geometry,raster,integer) line 10 at RETURN QUERY
SQL function "st_intersection" statement 1

I have tried the function suggested here: https://trac.osgeo.org/postgis/wiki/UsersWikiCleanPolygons

for cleaning geometries, the code I have used being:

UPDATE public.mytable
SET geom=cleangeometry(geom);

With the result:

ERROR:  GEOSisSimple: IllegalArgumentException: This method does not support GeometryCollection arguments

and also

UPDATE public.valid_mytable
SET geom=ST_MakeValid(geom);

This one works, but only if I first change my geometry column to geometry

ALTER TABLE public.mytable  ALTER COLUMN geom SET DATA TYPE geometry;

Which then leaves me with a table that no longer works with my other functions!

ERROR:  Relate Operation called with a LWGEOMCOLLECTION type.  This is unsupported.

I have tried changing the columns back to geometry(MultiPolygon)

ALTER TABLE public.my_table ALTER COLUMN geom SET DATA TYPE geometry(MultiPolygon);

But this fails

ERROR:  Geometry type (GeometryCollection) does not match column type (MultiPolygon)

I have tried going through PostGIS in Action (Second Ed) http://www.manning.com/obe/ but I can only find functions for finding invalid geometries, but my data-set is so big to fix this manually, I really need something that will fix them automatically.

---

I have been able to isolate the problem polygons, when I try and run ST_MakeValid() I get the result:

ERROR:  Geometry type (GeometryCollection) does not match column type      (MultiPolygon)
 ********** Error **********

 ERROR: Geometry type (GeometryCollection) does not match column type      (MultiPolygon)
SQL state: 22023

I did a type check on my geometry column, and it said the type was "MULTIPOLYGON"

Answer 1

You can try ST_CollectionExtract to extract [Multi]Polygons from GeometryCollections. Use ST_Multi to force them as MuliPolygons.

UPDATE public.valid_lcmsouthshapefile
  SET geom=ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3))
  WHERE NOT ST_IsValid(geom);

After your done, use a CHECK constraint to ensure they stay valid. See details here.

Answer 2

If you only want Polygons or Multipolygons from ST_MakeValid you can use ST_Dump to extract the constituent geometries and then test for the geometry type. ST_MakeValid will sometimes produce Points or LineStrings which is where the GeometryCollection is coming from. Try something like:

SELECT 
  g.geom, 
  row_number() over() AS gid,
FROM 
  (SELECT 
     (ST_DUMP(ST_MakeValid (geom))).geom FROM your_table
  ) AS g
WHERE ST_GeometryType(g.geom) = 'ST_MultiPolygon' 
   OR ST_GeometryType(g.geom) = 'ST_Polygon';

You could use an IN clause instead of OR condition, though the result and query plan would be the same. If you only want Multipolygons you can wrap ST_Dump in the ST_Multi function.

The row_number() over() will simply return you a unique id, starting from one, for each geometry returned from ST_Dump. You can also use the path element returned by ST_Dump, with the same result.

You will presumably want to combine this with a CREATE TABLE cleaned_geoms AS SELECT .... type statement, as a direct update is unlikely to work as ST_MakeValid wil notl generally (or always) produce a one to one mapping from I put to output.

This is untested as I have no means, currently, so there might be a misplaced parenthesis, but the general principle is sound. Hope this helps.