<jupyter_start><jupyter_text>Sup01-Quickly Visualize DEM Attributes with Datashader[Datashader](https://datashader.org/index.html) is a general-purpose tool for rasterizing (and re-rasterizing) data of many different types that certainly include DEM data. Moreover, Datashader provides a few geospatial-specific utilities (***datashader.geo***) that could be easily applied to calculate DEM Attributes such as:* [Slope](Slope)* [Aspect](Aspect).Datashader provides a good [user guide](https://datashader.org/user_guide/index.html) to explain key concepts of Datashader in detail. This notebook will follow the guide and apply the functions of datashader.geo to quickly calculate and visualize the DEM Attributes. To increase visuilzatioin effects, we will overlap ***Slope*** and ***Aspect*** onto DEM and its [Hillshade](Hillshade). It should be mentioned that Datashader is extremely powerful for creating meaningful representations of large datasets quickly and flexibly. In addition, it supports xarray.DataArray internally without having to convert data into numpy.array after reading DEM data.<jupyter_code>import numpy as np
import xarray as xr
from xarray import DataArray
import datashader as ds
import datashader.geo as dsgeo
from datashader.transfer_functions import shade, stack
from datashader.colors import Elevation

import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Open up the DEMSet all terrain values < 0 to nan.<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>shade(da_dem, cmap=Elevation, how='linear')<jupyter_output><empty_output><jupyter_text>Hillshade[Hillshade](https://en.wikipedia.org/wiki/Terrain_cartography) 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.<jupyter_code>illuminated = dsgeo.hillshade(da_dem)
shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')<jupyter_output><empty_output><jupyter_text>You can combine hillshading with elevation colormapping to convey differences in terrain with elevation:<jupyter_code>stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),
      shade(da_dem ,     cmap=Elevation,         alpha=128, how='linear'))<jupyter_output><empty_output><jupyter_text>Slope[Slope](https://en.wikipedia.org/wiki/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](http://wenatcheeoutdoors.org/2016/04/07/avalanche-abcs-for-snowshoers/).<jupyter_code>risky = dsgeo.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))<jupyter_output><empty_output><jupyter_text>AspectHorn (1981) calculates aspect as the direction of the maximum slope of the focal cell. The value returned is in Degrees. [Aspect](https://en.wikipedia.org/wiki/Aspect_(geography)) 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.<jupyter_code>north_faces = dsgeo.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))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10136651029981644, 'fraction_numerical': 0.014073016520497654, 'mean_word_length': 5.564926372155288, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '36636', 'license': 'MIT', 'n_tokens_mistral': 1557, 'n_tokens_neox': 1445, 'n_words': 496, 'path': 'Sup01-Quickly Visualize DEM Attributes with Datashader.ipynb', 'repo_name': 'royalosyin/Work-with-DEM-data-using-Python-from-Simple-to-Complicated', 'size': '4903'}