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
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
