In this tutorial, we create a Relative Elevation Model (REM) in Python using xarray.
A REM is created by subtracting the river's elevation from a Digital Elevation Model, displaying the height above river. This makes the river's meanders more visible.
# Install Git LFS
!curl -s https://packagecloud.io/install/repositories/github/git-lfs/script.deb.sh | sudo bash
!apt-get install git-lfs
!git lfs install
# Clone repo
!rm REM-xarray -rf
!git clone https://github.com/DahnJ/REM-xarray.git
# Change working directory and install requirements
%cd REM-xarray
!pip install -r requirements.txt
from pathlib import Path
from IPython.core.display import Video
import numpy as np
import pandas as pd
import geopandas as gpd # Vector data handling
import osmnx as ox # Downloading data from OSM
from shapely.geometry import box
from scipy.spatial import cKDTree as KDTree # For Inverse Distance Weight calculation
import xarray as xr
import xrspatial # Hillshading
import rioxarray # Working with geospatial data in xarray
import matplotlib.pyplot as plt
from datashader.transfer_functions import shade, stack
For this tutorial, we will use DEM of the Ivalo river from Finnish NLS.
Examples of other good sources:
Video("data/howto-download-ivalo-data.mp4", embed=True)
The downloaded data can be combined using the convenience functioncombine_by_coords
.
Here we use a prepared netCDF file.
path = Path('data').joinpath('ivalo', 'ivalo.nc')
dem = rioxarray.open_rasterio(path)
We will make the DEM smaller so that it's easier to work with
dem = dem.coarsen(x=3, boundary='trim').mean().coarsen(y=3, boundary='trim').mean()
Let's plot the DEM
dem.squeeze().plot.imshow()
<matplotlib.image.AxesImage at 0x7f7af631a430>
To get coordinates of the river, we can use freely available OpenStreetMap data.
Thanks to OSMnx, we can automatically obtain a GeoPandas DataFrame from OSM data.
We thus only need to find out the river's Way ID. Using the OSM editor, we find out it's 34406947.
Note: To get the id when the river is selected in the OSM editor, press ctrl
+i
.
Video("data/howto-osm-id.mp4", embed=True)
osm_id = 'W34406947'
river = ox.geocode_to_gdf(osm_id, by_osmid=True)
river = river.to_crs(dem.rio.crs)
Let's plot the river
river.plot()
<AxesSubplot:>
xmin, ymin, xmax, ymax = 5.16e5, 7.6095e6, 5.23e5, 7.6155e6
bounds = box(xmin, ymin, xmax, ymax)
river = river.clip(bounds)
river_geom = river.geometry.iloc[0]
dem = dem.sel(y=slice(ymax, ymin), x=slice(xmin, xmax))
Let's check if everything looks good
fig, ax = plt.subplots()
dem.squeeze().plot.imshow(ax=ax)
river.plot(ax=ax, color='red')
<AxesSubplot:title={'center':'band = 1, spatial_ref = 0'}, xlabel='x', ylabel='y'>
To calculate the REM, we need to
Extract coordinates as DataArray
xs, ys = river_geom.xy
xs, ys = xr.DataArray(xs, dims='z'), xr.DataArray(ys, dims='z')
Use xarray's interp
to extract the river's elevation.
sampled = dem.interp(x=xs, y=ys, method='nearest').dropna(dim='z')
Prepare data for the interpolation
# Sampled river coordinates
c_sampled = np.vstack([sampled.coords[c].values for c in ('x', 'y')]).T
# All (x, y) coordinates of the original DEM
c_x, c_y = [dem.coords[c].values for c in ('x', 'y')]
c_interpolate = np.dstack(np.meshgrid(c_x, c_y)).reshape(-1, 2)
# Sampled values
values = sampled.values.ravel()
Perform the interpolation. Here we use a simple implementation of IWD which averages the 5 nearest points, weighted by inverse distance.
tree = KDTree(c_sampled)
# IWD interpolation
distances, indices = tree.query(c_interpolate, k=5)
weights = 1 / distances
weights = weights / weights.sum(axis=1).reshape(-1, 1)
interpolated_values = (weights * values[indices]).sum(axis=1)
We create a DataArray
out of the inerpolated values
elevation_raster = xr.DataArray(
interpolated_values.reshape((len(c_y), len(c_x))).T, dims=('x', 'y'), coords={'x': c_x, 'y': c_y}
)
fig, ax = plt.subplots()
elevation_raster.transpose().plot.imshow(ax=ax)
river.plot(ax=ax, color='red')
<AxesSubplot:xlabel='x', ylabel='y'>
rem = dem - elevation_raster
Let's make some pretty REM visualizations!
colors = ['#f2f7fb', '#81a8cb', '#37123d']
shade(rem.squeeze(), cmap=colors, span=[0, 10], how='linear')
We can also visualize the DEM along with the REM
a = shade(xrspatial.hillshade(dem.squeeze(), angle_altitude=1, azimuth=310), cmap=['black', 'white'], how='linear')
b = shade(rem.squeeze(), cmap=colors, span=[0, 10], how='linear', alpha=200)
stack(a, b)