Zonal

Importing Packages

[1]:
import numpy as np
import datashader as ds
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.transfer_functions import dynspread
from datashader.transfer_functions import set_background
from datashader.colors import Elevation

import xrspatial

Generate Terrain Data

The rest of the geo-related functions focus on raster data (or rasterized data, after a previous Datashader step that returns an Xarray object). To demonstrate using these raster-based functions, let’s generate some fake terrain as an elevation raster:

[2]:
from xrspatial import generate_terrain

W = 800
H = 600

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))
terrain = generate_terrain(canvas=cvs)

shade(terrain, cmap=['black', 'white'], how='linear')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_6416/2412491959.py in <module>
      5
      6 cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))
----> 7 terrain = generate_terrain(canvas=cvs)
      8
      9 shade(terrain, cmap=['black', 'white'], how='linear')

TypeError: generate_terrain() got an unexpected keyword argument 'canvas'

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:

[3]:
shade(terrain, cmap=Elevation, how='linear')
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_6416/3138436693.py in <module>
----> 1 shade(terrain, cmap=Elevation, how='linear')

NameError: name 'terrain' is not defined

Zonal Statistics

Zonal statistics allows for calculating summary statistics for specific areas or zones within a datashader aggregate. Zones are defined by creating an integer aggregate where the cell values are zone_ids. The output of zonal statistics is a Pandas dataframe containing summary statistics for each zone based on a value raster.

Imagine the following scenario: - You are a hiker on a six-day-trip. - The path for each day is defined by a line segement. - You wish to calculate the max and min elevations for each hiking segment as a Pandas dataframe based on an elevation dataset.

[4]:
from xrspatial import hillshade
from datashader.colors import Set1
import pandas as pd

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))

terrain = generate_terrain(canvas=cvs)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')

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

zone_df = pd.DataFrame({
   'x': [-11, -5, 4, 12, 14, 18, 19],
   'y': [-5, 4, 10, 13, 13, 13, 10],
   'trail_segement_id': [11, 12, 13, 14, 15, 16, 17]
})

zones_agg = cvs.line(zone_df, 'x', 'y', ds.sum('trail_segement_id'))
zones_shaded = dynspread(shade(zones_agg, cmap=Set1), max_px=5)

stack(illuminated_shaded, terrain_shaded, zones_shaded)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_6416/685477265.py in <module>
      5 cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))
      6
----> 7 terrain = generate_terrain(canvas=cvs)
      8 terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')
      9

TypeError: generate_terrain() got an unexpected keyword argument 'canvas'
[5]:
from xrspatial import zonal_stats
zones_agg.values = np.nan_to_num(zones_agg.values, copy=False).astype(np.int)
zonal_stats(zones_agg, terrain)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_6416/1703978268.py in <module>
      1 from xrspatial import zonal_stats
----> 2 zones_agg.values = np.nan_to_num(zones_agg.values, copy=False).astype(np.int)
      3 zonal_stats(zones_agg, terrain)

NameError: name 'zones_agg' is not defined

Calculate custom stats for each zone

[6]:
custom_stats = dict(elevation_change=lambda zone: zone.max() - zone.min(),
                    elevation_min=np.min,
                    elevation_max=np.max)

zonal_stats(zones_agg, terrain, custom_stats)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_6416/2868890181.py in <module>
      3                     elevation_max=np.max)
      4
----> 5 zonal_stats(zones_agg, terrain, custom_stats)

NameError: name 'zones_agg' is not defined

Here the zones are defined by line segments, but they can be any spatial pattern, and in particular can be any region computable as a Datashader aggregate.