Creating Contour Lines on Folium Map with Python

Contour lines to add an extra dimension to you visualisations and maps

Photo by Annie Spratt on Unsplash

Credits:

• https://pypi.org/project/geojsoncontour/
• https://www.tjansson.dk/2018/10/contour-map-in-folium/
• Jorge Manchola for catching errors in the 4th block of the code

In many cases, we might encounter the need to put a third dimension such as elevation, rainfall rate, temperature, population and many others on the 2-D map. Creating contour lines is one way to do so. This article will guide you through step by step process in creating contour lines on Folium with Python.

Required libraries for this process are :

numpy, pandas, folium, branca, matplotlib.pyplot, scipy.interpolate, geojsoncontour, scipy, Raster2xyz, tiffile and rasterio. Some of them are those that commonly used, but the rest are special for mapping. Depending on what file format you have for your third dimension data, you might be able to skip some of the codes, or you might need to use different tools to read different file format. In this example, we will use a DEM (Digital Elevation Model) in TIFF format.

Reading DEM in TIFF and display it.

import tifffile as tiff
import matplotlib.pyplot as plt
import rasterio
filename = '../path/...TIF' # name the file path and name
tfile = tiff.imread(filename) # read the data
tfile.shape # arrange the data into an array
tiff.imshow(tfile) # display the map

The above codes results in the following figure. If this is the only thing you need, then you are all set.

Key West, Florida Map Produced from Digital Elevation Model

The figure above shows the map in x,y coordinates, and mainly shows elevations lower than mean sea level of zero. We might want to have more information out of the DEM, such as ground surface elevation. We might need to create contour lines to show more informative map. Let’s create contour lines and display it on Folium.

Contour Lines Generation

First, we would like to save our image in csv format as follows:

from raster2xyz.raster2xyz import Raster2xyz
input_raster = filename # set the input file name
out_csv = filename + ".csv" # set the output file name
rtxyz = Raster2xyz(input_raster) # transform raster to x,y,z
rtxyz.translate(input_raster, out_csv) # save to a csv file
DEM = pd.read_csv('Job504475_2007_FDEM_South_xyz.csv', header = None, names = ['longitude','latitude','elevation'])# read the x,y,z file and assign the column names
DEM.drop(DEM.head(1).index, inplace=True)# delete the first row that contain labels 'x,y,z'
DEM = DEM.replace(',','.', regex=True).astype(float) # set data type to float, they were saved as strings

In order to display our map on Folium, we need to transform our coordinates into longitude and latitude format. The loop below will do the job. However, keep in mind that depends on the file size, this process can be very time consuming. For a practice, use a small datasets to make sure that your codes work.

R = 6371 # Earth radius in km
for i in range(1,len(DEM)):
    DEM['latitude'][i] = np.degrees(np.arcsin(DEM['elevation'][i]/(3.280833*R*1000))) # our data is in ft.
    DEM['longitude'][i] = np.degrees(np.arctan2(DEM['latitude'][i]/(3.280833*1000), DEM['longitude'][i]/(3.280833*1000)))
DEM = DEM.to_csv("DEM.csv", index=True, header=True) # save to csv

Now we are ready to create the contour elevation of Key West Map. To do this, we need to import some libraries.You might have to pre-install some of them if you have not done so.

import folium
import branca
from folium import plugins
from scipy.interpolate import griddata
import geojsoncontour
import scipy as sp
import scipy.ndimage
%matplotlib inline
# Read the digital elevation  model in longitude, latitude and elevation format.
DEM = pd.read_csv('DEM.csv', header = None, names = ['longitude','latitude','elevation'])
# Set all negative elevations to zero to show that they will be under the mean sea level of zero
for i in range(0,len(DEM)):
    if DEM['elevation'][i] <= 0.:
        DEM['elevation'][i] = 0.
# Setup minimum and maximum values for the contour lines
vmin = DEM['elevation'].min() 
vmax = DEM['elevation'].max()
# Setup colormap
colors = ['blue','royalblue', 'navy','pink',  'mediumpurple',  'darkorchid',  'plum',  'm', 'mediumvioletred', 'palevioletred', 'crimson',
         'magenta','pink','red','yellow','orange', 'brown','green', 'darkgreen']
levels = len(colors)
cm     = branca.colormap.LinearColormap(colors, vmin=vmin, vmax=vmax).to_step(levels)
# Convertion from dataframe to array
x = np.asarray(DEM.longitude.tolist())
y = np.asarray(DEM.latitude.tolist())
z = np.asarray(DEM.elevation.tolist()) 
# Make a grid
x_arr          = np.linspace(np.min(x), np.max(x), 500)
y_arr          = np.linspace(np.min(y), np.max(y), 500)
x_mesh, y_mesh = np.meshgrid(x_arr, y_arr)
 
# Grid the elevation (Edited on March 30th, 2020)
z_mesh = griddata((x, y), z, (x_mesh, y_mesh), method='linear')
 
# Use Gaussian filter to smoothen the contour
sigma = [5, 5]
z_mesh = sp.ndimage.filters.gaussian_filter(z_mesh, sigma, mode='constant')
 
# Create the contour
contourf = plt.contourf(x_mesh, y_mesh, z_mesh, levels, alpha=0.5, colors=colors, linestyles='None', vmin=vmin, vmax=vmax)

The above code creates contour lines in MatPlotLib format. This is not compatible with Folium map. Therefore, we need to convert this contour into GeoJson format with the following code.

# Convert matplotlib contourf to geojson
geojson = geojsoncontour.contourf_to_geojson(
    contourf=contourf,
    min_angle_deg=3.0,
    ndigits=5,
    stroke_width=1,
    fill_opacity=0.1)
 
# Set up the map placeholdder
geomap1 = folium.Map([DEM.latitude.mean(), DEM.longitude.mean()], zoom_start=12, tiles="OpenStreetMap")
# Plot the contour on Folium map
folium.GeoJson(
    geojson,
    style_function=lambda x: {
        'color':     x['properties']['stroke'],
        'weight':    x['properties']['stroke-width'],
        'fillColor': x['properties']['fill'],
        'opacity':   0.5,
    }).add_to(geomap1)
 
# Add the colormap to the folium map for legend
cm.caption = 'Elevation'
geomap1.add_child(cm)
 
# Add the legend to the map
plugins.Fullscreen(position='topright', force_separate_button=True).add_to(geomap1)
geomap1

The following figure is our contour elevation map for Key West, Florida.

Contour Map of Key West, Florida

Enjoy creating maps, I hope I can share more mapping tricks with you.

Data resources: NOAA (National Oceanography and Atmospheric Agency)