Classification#

The classification tools let you reclassify the values in an xarray DataArray into a new set of values based on set bins.

Importing Packages#

[1]:

import numpy as np
import xarray as xr


import xrspatial

[2]:

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import numpy as np

def _shade(arr, cmap=None, alpha=None, min_alpha=None, span=None, how='linear'):
    arr = np.asarray(arr, dtype=np.float64)
    finite = np.isfinite(arr)
    if span is not None:
        lo, hi = span
    elif finite.any():
        lo = np.nanmin(arr); hi = np.nanmax(arr)
    else:
        lo = hi = 0.0
    norm = np.zeros_like(arr, dtype=np.float32)
    if hi > lo:
        norm = np.where(finite, (arr - lo) / (hi - lo), 0.0).astype(np.float32)
    if isinstance(cmap, (list, tuple)):
        cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
    elif cmap is None:
        cmap = plt.get_cmap('terrain')
    rgba = cmap(norm, bytes=True)
    a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
    if min_alpha is not None:
        a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
    rgba[..., 3] = a
    return rgba

def _stack(*layers):
    base = layers[0].astype(np.float64)
    for top in layers[1:]:
        a = (top[..., 3:4] / 255.0)
        rgb = top[..., :3] * a + base[..., :3] * (1 - a)
        out_a = np.maximum(top[..., 3], base[..., 3])
        base = np.concatenate([rgb, out_a[..., None]], axis=-1)
    return np.clip(base, 0, 255).astype(np.uint8)

def _show(rgba, bg=None, figsize=(8, 6)):
    plt.figure(figsize=figsize)
    plt.imshow(rgba)
    plt.axis('off')
    plt.show()

Elevation = LinearSegmentedColormap.from_list(
    'elevation',
    [(0.0, (0.0, 0.3, 0.6)), (0.2, (0.2, 0.6, 0.85)), (0.3, (0.55, 0.78, 0.92)),
     (0.32, (0.4, 0.73, 0.45)), (0.45, (0.3, 0.55, 0.25)), (0.6, (0.5, 0.5, 0.3)),
     (0.75, (0.65, 0.55, 0.4)), (0.9, (0.8, 0.8, 0.78)), (1.0, (1.0, 1.0, 1.0))],
    N=256)
Set1 = plt.get_cmap('Set1')
Set3 = plt.get_cmap('Set3')

Generate Terrain Data#

To test out the classification functions, we’ll need some rasterized data. We can generate an artificial digital elevation model (dem), or terrain, with xarray-spatial’s generate_terrain, with the help of matplotlib’s Canvas for the aggregation of values into a raster.

[3]:

from xrspatial import generate_terrain
from xrspatial import hillshade

W = 800
H = 600

template_terrain = xr.DataArray(np.zeros((H, W)))
x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)

terrain = generate_terrain(
    template_terrain, x_range=x_range, y_range=y_range, seed=1, zfactor=1000
)
_show(_stack(
    _shade(hillshade(terrain), cmap=["grey", "white"]),
    _shade(terrain, cmap=Elevation, alpha=128),
))
../_images/user_guide_classification_7_0.png

Reclassify#

Quantile Reclassify#

One method of reclassification is by quantile. In this method, the set of all values is divided into bins such that each bin contains the same number of cell values. This method works better for non-evenly distributed value sets.

Let’s try it out.

[4]:

from xrspatial import hillshade
from xrspatial import quantile

qcut_agg = quantile(terrain, k=15)

_show(_stack(
    _shade(hillshade(qcut_agg), cmap=["gray", "white"], alpha=255, how="linear"),
    _shade(qcut_agg, cmap=Elevation, alpha=128, how="linear"),
))
../_images/user_guide_classification_11_0.png

Equal Interval Reclassify#

Another method of reclassification is equal interval. This simply divides the entire values range by the given number of bins and assigns the vlaues into those bins based on where the value lies.

As you can see below, for our terrain, this flattens out a lot of our values since our data is not evenly distributed.

[5]:

from xrspatial.classify import equal_interval
from xrspatial import hillshade

equal_interval_agg = equal_interval(terrain, k=15)

_show(_stack(
    _shade(
        hillshade(equal_interval_agg), cmap=["gray", "white"], alpha=255, how="linear"
    ),
    _shade(equal_interval_agg, cmap=Elevation, alpha=128, how="linear"),
))
../_images/user_guide_classification_14_0.png

Natural Breaks (Jenks) Reclassify#

This is another non-linear classification method that is best for non-evenly distributed data that does not skew only towards the high or low range.

[6]:

from xrspatial.classify import natural_breaks
from xrspatial import hillshade

natural_breaks_agg = natural_breaks(terrain, num_sample=1000, k=15)

_show(_stack(
    _shade(
        hillshade(natural_breaks_agg), cmap=["gray", "white"], alpha=255, how="linear"
    ),
    _shade(natural_breaks_agg, cmap=Elevation, alpha=128, how="linear"),
))
../_images/user_guide_classification_17_0.png

Regions: Groupby Pixel-Value Connectivity#

Xarray-spatial’s regions function creates a raster with unique regions based on connected pixel areas with the same value. Regions assigns each such area a unique integer value.

For our reclassified rasters, this can easily be applied to separate out all of the equal-elevation ‘plateaus’ and name each one.

Connectivity can be set to either 4 or 8-pixel neighborhoods, with the default being 4.

[7]:

from xrspatial.zonal import regions
from xrspatial import hillshade

regions_agg = regions(equal_interval_agg, neighborhood=4)

_show(_stack(
    _shade(hillshade(regions_agg), cmap=["gray", "white"], alpha=255, how="linear"),
    _shade(regions_agg, cmap=Set1, alpha=128, how="eq_hist"),
))
../_images/user_guide_classification_20_0.png

References#