API: how to handle Polygon vs MultiPolygon in Shapefile reader

[jorisvandenbossche]

Moving this from a TST test update PR in #665 to an actual issue for a bit more visibility.

Summary of the upstream GDAL change (OSGeo/gdal#14602 / OSGeo/gdal#14662): when reading a Shapefile with a geometry field with a shape type of polylines or polygons, GDAL historically reported (when querying the geometry type of the field) those as LineString and Polygon, respectively, even though for the Shapefile specification those types can hold both LineStrings and MultiLineStrings or Polygons and MultiPolygons. When reading the actual geometry, GDAL would return LineString/Polygon or MultiLineString/MultiPolygon depending on the number of parts of the individual geometry (thus, when reading multiple rows, you can end up with a mixture of single and multi types).
With the abovementioned PRs, GDAL 3.14+ will now always report those fields as MultiLineString and MultiPolygon, regardless of actual geometries, since those types can encompass both single and multi geoms, and on reading also return all geometries as the Multi variant (but provides an PROMOTE_TO_MULTI=NO option to preserve the old behaviour).

While this change certainly makes sense for GDAL (e.g. for translating to a different format), for us where the main use case is reading the data into memory, we might want to keep the mix of single and multi types (e.g. geopandas/shapely can handle the mixture just fine).

Because the new behaviour means that you will now always get only MultiPolygons, even when storing only polygons. In comparison, for e.g. GPKG we default to "promote to multi" when writing a GPKG file if there is a mixture if single and multi polygons, and so in that case you will also read it as all MultiPolygons, but a dataset with all single Polygons will still roundtrip as such.

---

Small illustration: roundtripping (writing and reading again) a small dataframe with either two polygons, two multipolygons, or one polygon and one multipolygon (mixed). The table below then shows the geometry types of the roundtripped geometries per file format. The GeoJSON case preserves the geometry types exactly, so that columns illustrates the original data:

Details

import pyogrio
import shapely
import geopandas
import pandas as pd

geom = shapely.box(0, 0, 1, 1)
poly = shapely.box(0, 0, 1, 1)
poly1 = shapely.box(0, 0, 1, 1)
poly2 = shapely.box(2, 2, 3, 3)
mpoly = shapely.MultiPolygon([poly1, poly2])

gdf_single = geopandas.GeoDataFrame({"col": [0, 1]}, geometry=[poly1, poly2], crs=4326)
gdf_multi = geopandas.GeoDataFrame({"col": [0, 1]}, geometry=[mpoly, mpoly], crs=4326)
gdf_mixed = geopandas.GeoDataFrame({"col": [0, 1]}, geometry=[poly1, mpoly], crs=4326)

for ext in ["geojson", "shp", "gpkg"]:
    pyogrio.write_dataframe(gdf_single, f"test-poly-single.{ext}")
    pyogrio.write_dataframe(gdf_multi, f"test-poly-multi.{ext}")
    pyogrio.write_dataframe(gdf_mixed, f"test-poly-mixed.{ext}")

# roundtrip and check geometry types
result = defaultdict(dict)
for ext in ["geojson", "shp", "gpkg"]:
    for geom in ["single", "multi", "mixed"]:
        result[ext][geom] = pyogrio.read_dataframe(f"test-poly-{geom}.{ext}").geometry.geom_type.tolist()

print(pd.DataFrame(result).to_markdown())

	geojson	shp	gpkg
single	['Polygon', 'Polygon']	['Polygon', 'Polygon'] (-> ['MultiPolygon', 'MultiPolygon'])	['Polygon', 'Polygon']
multi	['MultiPolygon', 'MultiPolygon']	['MultiPolygon', 'MultiPolygon']	['MultiPolygon', 'MultiPolygon']
mixed	['Polygon', 'MultiPolygon']	['Polygon', 'MultiPolygon'] (-> ['MultiPolygon', 'MultiPolygon'])	['MultiPolygon', 'MultiPolygon']

The above is with released GDAL, with the GDAL 3.14 change highlighted in bold. So the change for Shapefile would mean that the full column would become ['MultiPolygon', 'MultiPolygon']. For the mixed case this is certainly fine I think (since we get that for GPKG by default as well). But so also the first row for the all-single-Polygons case would become MultiPolygons, while we do "preserve" the geometry type here for all other file formats.

Thoughts or preferences? Preserve the current reading behaviour for Shapefiles, or bite the bullet and align with GDALs change? (in shapely, a single vs multi Geometry has slightly different APIs, so this can cause some breakage)

cc @brendan-ward @martinfleis @ljwolf

[jorisvandenbossche]

There is also already a bit of discussion between @theroggy and me in #665 (comment). If we would deviate from GDALs default, there are two places to decide about this: the geometry type reported in metadata, and the actual geometry read (the WKB type).

For the geometry type reported in metadata, there we could align with GDAL to always say "MultiLineString" and "MultiPolygon", since we don't know the exact type before reading the full data (unless we would use a different API to inspect the actual geometry types, but that would make this slower for Shapefiles).
And we could then do this also for older GDAL versions as well, so this is a change in a new pyogrio version regardless of the GDAL version used.

---

I think I would then lean towards indeed reporting the multi-type (i.e. new GDAL behaviour), but when reading the actual data to preserve the single/multi distinction (i.e. old GDAL behaviour), but not a very strong opinion.

[martinfleis]

This is a bit pain. I can foresee a lot of breakage when we start reading stuff as Multi types given those do not have stuff like exterior and so on. But I am also not super keen on deviating from GDAL.

I'd be keen to go with your suggestion to follow GDAL in metadata but stick to old behaviour in geometry reading.

PS: http://switchfromshapefile.org

[theroggy]

As already mentioned in the PR, I would also go for reporting multi-type and reading the data without promoting everything to multi...