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.

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.

[1]:
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.

[2]:
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)
[2]:
{'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#

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

imgs = Images(*shaded)
imgs.num_cols = 2
imgs
[3]:
blue