shapefile writer: issue when writing non-coplanar 3D multipolygons

[rouault]

From https://lists.osgeo.org/pipermail/gdal-dev/2022-February/055443.html

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.

[agiudiceandrea]

Hi @rouault,
I've found a similar issue that does occur also using GDAL 3.6 with a very simple geometry.

For the MultiPolygonZ GeoPackage layer multipz.gpkg (in multipz_gpkg.zip)

ogrinfo multipz.gpkg multipz -geom=SUMMARY reports:

Feature Count: 1
...
  MULTIPOLYGON : 1 geometries:
POLYGON : 4 points, 1 inner rings (4 points)

ogrinfo multipz.gpkg -sql "select st_isvalid(geom) from multipz" -al -q reports:

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

The outer (first) ring (biggest triangle) is CCW and the inner (second) ring (smallest triangle) is CW.

QGIS displays such layer as a triangular polygon with a triangular hole:

The Cartesian area calculated for such geometry (5243 m²) is equivalent to the area enclosed by the outer ring minus the area enclosed by the inner ring.

The QGIS "Check validity" (GEOS) tool reports the geometry as valid.

---

Converting the GeoPackage layer to an ESRI Shapefile layer using ogr2ogr multipz.shp multipz.gpkg or QGIS, the layer multipz.shp (in multipz_shp.zip) is created.

For such layer
ogrinfo multipz.shp multipz -geom=SUMMARY reports

Feature Count: 1
...
  MULTIPOLYGON : 2 geometries:
POLYGON : 4 points
POLYGON : 4 points

Both the outer (first) ring (biggest triangle) and the inner (second) ring (smallest triangle) are CW.

QGIS displays such layer as a two overlapping triangular polygons parts:

The Cartesian area calculated for such geometry (8267 m²) is equivalent to the area enclosed by the outer ring plus the area enclosed by the inner ring.

The QGIS "Check validity" (GEOS) tool reports the geometry as invalid ("Nested shells").

---

Converting the GeoPackage layer to an ESRI Shapefile layer using the command ogr2ogr multipz2.shp multipz.gpkg -sql "select ExtractMultiPolygon(st_makevalid(geom)) as geom, * from multipz" -dialect ind creates the same ESRI Shapefile layer as before and thus doesn't fix the issue.

---

Converting the GeoPackage layer to an ESRI Shapefile layer using the command ogr2ogr multipz3.shp multipz.gpkg --config SHAPE_REWIND_ON_WRITE OFF, the layer multipz3.shp (in multipz3_shp.zip) is created.

For such layer
ogrinfo multipz3.shp multipz3 -geom=SUMMARY reports

Feature Count: 1
...
  POLYGON : 4 points, 1 inner rings (4 points)

The inner (first) ring (smallest triangle) is CW and the outer (second) ring (biggest triangle) is CCW.

QGIS displays such layer as a triangular polygon with a triangular hole:

The Cartesian area calculated for such geometry is negative (-5243 m²).

The QGIS "Check validity" (GEOS) tool reports the geometry as invalid ("Hole lies outside shell").

[jratike80]

The test triangle as WKT

MULTIPOLYGON Z(((516113.631069 5041435.137874 137.334, 516141.2239 5041542.465874 137.614, 515998.390418 5041476.527121 137.288, 516113.631069 5041435.137874 137.334), (516041.808551 5041476.527121 137.418, 516111.602184 5041505.337284 137.322, 516098.617322 5041456.644051 137.451, 516041.808551 5041476.527121 137.418)))

[rouault]

I've done the maths with the above geometry, and the "inner ring" is not coplanar with the "outer ring". So as far as I can see, GDAL is doing the right thing here considering this is actually not a inner ring
You cannot trust GEOS based utilities regarding 3D geometries as GEOS ignores the 3D dimension.

[rouault]

hum just thinking that the Shapefile writer should actually not call the fragile logic of SHPRewindObject() but analyze the structure of the OGRPolygon to get which rings are outer & inner rings (even if they don't really make sense with respect to the Simple Features criteria like here)

[jratike80]

I've done the maths with the above geometry, and the "inner ring" is not coplanar with the "outer ring". So as far as I can see, GDAL is doing the right thing here considering this is actually not a inner ring

Hmm, I think that lots of 2.5D polygons I have seen in GIS data are not coplanar. The z-coordinate of vertices of cadastral parcels etc. present usually the height at each point and polygons are like rubber sheets spread over the terrain. The outer ring can have whatever Z values and be anything but coplanar and GDAL accepts it. I think that the same lenient rule should apply to inner rings as well.

[rouault]

Fix in #7157