Display a polygonized raster in a Django Webapp (1/2): Polygonize a raster and save polygons to the database

This tutorial requires that a Django project has already been set-up with a PostGIS database. If you’re not sure about the steps needed to configure Django to serve geospatial data, this tutorial is a nice way to get started.

In this tutorial, we will see how to polygonize a raster image using rasterio to turn it into a vector. This vector can be subsequently either saved on the disk with Fiona or on the PostGIS database using Django (as will be done in this tutorial). We can then visualize the produced polygons one by one on the Django admin page.

Rasterio basically relies on GDAL and offers the same functionalities found in the python-GDAL wrapper, but with a more pythonic syntax, and hopefully fixes or at least avoids the common issues encountered with GDAL.

Rasterio can be installed in Linux with a pip:

However if on Windows, you need to install it with gdal from the wheel files built by Christoph Gohlke and available here.

The GeoTiff image we will polygonize is a subset of a binary flood map I’ve produced. It shows the town of Tewkesbury in South-West England which was flooded in 2007.

I’m going to assume in what follows that a Django app (let’s call it “maps”) has already been created and configured in “<project>/settings.py”. You can check out this tutorial if you’re not sure how to proceed. Inside the “maps” app, we will add the code below to “maps/models.py” that will result in the creation of a database table called “maps_flood” later:

The database table that will be generated after executing the next two shell commands will contain a single geometry column (besides the id obviously), corresponding to the PolygonField in the model class above:

Now we’ll move on to the “maps/views.py” file where the actual vectorization of the raster is carried out. We begin by importing the libraries needed:

Then, we create a view function that we’ll call “vectorize()” and which gets as an argument the http request. The code that we’ll explain below (reading the raster image & vectorizing it) will go inside this function:

We read the GeoTiff image by passing the path to it as an argument to rasterio.open(). The raster variable below contains the NumPy array image, while we keep copies of the projection in the crs variable and the geocoordinates of the upper-left corner of the image in the transform variable. The crs and the transform variables will prove useful later to get a georeferenced vector. We’ll actually need during the insertion in the database table the EPSG, which can be extracted from the crs dictionary in a quite straightforward way:

The polygonization of the raster comes next and deserves a thorough explanation. In the following for loop, we’re going through the pair of polygons, values returned by rasterio.features.shapes(). As mentioned above I’ve passed the transform variable to this function, in order to keep the same georeferencing of the input raster in the resulting vector. We also give as input to this function, besides the raster to polygonize, the same raster but as a boolean mask, so that the polygonization algorithm isn’t performed on the zero pixels in the raster.

Enough with the loop header. The for loop contains a single instruction which creates a polygon entry in the “maps_flood” database table in every iteration. We give as input to the create method the polygon object in the standard projection EPSG:4326 (the default one). The geometry argument accepts a GEOSGeometry object, which itself received as an input the polygon in the WKT format (hence the trick below):

We close the view function by returning an http response having as a content an h1 html tag, to show the user on the browser that the vectorization process is done:

Afterwards, just like in the previous tutorial we make the view created browsable by adding this code to “<project>/urls.py”:

And this other piece of code to “maps/urls.py”, to make the view function created above accessible from “localhost:8000/”:

If you browse to localhost:8000 the vectorization of the raster should be carried out and the obtained polygons should be saved in the PostGIS database. Note that this operations might take a few seconds, because polygonizing a raster and saving each polygon created to the database can be resource-intensive, especially if the number of polygons is important (429 resulting polygons in my case). The result can be visualized by going to the Django Admin section (localhost:8000/admin). Django-leaflet is needed to visualize the blue polygon on the map below on the admin section. As before this tutorial can come in handy to set it up:

In a future tutorial, I will show how we can display all the polygons on a Mapbox map in this same Django web app.

Don’t forget to clap if this story has been useful to you.