[gdal-dev] Strange (expected?) behaviour exporting multipart geometry to ESRI Shapefile format

[Even Rouault]

Andrea,

your [2] link contains the shapefile, not the GeoPackage file.

That said I've tried replicated using 
https://github.com/qgis/QGIS/files/8058606/Issue47288_Polygons_GPKG.zip 
from the QGIS ticket and converting it to shapefile. The issue here is 
that the original geometry is not (considered as) valid:

$ ogrinfo Issue47288_Polygons.gpkg -sql "select st_isvalid(geom) from 
Issue47288_Polygons" -al -q
GEOS warning: Self-intersection at or near point 352107.37400937825 
5662263.0605139844 54.963999999999999

Layer name: SELECT
OGRFeature(SELECT):0
   st_isvalid(geom) (Integer) = 0

When writing polygons to shapefile, the OGR shapefile driver takes all 
the rings of a (multi)polygons and detect which ones are inner rings of 
which outer rings, to apply the expected winding order of the shapefile 
specification. That analysis assumes that the geometries are valid, 
which simplifies (=speed up) the topological analysis. If they aren't, 
then the algorithm might produce "unexpected" results (but there's no 
real expected result when the geometry is not valid)

One possibility is to force the validity during the conversion, but this 
will cause one of the parts to be discarded:

ogr2ogr out2.shp Issue47288_Polygons.gpkg -sql  "select 
ExtractMultiPolygon(st_makevalid(geom)) as geom, * from 
Issue47288_Polygons" -dialect indirect_sqlite

Or you can use the SHAPE_REWIND_ON_WRITE configuration option mentioned 
at https://gdal.org/drivers/vector/shapefile.html#advanced-topics :

ogr2ogr out3.shp Issue47288_Polygons.gpkg --config SHAPE_REWIND_ON_WRITE OFF

This will keep the original geometry mostly untouched, possibly creating 
a non-conformant shapefile. And on the reading side, the OGR shapefile 
reader might also have issues reconstructing a single-feature geometry. 
On that particular case, it seems to produce the "expected" result

I think the gist of the issue is that we have here 3D geometries and 
that the algorithms in the shapefile driver or in GEOS that check 
validity ignore the Z dimension. Reading the OGC  Simple feature access 
- Part 1: Common architecture  spec 
(https://portal.ogc.org/files/?artifact_id=25355) , it is not entirely 
clear to me what the validity rules are for a multipolygon whose 
polygons are not co-planar

In § 6.1.10.1 (Surface description) it i said that "A single such 
Surface patch in 3-dimensional space is isometric to planar Surfaces"

As far as I can see in this multipolygon, there are 2 issues: one of the 
part is repeated, and the 2 "small" rectangular polygons that share the 
352107.37400937825 5662263.0605139844 54.963999999999999 points are 
likely non-planar.

In § 6.1.13.1 (MultiSurface description) it is also said that "The 
boundaries of any two coplanar elements in a MultiSurface may intersect, 
at most, at a finite number of Points", which seems to imply that the 
parts of a MultiSurface might be non coplanar.

AFAICS, the shapefile specification is silent on that aspect of validity 
rules regarding 3D geometries.

My inclination would be to consider that in the case where we have a 
multipolygon with non -coplanar parts we should probably be avoiding 
considering the ones that are non-coplanar are inner/outer rings of others.

I've ticketed that in https://github.com/OSGeo/gdal/issues/5315

Even

Le 16/02/2022 à 08:07, Andrea Giudiceandrea via gdal-dev a écrit :
> Hi GDAL devs,
> trying find the root cause of a bug [1] reported in the QGIS GitHub 
> repository, I've noticed a strange behaviour when a multipart geometry 
> is exported to an ESRI Shapefile layer.
>
> I've created a minimal GeoPackage layer file [2] which contains 1 
> multipart polygon geometry (MultiPolygon) feature consisting of 3 parts.
> The 3 parts are 3 clockwise polygons external rings.
>
> The following command for such GeoPackage layer:
> ogrinfo Polygon_3.gpkg Polygon_3 -geom=SUMMARY
>
> reports:
> OGRFeature(Polygon_3):0
>   MULTIPOLYGON : 3 geometries:
> POLYGON : 7 points
> POLYGON : 9 points
> POLYGON : 21 points
>
> It seems to me something strange happens converting the GeoPackage 
> layer to an ESRI Shapefile layer using e.g:
> ogr2ogr -f "ESRI Shapefile" Polygon_3.shp Polygon_3.gpkg
>
> ogrinfo for the resulting ESRI Shapefile layer reports:
> OGRFeature(Polygon_3):0
>   FID (Integer64) = 0
>   MULTIPOLYGON : 2 geometries:
> POLYGON : 9 points
> POLYGON : 21 points, 1 inner rings (7 points)
>
> That is, the 7 points part has been written in the Shapefile layer as 
> a counter clockwise polygon and "attached" to the 21 points part 
> (which remains a clockwise polygon external ring) as his inner ring. 
> The 9 points part remains a clockwise polygon external ring.
>
> On the contrary, converting the GeoPackage layer to GML or FlatGeobuf 
> or GeoJSON or Spatialite formats, the resulting layer correctly 
> contains 1 MultiPolygon feature consisting of 3 clockwise Polygon 
> parts like the original layer and ogrinfo reports:
> OGRFeature(Polygon_3):0
>   MULTIPOLYGON : 3 geometries:
> POLYGON : 7 points
> POLYGON : 9 points
> POLYGON : 21 points
>
> I've checked the orientation of the polygons parts using QGIS 3.22.3 
> for all the layers and also ArcMap/ArcGIS 9.3.1 for the Shapefile layer.
>
> QGIS reports 3 parts for the geometry in the GeoPackage, GML, 
> FlatGeobuf, GeoJSON or Spatialite layers, while it reports 2 parts for 
> the ESRI Shapefile layer.
>
> Is this an expected behaviour of the ESRI Shapefile driver writing 
> multipart features geometries?
>
> Best regards.
>
> Andrea Giudiceandrea
>
> [1] https://github.com/qgis/QGIS/issues/47288
> [2] https://drive.google.com/file/d/10HiCARTYJqhAXUckK43nU4MzV1wUbLg9
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.