Visualize DEM in An Interactive Map

Generally, a raster layer could easily be overlaid on an interactive map colorfully using colormaps from matplotlib.cm through the class of folium.raster_layers.ImageOverlay, when there are not nodata values. However, when the raster data come with nodata values, the overlaid layer would look like a little bit patchy. Sometimes, the nodata value is set to a very big negative value such as -32767. Although we can replace them with numpy.nan, the overlaid layer would look still weird.

This tutorial walks you through how to visualize a DEM data on an interactive map created with the python package of folium. At the same time, we will deal with the no-data issue using a self-defined function to map DEM values to colors and filter the nodata values.

The default coordinate system and projection for web-based basemaps is WGS84 Web Mercator. To overlay data (here DEM) on web-based basemaps, the overlay data needs to be in the WGS84 coordinate system. Refer to previous tutorials to learn how to reproject a DEM data.

• folium

folium builds on the data wrangling strengths of the Python ecosystem and the mapping strengths of the leaflet.js library. Manipulate your data in Python, then visualize it in on a Leaflet map via folium.

%matplotlib inline

import numpy as np
import xarray as xr
import numpy.ma as ma
import folium

from IPython.display import IFrame
import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings("ignore")

Read and preprocess data

The demo data must contain some grids with no-data, which will show how we deal with raster with NAN values.

• Our demo data has been reprojected to EPSG:4326.
• Set them to np.nan.

infile  = "data/es_dem/pre_DTM_EPSG4326.tif"
da_dem  = xr.open_rasterio(infile).drop('band')[0].rename({'x':'longitude', 'y':'latitude'})
nodata  = da_dem.nodatavals[0]
da_dem  = da_dem.where(da_dem>nodata, np.nan)
arr_dem = da_dem.values

Calucate center and bounds - These will be used for folium.Map

clat = da_dem.latitude.values.mean()
clon = da_dem.longitude.values.mean()

mlat = da_dem.latitude.values.min()
mlon = da_dem.longitude.values.min()

xlat = da_dem.latitude.values.max()
xlon = da_dem.longitude.values.max()

print(clat, clon, mlat, mlon, xlat, xlon)

Map values to colors

function for mapping values to colors

def colorize(array, cmap='viridis'):
    normed_data = (array - array.min()) / (array.max() - array.min())    
    cm = plt.cm.get_cmap(cmap)    
    return cm(normed_data)  

Tip - The key point is that DEM values must be masked. So the grids with no-data would be set totally transprant (i.e., alpha=0).

data         = ma.masked_invalid(arr_dem)
colored_data = colorize(data, cmap='viridis_r')

Visualize interactively

It is a little bit weried that the DEM image could not be displayed in folium directly. Perhaps, it is too big? However, we can save the map to a HTML file and then display it in a IPython.display.IFrame.

mapa = folium.Map(location=[clat, clon], 
                  tiles="Stamen Terrain",
                  zoom_start=14) 

folium.raster_layers.ImageOverlay(colored_data,
                                  [[mlat, mlon], [xlat, xlon]],
                                  opacity=0.7).add_to(mapa)

html_file = 'index.html'
mapa.save(html_file)
IFrame(src=html_file, width=900, height=400)

References

Travis E, Oliphant. A guide to NumPy, USA: Trelgol Publishing, (2006).

Stéfan van der Walt, S. Chris Colbert and Ga?l Varoquaux. The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science & Engineering, 13, 22-30 (2011), DOI:10.1109/MCSE.2011.37

Fernando Pérez and Brian E. Granger. IPython: A System for Interactive Scientific Computing, Computing in Science & Engineering, 9, 21-29 (2007), DOI:10.1109/MCSE.2007.53

John D. Hunter. Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 9, 90-95 (2007), DOI:10.1109/MCSE.2007.55

https://www.dhirubhai.net/pulse/overlay-raster-nodata-interactive-map-folium-chonghua-yin/

https://python-visualization.github.io/folium/modules.html#module-folium.raster_layers

https://www.earthdatascience.org/courses/earth-analytics-python/lidar-raster-data/interactive-maps/

https://github.com/python-visualization/folium/tree/master/examples?1560294699938

https://rasterio.readthedocs.io/en/stable/index.html

https://www.dhirubhai.net/pulse/postprocess-swat-simulations-5-fancy-interactive-map-chonghua-yin/?lipi=urn%3Ali%3Apage%3Ad_flagship3_profile_view_base_post_details%3BOD%2F3r0YESPOZwq%2FgVg03sQ%3D%3D

https://qingkaikong.blogspot.com/2016/06/using-folium-5-image-overlay-overlay.html