[raster] st_dumpaspolygons makes not valid polygons

[jomarlla]

Dear PostGIS team,

ST_Dumpaspolygons is building not valid polygons.

select count(*) from (select st_dumpaspolygons (rast) as vect from co1) as tabla where not st_isvalid ((vect).geom);

NOTICE: Ring Self-intersection at or near point 503582.96553112275 4420706.2469141362 NOTICE: Ring Self-intersection at or near point 505882.96553112275 4417206.2469141362 ...282 more errors.

Some cases:: POLYGON((462682.965531123 4376006.24691414,462682.965531123 4375706.24691414,462982.965531123 4375706.24691414,462982.965531123 4375606.24691414,463182.965531123 4375606.24691414,463182.965531123 4375706.24691414,463082.965531123 4375706.24691414,463082.965531123 4375806.24691414,463182.965531123 4375806.24691414,463182.965531123 4375706.24691414,463282.965531123 4375706.24691414,463282.965531123 4375806.24691414,463382.965531123 4375806.24691414,463482.965531123 4375806.24691414,463482.965531123 4375906.24691414,463282.965531123 4375906.24691414,463282.965531123 4376006.24691414,462682.965531123 4376006.24691414))

POLYGON((500282.965531123 4380006.24691414,500282.965531123 4379606.24691414,500382.965531123 4379606.24691414,500382.965531123 4379406.24691414,500582.965531123 4379406.24691414,500582.965531123 4379306.24691414,500682.965531123 4379306.24691414,500682.965531123 4379206.24691414,500782.965531123 4379206.24691414,500782.965531123 4379506.24691414,500582.965531123 4379506.24691414,500582.965531123 4379606.24691414,500482.965531123 4379606.24691414,500482.965531123 4379706.24691414,500382.965531123 4379706.24691414,500382.965531123 4379806.24691414,500482.965531123 4379806.24691414,500482.965531123 4379706.24691414,500582.965531123 4379706.24691414,500582.965531123 4379906.24691414,500382.965531123 4379906.24691414,500382.965531123 4380006.24691414,500282.965531123 4380006.24691414))

POLYGON((516082.965531123 4364206.24691414,516082.965531123 4363806.24691414,516382.965531123 4363806.24691414,516382.965531123 4363906.24691414,516282.965531123 4363906.24691414,516282.965531123 4364006.24691414,516182.965531123 4364006.24691414,516182.965531123 4364106.24691414,516282.965531123 4364106.24691414,516282.965531123 4364006.24691414,516382.965531123 4364006.24691414,516382.965531123 4363906.24691414,516582.965531123 4363906.24691414,516582.965531123 4364006.24691414,516682.965531123 4364006.24691414,516782.965531123 4364006.24691414,516782.965531123 4364206.24691414,516582.965531123 4364206.24691414,516582.965531123 4364106.24691414,516482.965531123 4364106.24691414,516482.965531123 4364206.24691414,516082.965531123 4364206.24691414))

2 notes:

ST_Makevalid fixes all of them and I think in a right way.

Any of the fixed polygons become a multipolygon after fixing. Regards,

[pracine]

We could simply add a ST_MakeValid() call to the plpgsql st_dumpaspolygons function (at a little performance cost). Users can also easily do it with the geom part of the geomval.

any good reason why we would like to do it at the source?

What did the fact that the returned geometry were invalid prevented you to do?

[jomarlla]

ST_Makevalid is making st_dumpaspolygons much slower. Dont know if it is worth it to do it at the source. I was wondering if the algorithm can be changed for not making invalid polygons instead of using st_makevalid without a lot of work.

I just dont like invalid polygons because they just bother me about any future spatial analysis about unpredictable results. As you said users can do it easily but in that case maybe the docs should point out that st_dumpaspolygons can return invalid geometries, then the user can take care about that.

Jose

[pracine]

Replying to jomarlla:

ST_Makevalid is making st_dumpaspolygons much slower. Dont know if it is worth it to do it at the source. I was wondering if the algorithm can be changed for not making invalid polygons instead of using st_makevalid without a lot of work.

I guess the amount of work for producing a valid geometry from what GDAL is returning us would be equivalent to ST_MakeValid but I might be wrong.

I just dont like invalid polygons because they just bother me about any future spatial analysis about unpredictable results. As you said users can do it easily but in that case maybe the docs should point out that st_dumpaspolygons can return invalid geometries, then the user can take care about that.

If there is no simple low cost way to make them valid we will add a note. It's the first time someone is complaining about that so I guess the type of invalidity does not affect most analyses. Would be great if you could find an example of invalid post st_dumpaspolygon analysis.

[Bborie Park]

After doing some thinking about GDALPolygonize() and GDALFPolygonize() (both used by ST_DumpAsPolygons), it is possible/probably for invalid polygons to be generated. From the GDAL API docs for GDALPolygonize:

This function creates vector polygons for all connected regions of pixels in the raster sharing a common pixel value

Two pixels with the same value with one common corner may cause a bow-tie...

   +--+
   |  |
   |  |
+--+--+
|  |
|  |
+--+

I don't see any way to make GDAL generate valid polygons but we could add a validation routine to the appropriate C code in rt_api.c.

[pracine]

There is already a ticket (#637) to convert wkt geometries to wkb geometries. We could probably call some the right GEOS function to make it valid afterward. So first make it a true geometry second make it a valid one all in C. That should reduce the performance impact. I would set this one to 2.0.1.

[Bborie Park]

So for a simple bow-tie situation, ST_DumpAsPolygons returns two polygons

WITH foo AS (
	SELECT
		ST_SetValue(
			ST_SetValue(
				ST_AddBand(
					ST_MakeEmptyRaster(2, 2, 0, 0, 1, 1, 0, 0, 0)
					, '8BUI'::text, 0, 0
				)
				, 1, 1, 1
			)
			, 2, 2, 1
		) AS rast
)
SELECT
	ST_AsText((ST_DumpAsPolygons(rast)).geom)
FROM foo

Need to find an example of a raster dumping invalid polygons.

[Bborie Park]

Success! Found one...

WITH foo AS (
	SELECT (ST_DumpAsPolygons(
		ST_SetValue(
			ST_SetValue(
				ST_AddBand(
					ST_MakeEmptyRaster(4, 4, 0, 0, 1, 1, 0, 0, 0)
					, '8BUI'::text, 1, 0
				)
				, 2, 2, 0
			)
			, 3, 3, 0
		)
	)).geom AS geom
)
SELECT
	ST_IsValid(geom),
	ST_AsText(geom)
FROM foo

[Bborie Park]

Fixed in r9280