Sup04-Quickly Visualize DEM Attributes with xarray-spatial

xarray-spatial is a novel python package that implements common raster analysis functions using Numba and provides an easy-to-install, easy-to-extend codebase for raster analysis. It supports xarray.DataArray internally without having to convert data into numpy.array after reading DEM data.

xarray-spatial grew out of the Datashader project, which provides fast rasterization of vector data (points, lines, polygons, meshes, and rasters) for use with xarray-spatial. Therefore, we could easily apply it to calculate DEM Attributes such as:

• Slope
• Aspect.

xarray-spatial provides a good user guide to explain key functions of xarray-spatial in detail. This notebook will follow the guide and apply the spatial functions of xarray-spatial to quickly calculate the DEM Attributes and visualize them with Datashader. It should be mentioned that Datashader is extremely powerful for creating meaningful representations of large datasets quickly and flexibly.

import xarray as xr
from xrspatial import hillshade
from xrspatial import slope
from xrspatial import aspect
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.colors import Elevation

import warnings
warnings.filterwarnings("ignore")

Open up the DEM

Set all terrain values < 0 to nan.

infile  = "data/es_dem/pre_DTM.tif"
da_dem = xr.open_rasterio(infile).drop('band')[0]

# Have to check the res property.
# No support for res in both directioins of (x, y) or (lat, lon).
da_dem.attrs['res'] = da_dem.attrs['res'][0]
nodata = da_dem.nodatavals[0]
da_dem = da_dem.where(da_dem>nodata, np.nan)

shade(da_dem, cmap=['black', 'white'], how='linear')

The grayscale value above shows the elevation linearly in intensity (with the large black areas indicating low elevation), but it will look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops:

shade(da_dem, cmap=Elevation, how='linear')

Hillshade

Hillshade is a technique used to visualize terrain as shaded relief, illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which is based on slope and aspect.

illuminated = hillshade(da_dem)
shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')

You can combine hillshading with elevation colormapping to convey differences in terrain with elevation:

stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),
      shade(da_dem ,     cmap=Elevation,         alpha=128, how='linear'))

Slope

Slope is the inclination of a surface. In geography, slope is amount of change in elevation of a terrain regarding its surroundings. Horn (1981) calculates the slope of a focal cell by using a central difference estimation of a surface fitted to the focal cell and its neighbours. The slope chosen is the maximum of this surface and can be returned in several formats.

Datashader's slope function returns slope in degrees. Below we highlight areas at risk for avalanche by looking at slopes around 38 degrees.

risky = slope(da_dem)
risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)

stack(shade(da_dem,      cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(risky,       cmap='red',              how='linear', alpha=200))

Aspect

Horn (1981) calculates aspect as the direction of the maximum slope of the focal cell. The value returned is in Degrees. Aspect) is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.

Below, we look to find slopes that face close to North.

north_faces = aspect(da_dem)
north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,
                                          north_faces.data < 10), 1, np.nan)

stack(shade(da_dem,      cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(north_faces, cmap=['aqua'],           how='linear', alpha=100))

xarray-spatial supports the following spatial functions:

• Slope
• Aspect
• Curvature
• Hillshade
• Normalized Difference Vegetation Index (NDVI)
• Focal Statistics
• Zonal Statistics
• Zonal Cross Tabulate
• Viewshed
• Proximity
• Bump Mapping
• Perlin Noise
• Procedural Terrain Generation

If you are interested, it is good to have a look at the great user-guide.

References

https://datashader.org/

https://github.com/makepath/xarray-spatial

https://github.com/makepath/xarray-spatial/blob/master/examples/user-guide.ipynb

https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-slope-works.htm

Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69, 14–47. doi:10.1109/PROC.1981.11918

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