Add a method to compute overlaps between grid cells and polygons

[Huite]

We can triangulate the polygons using earcut, build a grid, and use an OverlapRegridder to get statistics:

This is the basic idea:

vertices, index = shapely.get_coordinates(gdf.geometry.to_numpy(), return_index=True)

split_at = np.flatnonzero(np.diff(index)) + 1
polygon_verts = np.split(vertices, split_at)

result = []
increment = 0
for i, ring in enumerate(polygon_verts):
    print(i)
    connectivity = earcut.triangulate_float64(ring, [len(ring)])
    result.append(connectivity + increment)    
    increment += len(ring)
out = np.concatenate(result).reshape((-1, 3))
grid = xu.Ugrid2d(*vertices.T, -1, out)

xmin, ymin, xmax, ymax = grid.bounds
lhm_grid = imod.util.empty_2d(xmin=xmin, xmax=xmax, dx=250.0, ymin=ymin, ymax=ymax, dy=-250.0)
regridder = xu.OverlapRegridder(source=grid, target=lhm_grid)

# %%

ds = regridder.to_dataset()
i = ds["__regrid_indptr"]
j = ds["__regrid_indices"]
v = ds["__regrid_data"].to_numpy()

from scipy import sparse

matrix = sparse.csr_matrix((v, j, i), shape=(ds["__regrid_n"].item(), ds["__regrid_m"].item()))
total_area = np.asarray(matrix.sum(axis=1)).ravel()
out = lhm_grid.copy(data=total_area.reshape(lhm_grid.shape)[::-1, :])

It appears quite performant; I did run into two easily fixable issues: https://github.com/Deltares/xugrid/issues/233 https://github.com/Deltares/numba_celltree/issues/13