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 0x7ff23fc31ed0>

Polygonize the lakes

polygonized_lakes = imod.prepare.polygonize(lake_grid)

polygonized_lakes.head(5)
value geometry
0 1.0 POLYGON ((94950.000 467550.000, 94950.000 4671...
1 1.0 POLYGON ((95150.000 466950.000, 95150.000 4667...
2 1.0 POLYGON ((103150.000 467550.000, 103150.000 46...
3 1.0 POLYGON ((109150.000 467550.000, 109950.000 46...
4 1.0 POLYGON ((102550.000 465950.000, 102950.000 46...


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.000 467550.000, 94950.000 4671...
1 1.0 POLYGON ((95150.000 466950.000, 95150.000 4667...
2 1.0 POLYGON ((103150.000 467550.000, 103150.000 46...
3 1.0 POLYGON ((109150.000 467550.000, 109950.000 46...
4 1.0 POLYGON ((102550.000 465950.000, 102950.000 46...


Plotted, we see a similar picture to the plotted raster

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

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
<AxesSubplot: >

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

Gallery generated by Sphinx-Gallery