Multispectral

Xarray-spatial’s Multispectral tools provide a range of functions pertaining to remote sensing data such as satellite imagery. A range of functions are available to calculate various vegetation and environmental parameters from the range of band data available for an area. These functions accept and output data in the form of xarray.DataArray rasters.

• Generate terrain

• Bump

• NDVI

Load data

To get started, we’ll import some basic packages, along with several handy datashader functions, mainly for rendering.

To download the examples data, run the command xrspatial examples in your terminal. All the data will be stored in your current directory inside a folder named xrspatial-examples.

import datashader as ds
from datashader.colors import Elevation
import datashader.transfer_functions as tf
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.transfer_functions import Images, Image
from datashader.utils import orient_array
import numpy as np
import xarray as xr

The following functions apply to image data with bands in different parts of the UV/Visible/IR spectrum (multispectral), so we’ll bring in some multispectral satellite image data to work with.

Below, we loaded all of the images and transformed them into the form of an xarray DataArray to use in the Xarray-spatial functions.

SCENE_ID = "LC80030172015001LGN00"
EXTS = {
    "blue": "B2",
    "green": "B3",
    "red": "B4",
    "nir": "B5",
}

cvs = ds.Canvas(plot_width=1024, plot_height=1024)
layers = {}
for name, ext in EXTS.items():
    layer = xr.open_rasterio(f"data/{SCENE_ID}_{ext}.tiff").load()[0]
    layer.name = name
    layer = cvs.raster(layer, agg="mean")
    layer.data = orient_array(layer)
    layers[name] = layer
layers

/home/giancastro/miniconda3/envs/xrspatial/lib/python3.8/site-packages/numba/np/ufunc/parallel.py:365: NumbaWarning: The TBB threading layer requires TBB version 2019.5 or later i.e., TBB_INTERFACE_VERSION >= 11005. Found TBB_INTERFACE_VERSION = 9107. The TBB threading layer is disabled.
  warnings.warn(problem)

{'blue': <xarray.DataArray (y: 1024, x: 1024)>
 array([[7595, 7136, 7328, ..., 7543, 7628, 7602],
        [7320, 7404, 7834, ..., 7529, 7536, 7548],
        [7185, 7255, 7760, ..., 7931, 7656, 7536],
        ...,
        [7312, 7293, 7285, ..., 7510, 7519, 7542],
        [7287, 7273, 7280, ..., 7523, 7535, 7541],
        [7281, 7280, 7288, ..., 7546, 7563, 7566]], dtype=uint16)
 Coordinates:
   * x        (x) float64 6.128e+05 6.129e+05 6.129e+05 ... 6.521e+05 6.521e+05
   * y        (y) float64 6.853e+06 6.853e+06 6.853e+06 ... 6.82e+06 6.82e+06
 Attributes:
     res:      30.0
     nodata:   nan,
 'green': <xarray.DataArray (y: 1024, x: 1024)>
 array([[6943, 6270, 6521, ..., 6628, 6695, 6674],
        [6474, 6460, 7015, ..., 6592, 6581, 6600],
        [6255, 6531, 7203, ..., 7119, 6726, 6566],
        ...,
        [6280, 6270, 6269, ..., 6479, 6478, 6493],
        [6265, 6254, 6261, ..., 6478, 6497, 6498],
        [6247, 6257, 6271, ..., 6496, 6524, 6535]], dtype=uint16)
 Coordinates:
   * x        (x) float64 6.128e+05 6.129e+05 6.129e+05 ... 6.521e+05 6.521e+05
   * y        (y) float64 6.853e+06 6.853e+06 6.853e+06 ... 6.82e+06 6.82e+06
 Attributes:
     res:      30.0
     nodata:   nan,
 'red': <xarray.DataArray (y: 1024, x: 1024)>
 array([[7411, 6202, 6725, ..., 6769, 6880, 6848],
        [6552, 6618, 7775, ..., 6704, 6687, 6743],
        [6165, 6773, 7959, ..., 7638, 6921, 6644],
        ...,
        [6128, 6121, 6132, ..., 6461, 6452, 6472],
        [6109, 6103, 6115, ..., 6457, 6472, 6479],
        [6098, 6112, 6122, ..., 6480, 6520, 6527]], dtype=uint16)
 Coordinates:
   * x        (x) float64 6.128e+05 6.129e+05 6.129e+05 ... 6.521e+05 6.521e+05
   * y        (y) float64 6.853e+06 6.853e+06 6.853e+06 ... 6.82e+06 6.82e+06
 Attributes:
     res:      30.0
     nodata:   nan,
 'nir': <xarray.DataArray (y: 1024, x: 1024)>
 array([[7952, 6157, 7107, ..., 7031, 7165, 7318],
        [6720, 6876, 8937, ..., 6891, 6838, 6922],
        [5965, 7029, 9140, ..., 8479, 7247, 6733],
        ...,
        [5803, 5792, 5801, ..., 6356, 6331, 6349],
        [5787, 5780, 5784, ..., 6326, 6344, 6379],
        [5769, 5781, 5793, ..., 6384, 6454, 6454]], dtype=uint16)
 Coordinates:
   * x        (x) float64 6.128e+05 6.129e+05 6.129e+05 ... 6.521e+05 6.521e+05
   * y        (y) float64 6.853e+06 6.853e+06 6.853e+06 ... 6.82e+06 6.82e+06
 Attributes:
     res:      30.0
     nodata:   nan}

Let’s do a quick visualization to see what these images look like

shaded = []
for name, raster in layers.items():
    img = shade(raster)
    img.name = name
    shaded.append(img)

imgs = Images(*shaded)
imgs.num_cols = 2
imgs

blue	green
red	nir

True Color

Now we’re ready to apply some xarray-spatial functions.

To start, we can apply true_color to the red, green, and blue bands from above to generate a real-looking image.

import xrspatial.multispectral as ms

ms.true_color(layers["red"], layers["green"], layers["blue"])

<xarray.DataArray 'true_color' (y: 1024, x: 1024, band: 4)>
array([[[126, 149, 186, 255],
        [ 82, 104, 151, 255],
        [100, 120, 166, 255],
        ...,
        [102, 128, 182, 255],
        [106, 132, 188, 255],
        [105, 131, 186, 255]],

       [[ 94, 117, 166, 255],
        [ 97, 116, 172, 255],
        [140, 154, 201, 255],
        ...,
        [100, 125, 181, 255],
        [ 99, 124, 182, 255],
        [101, 126, 183, 255]],

       [[ 81, 103, 155, 255],
        [102, 121, 161, 255],
        [147, 166, 197, 255],
        ...,
...
        ...,
        [ 91, 118, 180, 255],
        [ 91, 117, 181, 255],
        [ 91, 118, 182, 255]],

       [[ 79, 103, 163, 255],
        [ 79, 103, 162, 255],
        [ 79, 103, 163, 255],
        ...,
        [ 91, 117, 181, 255],
        [ 91, 119, 182, 255],
        [ 92, 119, 182, 255]],

       [[ 79, 102, 163, 255],
        [ 79, 103, 163, 255],
        [ 80, 104, 163, 255],
        ...,
        [ 92, 119, 183, 255],
        [ 93, 121, 184, 255],
        [ 93, 121, 184, 255]]], dtype=uint8)
Coordinates:
  * y        (y) float64 6.853e+06 6.853e+06 6.853e+06 ... 6.82e+06 6.82e+06
  * x        (x) float64 6.128e+05 6.129e+05 6.129e+05 ... 6.521e+05 6.521e+05
  * band     (band) int64 0 1 2 3
Attributes:
    res:      30.0
    nodata:   nan

NDVI

The Normalized Difference Vegetation Index (NDVI) is a metric designed to detect regions with vegetation by measuring the difference between near-infrared (NIR) light (which vegetation reflects) and red light (which vegetation absorbs).

The NDVI ranges over [-1,+1], where -1 means more “Red” radiation while +1 means more “NIR” radiation. NDVI values close to +1.0 suggest areas dense with active green foliage, while strongly negative values suggest cloud cover or snow, and values near zero suggest open water, urban areas, or bare soil.

For our synthetic example here, we don’t have access to NIR measurements, but we can approximate the results for demonstration purposes by using the green and blue channels of a colormapped image, as those represent a difference in wavelengths similar to NIR vs. Red.

Let’s start by applying xrspatial.ndvi to the satellite band images from above.

import xrspatial.multispectral as ms
from xrspatial.multispectral import ndvi
from xrspatial.multispectral import savi

nir = layers["nir"]

red = layers["red"]


nir_img = shade(nir, cmap=["purple", "black", "green"])
nir_img.name = "nir"

red_img = shade(red, cmap=["purple", "black", "green"])
red_img.name = "red"

ndvi_img = ndvi(nir_agg=nir, red_agg=red)
ndvi_img = shade(ndvi_img, cmap=["purple", "black", "green"])
ndvi_img.name = "ndvi"

Images(nir_img, red_img, ndvi_img)

nir

red

ndvi

Now, substituting the blue and green bands, we get the following image.

tf.shade(
    ndvi(nir_agg=layers["green"], red_agg=layers["blue"]),
    how="eq_hist",
    cmap=["purple", "black", "green"],
)

As you can see, we get a similar image as before, though it is not as well-defined.

SAVI

xrspatial.savi also computes the vegetation index from the red and nir bands, but it applies a correction factor for the soil brightness.

Let’s try applying that to our bands from above.

shade(savi(layers["nir"], layers["red"]))

For the next few functions, we’ll experiment with an artificial terrain. We’ll use xarray-spatial’s generate_terrain along with datashader’s Canvas to smooth thing

Generate Terrain

from xrspatial import generate_terrain

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)
terrain.attrs["unit"] = "meter"

shade(terrain, cmap=["black", "white"], how="linear")

The grayscale values in the image above show the elevation, scaled linearly in intensity (with the large black areas indicating low elevation). This is good, but it would look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops. We can use the Elevation colormap for this.

shade(terrain, cmap=Elevation, how="linear")

Now we can generate the rgba PIL image, extract the green and blue bands, and input those into ndvi.

The result is displayed below.

rgba = stack(shade(terrain, cmap=Elevation, how="linear")).to_pil()
r, g, b, a = [
    xr.DataArray(np.flipud(np.asarray(rgba.getchannel(c)))) / 255.0
    for c in ["R", "G", "B", "A"]
]

ndvi_img = ndvi(nir_agg=g, red_agg=b)
shade(ndvi_img, cmap=["purple", "black", "green"], how="linear")

Bump

Bump mapping is a cartographic technique that can be used to create the appearance of trees or other land features, which is useful when synthesizing human-interpretable images from source data like land use classifications.

xrspatial.bump will produce a bump aggregate for adding detail to the terrain.

In this example, we will pretend the bumps are trees, and shade them with green. We’ll also use the elevation data to modulate whether there are trees and if so how tall they are.

• First, we’ll define a custom height function to return tree heights suitable for the given elevation range

• xrspatial.bump accepts a function with only a single argument (locations), so we will use functools.partial to provide values for the other arguments.

• Bump mapping isn’t normally a performance bottleneck, but if you want, you can speed it up by using Numba on your custom height function (from xrspatial.utils import ngjit, then put @ngjit above def heights(...)).

from functools import partial

from xrspatial import bump, hillshade


def heights(locations, src, src_range, height=20):
    num_bumps = locations.shape[0]
    out = np.zeros(num_bumps, dtype=np.uint16)
    for r in range(0, num_bumps):
        loc = locations[r]
        x = loc[0]
        y = loc[1]
        val = src[y, x]
        if val >= src_range[0] and val < src_range[1]:
            out[r] = height
    return out


T = 300000  # Number of trees to add per call
src = terrain.data
%time trees  = bump(W, H, count=T,    height_func=partial(heights, src=src, src_range=(1000, 1300), height=5))
trees += bump(
    W,
    H,
    count=T // 2,
    height_func=partial(heights, src=src, src_range=(1300, 1700), height=20),
)
trees += bump(
    W,
    H,
    count=T // 3,
    height_func=partial(heights, src=src, src_range=(1700, 2000), height=5),
)

tree_colorize = trees.copy()
tree_colorize.data[tree_colorize.data == 0] = np.nan
hillshaded = hillshade(terrain + trees)

stack(
    shade(terrain, cmap=["black", "white"], how="linear"),
    shade(hillshaded, cmap=["black", "white"], how="linear", alpha=128),
    shade(tree_colorize, cmap="limegreen", how="linear"),
)

CPU times: user 614 ms, sys: 0 ns, total: 614 ms
Wall time: 613 ms