imod.prepare.polygonize fails to handle island polygons correctly

[TeunvanWoerkom]

When using imod.prepare.polygonize, island polygons are not cut-out correctly from its 'sea' polygon. The polygons itself are correctly created, but no hole is made in its surrounding polygon.

imod version '1.0.0', rasterio version 1.4.4

Raster snippet:

Shape snippet:

[Huite]

I think the problem is the indexing here:

sg.Polygon(geom["coordinates"][0])

This seems a rather indefensible choice (... on my part, years ago). It should be straighforward to construct a Polygon using the first coordinate set as the exterior and the remainder as the interiors.

def polygonize(da):
    """
    Polygonize a 2D-DataArray into a GeoDataFrame of polygons.

    Parameters
    ----------
    da : xr.DataArray

    Returns
    -------
    polygonized : geopandas.GeoDataFrame
    """
    import geopandas as gpd
    import shapely.geometry as sg

    if da.dims != ("y", "x"):
        raise ValueError('Dimensions must be ("y", "x")')

    values = da.values
    if values.dtype == np.float64:
        values = values.astype(np.float32)

    transform = imod.util.transform(da)
    shapes = rasterio.features.shapes(values, transform=transform)

    geometries = []
    colvalues = []
    for geom, colval in shapes:
        geometries.append(sg.Polygon(geom["coordinates"][0]))
        colvalues.append(colval)

    gdf = gpd.GeoDataFrame({"value": colvalues, "geometry": geometries})
    gdf.crs = da.attrs.get("crs")
    return gdf

[JoerivanEngelen]

@Huite, your fix doesn't seem to change the code?

I guess you want to do something like:

...
        geometries.append(sg.Polygon(geom["coordinates"][0]), holes=geom["coordinates"][1:])
...

?

[TeunvanWoerkom]

Huite mentioned:
exterior = coordinates[0]
if len(coordinates) > 1:
interiors = coordinates[1:]
else:
interiors = None
polygon = sg.Polygon(exterior, interiors) and 'shapely noemt het shell en holes'. Not sure where shell and holes are located.

[JoerivanEngelen]

Thanks! Could you report if adapting above code with the proposed changes fixes your usecase? If so, I'll create a PR