How to convert a Raster to Contours with Python and GDAL - Tutorial

On the way to discover the analogy of the common desktop GIS procedures to Python we assumed that the process of extracting contours from raster was well documented or there were many tutorials on the topic. We found that there weren't many options to complete this process successfully or with few amount of pain, therefore we have done a complete tutorial on the process to create contours from a elevation raster with Python and GDAL that includes input data, scripts, and a discussion over the main steps of processing.

There are two ways to represent the resulting contour shapefiles, the first one is in QGIS and the second is inside Jupyter with Ipyleaflet.

You have three ways to run this code:

1. Using Anaconda with GDAL installed

2. From the Hakuchik server: hakuchik.hatarilabs.com (You need a Github account)

3. From the Hakuchik image for Docker: docker run -it -p 8888:8888 hatarilabs/hakuchik:74a2c9c18727

Tutorial

Input data

You download the input data from this link. Or you can just clone the repository with:

git clone https://github.com/SaulMontoya/rasterToContourwithPythonandGDAL.git

Python code

This is the complete Python code used for the tutorial

import numpy as np
from osgeo import osr
from osgeo import ogr
from osgeo import gdal

#Open tif file as select band
rasterDs = gdal.Open('../rst/20190109125130_1063922483_WGS13N.tif')
rasterBand = rasterDs.GetRasterBand(1)
proj = osr.SpatialReference(wkt=rasterDs.GetProjection())

#Get elevation as numpy array
elevArray = rasterBand.ReadAsArray()
print(elevArray[:4,:4])

#define not a number
demNan = -32768

#get dem max and min
demMax = elevArray.max()
demMin = elevArray[elevArray!=demNan].min()
print("Maximun dem elevation: %.2f, minimum dem elevation: %.2f"%(demMax,demMin))

[[-32768 -32768 -32768 -32768]
 [-32768 -32768 -32768 -32768]
 [-32768 -32768 -32768 -32768]
 [-32768 -32768 -32768 -32768]]
Maximun dem elevation: 2728.00, minimum dem elevation: 1247.00

contourPath = '../shp/contoursIncremental.shp'
contourDs = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(contourPath)

#define layer name and spatial 
contourShp = contourDs.CreateLayer('contour', proj)

#define fields of id and elev
fieldDef = ogr.FieldDefn("ID", ogr.OFTInteger)
contourShp.CreateField(fieldDef)
fieldDef = ogr.FieldDefn("elev", ogr.OFTReal)
contourShp.CreateField(fieldDef)

#Write shapefile using noDataValue
#ContourGenerate(Band srcBand, double contourInterval, double contourBase, int fixedLevelCount, int useNoData, double noDataValue, 
#                Layer dstLayer, int idField, int elevField
gdal.ContourGenerate(rasterBand, 50.0, 1250.0, [], 1, -32768., 
                     contourShp, 0, 1)

contourDs.Destroy()

contourPath = '../shp/contoursDefined.shp'
contourDs = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(contourPath)

#define layer name and spatial 
contourShp = contourDs.CreateLayer('contour', proj)

#define fields of id and elev
fieldDef = ogr.FieldDefn("ID", ogr.OFTInteger)
contourShp.CreateField(fieldDef)
fieldDef = ogr.FieldDefn("elev", ogr.OFTReal)
contourShp.CreateField(fieldDef)

#define number of contours and range
conNum = 50
conList =[int(x) for x in np.linspace(demMin,demMax,conNum)]

#Write shapefile using noDataValue
#ContourGenerate(Band srcBand, double contourInterval, double contourBase, int fixedLevelCount, int useNoData, double noDataValue, 
#                Layer dstLayer, int idField, int elevField
gdal.ContourGenerate(rasterBand, 0, 0, conList, 1, -32768., 
                     contourShp, 0, 1)

contourDs.Destroy()

from ipyleaflet import Map, GeoData, LayersControl
import geopandas as gpd

m = Map(center=(27.8840, -107.8895), zoom=10)

contourDf = gpd.read_file('../shp/contoursDefined.shp')
contourDfWgs84 = contourDf.to_crs(epsg=4326)

geo_data = GeoData(geo_dataframe = contourDfWgs84 )

m.add_layer(geo_data)
m.add_control(LayersControl())

m

Map(center=[27.884, -107.8895], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'z…