pyogrio icon indicating copy to clipboard operation
pyogrio copied to clipboard

Bug: write_dataframe does not properly handle Z geometry types

Open brendan-ward opened this issue 2 years ago • 2 comments

For example:

df = geopandas.GeoDataFrame(geometry=[pygeos.points([0,0,0])])
write_dataframe(df, '/tmp/test.fgb')

fails because FGB is picky about mismatched geometry types and our geometry type inference resolves this to Point instead of Point Z, but the WKB data passed to GDAL contain Z values. This is because the GeoSeries reports this as Point:

>>> df.geometry.type.unique()
array(['Point'], dtype=object)

Which happens because pygeos reports this as Point type (Z values are not given specific type codes):

>>> pygeos.get_type_id(df.geometry.values.data)
array([0], dtype=int32)

(POINT is 0)

Shapely does the same:

>>> df.geometry.values[0].type
'Point'

I think the upshot is that we also have to intersect the geometry types we obtain from the GeoSeries with " Z" if the geometry has Z values (@jorisvandenbossche is this something that should be handled in GeoPandas instead, so that the df.geometry.type reports Point Z instead of Point?)

Specifying the geometry type directly:

write_dataframe(geopandas.GeoDataFrame(geometry=[pygeos.points([0,0,0])]), '/tmp/test.fgb', layer_geometry_type="Point Z")

fails, because we don't have Point Z in _geometry.pyx::GEOMETRY_TYPE_CODES; we only have the 2.5D variants (which is technically how GDAL would store them, but we should consider using "<Geometry> Z" aliases to them because I think that all GEOS geometry Z types are equivalent to GDAL 2.5D geometry types).

brendan-ward avatar Jun 07 '22 12:06 brendan-ward

The latter case is resolved in 8679f3e2a07b62940b8bd9f55d1b4c690acc76ef (whoops, meant to make that a PR)

brendan-ward avatar Jun 07 '22 16:06 brendan-ward

is this something that should be handled in GeoPandas instead, so that the df.geometry.type reports Point Z instead of Point?)

We can't just change df.geometry.type, I think, but it might make sense to add some helper in geopandas to get this (combining the output of df.geometry.type with df.geometry.has_z). We should probably start with just adding this in pyogrio (as we can't rely on it being present of geopandas on the short term anyway)

jorisvandenbossche avatar Jun 07 '22 18:06 jorisvandenbossche