Polygonize raster#

iMOD Python also provides convenience functions to polygonize rasters.

import matplotlib.pyplot as plt
import imod

We’ll start off by creating an example raster lake_grid to convert to polygons. This is similar to the Rasterize shapefiles example.

temp_dir = imod.util.temporary_directory()
lakes = imod.data.lakes_shp(temp_dir)

# Create dummy grid
xmin = 90950.0
xmax = 115650.0
dx = 200

ymin = 445850.0
ymax = 467550.0
dy = -200.0

like_2d = imod.util.empty_2d(dx, xmin, xmax, dy, ymin, ymax)

# Rasterrize the shapes
lake_grid = imod.prepare.rasterize(lakes, like=like_2d)

Our raster looks like this:

fig, ax = plt.subplots()
lake_grid.plot(ax=ax)
dx = 200.0, dy = -200.0
<matplotlib.collections.QuadMesh object at 0x00000289AFD4BF10>

Polygonize the lakes

polygonized_lakes = imod.prepare.polygonize(lake_grid)

polygonized_lakes.head(5)
value geometry
0 1.0 POLYGON ((94950 467550, 94950 467150, 95350 46...
1 1.0 POLYGON ((95150 466950, 95150 466750, 95350 46...
2 1.0 POLYGON ((103150 467550, 103150 466950, 103350...
3 1.0 POLYGON ((109150 467550, 109150 467350, 108750...
4 1.0 POLYGON ((102550 465950, 102550 465750, 102150...


This also polygonized the areas with np.nan. So we can drop those, using regular pandas functionality

polygonized_lakes = polygonized_lakes.dropna()

polygonized_lakes.head(5)
value geometry
0 1.0 POLYGON ((94950 467550, 94950 467150, 95350 46...
1 1.0 POLYGON ((95150 466950, 95150 466750, 95350 46...
2 1.0 POLYGON ((103150 467550, 103150 466950, 103350...
3 1.0 POLYGON ((109150 467550, 109150 467350, 108750...
4 1.0 POLYGON ((102550 465950, 102550 465750, 102150...


Plotted, we see a similar picture to the plotted raster

fig, ax = plt.subplots()
polygonized_lakes.plot(ax=ax)
polygonize raster
<Axes: >

Since it is a GeoDataFrame, we can now do vector operations. Like computing the centroids and plotting them as points.

centroids = polygonized_lakes.centroid

fig, ax = plt.subplots()
polygonized_lakes.plot(ax=ax)
centroids.plot(ax=ax, color="k")
polygonize raster
<Axes: >

Total running time of the script: (0 minutes 0.745 seconds)

Gallery generated by Sphinx-Gallery