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. They also accept an xr.Dataset as the first argument with band-name keyword arguments to map variables to bands (e.g.
ndvi(ds, nir='B5', red='B4')).
Load data#
To get started, we’ll import some basic packages, along with matplotlib 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 numpy as np
import xarray as xr
from xrspatial.geotiff import open_geotiff
[2]:
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap, to_rgba
import numpy as np
import xarray as xr
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)):
if bg is not None:
h, w = rgba.shape[:2]
bgc = to_rgba(bg)
bg_rgba = np.zeros((h, w, 4), dtype=np.uint8)
bg_rgba[..., :3] = (np.array(bgc[:3]) * 255).astype(np.uint8)
bg_rgba[..., 3] = 255
rgba = _stack(bg_rgba, rgba)
plt.figure(figsize=figsize)
plt.imshow(rgba)
plt.axis('off')
plt.show()
def _show_grid(layers, titles=None, ncols=2, figsize=(12, 8)):
"""Display a list of (name, rgba-or-DataArray) in a grid."""
n = len(layers)
nrows = (n + ncols - 1) // ncols
fig, axes = plt.subplots(nrows, ncols, figsize=figsize, squeeze=False)
for idx, item in enumerate(layers):
name, arr = item if isinstance(item, tuple) else (None, item)
rgba = arr if arr.ndim == 3 and arr.shape[-1] == 4 else _shade(arr)
r, c = divmod(idx, ncols)
axes[r][c].imshow(rgba)
if name: axes[r][c].set_title(name)
axes[r][c].axis('off')
for idx in range(n, nrows * ncols):
r, c = divmod(idx, ncols)
axes[r][c].axis('off')
plt.tight_layout()
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)
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. Note: you can also load bands into an xr.Dataset and pass it directly to multispectral functions with band-name keyword arguments (e.g. ndvi(ds, nir='nir', red='red')).
[3]:
SCENE_ID = "LC80030172015001LGN00"
EXTS = {
"blue": "B2",
"green": "B3",
"red": "B4",
"nir": "B5",
}
# Resample each band to a 1024x1024 grid with mean aggregation
# (replaces ds.Canvas().raster(layer, agg='mean')). We coarsen when the
# source dims divide evenly, otherwise we interpolate.
TARGET = 1024
layers = {}
for name, ext in EXTS.items():
layer = open_geotiff(f"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff", band=0)
layer.name = name
ny, nx = layer.shape
if ny >= TARGET and nx >= TARGET and ny % (ny // TARGET) == 0 and nx % (nx // TARGET) == 0:
cy, cx = ny // TARGET, nx // TARGET
layer = layer.coarsen({'y': cy, 'x': cx}, boundary='trim').mean()
else:
new_y = np.linspace(layer.y[0], layer.y[-1], TARGET)
new_x = np.linspace(layer.x[0], layer.x[-1], TARGET)
layer = layer.interp(y=new_y, x=new_x, method='linear')
layers[name] = layer
layers
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[3], line 15
11 # source dims divide evenly, otherwise we interpolate.
12 TARGET = 1024
13 layers = {}
14 for name, ext in EXTS.items():
---> 15 layer = open_geotiff(f"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff", band=0)
16 layer.name = name
17 ny, nx = layer.shape
18 if ny >= TARGET and nx >= TARGET and ny % (ny // TARGET) == 0 and nx % (nx // TARGET) == 0:
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/xrspatial/geotiff/__init__.py:1214, in open_geotiff(source, dtype, window, bbox, overview_level, band, default_name, name, chunks, gpu, max_pixels, max_cloud_bytes, on_gpu_failure, missing_sources, allow_rotated, allow_unparseable_crs, allow_invalid_nodata, stable_only, allow_experimental_codecs, allow_internal_only_jpeg, band_nodata, masked, mask_nodata, unpack, mask_and_scale, parse_coordinates, lock, cache)
1207 kwargs['max_cloud_bytes'] = max_cloud_bytes
1209 # ``read_to_array`` validates ``window`` against the selected IFD's
1210 # extent and raises ``ValueError`` for out-of-bounds windows with
1211 # the same message format as the dask path's pre-flight validator
1212 # in :func:`_read_geotiff_dask`. That keeps the two backends in sync
1213 # on the contract without forcing a second metadata parse here.
-> 1214 arr, geo_info = _read_to_array(
1215 source, window=window,
1216 overview_level=overview_level, band=band,
1217 allow_rotated=allow_rotated,
1218 allow_invalid_nodata=allow_invalid_nodata,
1219 allow_experimental_codecs=allow_experimental_codecs,
1220 allow_internal_only_jpeg=allow_internal_only_jpeg,
1221 **kwargs,
1222 )
1224 if default_name is None:
1225 # Derive from source path. File-like buffers don't have a path,
1226 # so leave name unset rather than fabricating one.
1227 if isinstance(source, str):
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/xrspatial/geotiff/_reader.py:186, in _read_to_array(source, window, overview_level, band, max_pixels, max_cloud_bytes, allow_rotated, allow_invalid_nodata, allow_experimental_codecs, allow_internal_only_jpeg)
176 raise CloudSizeLimitError(
177 f"Cloud source {source!r} is {size:,} bytes, which "
178 f"exceeds max_cloud_bytes={cloud_budget:,}. Eager "
(...) 183 f"the check, or use chunks=... for a windowed dask "
184 f"read.")
185 else:
--> 186 src = _FileSource(source)
188 sidecar = None
189 # Wrap source lifetime in the try/finally immediately after
190 # construction so ``src.close()`` runs even when ``read_all()``
191 # raises (e.g. a fsspec network failure mid-download, a transient
(...) 195 # Mirrors the close-on-error contract that ``_read_cog_http``
196 # already enforces.
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/xrspatial/geotiff/_sources.py:396, in _FileSource.__init__(self, path)
394 def __init__(self, path: str):
395 self._path = path
--> 396 self._mm, self._size, self._entry = _mmap_cache.acquire(path)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/xrspatial/geotiff/_sources.py:220, in _MmapCache.acquire(self, path)
217 self._entries.move_to_end(real)
218 return entry[1], entry[2], entry
--> 220 fh = open(real, 'rb')
221 fh.seek(0, 2)
222 size = fh.tell()
FileNotFoundError: [Errno 2] No such file or directory: '/home/docs/checkouts/readthedocs.org/user_builds/xarray-spatial/checkouts/stable/xrspatial-examples/data/LC80030172015001LGN00_B2.tiff'
Let’s do a quick visualization to see what these images look like#
[4]:
_show_grid([(name, raster) for name, raster in layers.items()], ncols=2)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[4], line 1
----> 1 _show_grid([(name, raster) for name, raster in layers.items()], ncols=2)
Cell In[2], line 55, in _show_grid(layers, titles, ncols, figsize)
51 def _show_grid(layers, titles=None, ncols=2, figsize=(12, 8)):
52 """Display a list of (name, rgba-or-DataArray) in a grid."""
53 n = len(layers)
54 nrows = (n + ncols - 1) // ncols
---> 55 fig, axes = plt.subplots(nrows, ncols, figsize=figsize, squeeze=False)
56 for idx, item in enumerate(layers):
57 name, arr = item if isinstance(item, tuple) else (None, item)
58 rgba = arr if arr.ndim == 3 and arr.shape[-1] == 4 else _shade(arr)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/pyplot.py:1887, in subplots(nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw, **fig_kw)
1884 _raise_if_figure_exists(fig_kw.get('num'), "subplots", fig_kw.get('clear'))
1886 fig = figure(**fig_kw)
-> 1887 axs = fig.subplots(nrows=nrows, ncols=ncols, sharex=sharex, sharey=sharey,
1888 squeeze=squeeze, subplot_kw=subplot_kw,
1889 gridspec_kw=gridspec_kw, height_ratios=height_ratios,
1890 width_ratios=width_ratios)
1891 return fig, axs
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/figure.py:925, in FigureBase.subplots(self, nrows, ncols, sharex, sharey, squeeze, width_ratios, height_ratios, subplot_kw, gridspec_kw)
921 raise ValueError("'width_ratios' must not be defined both as "
922 "parameter and as key in 'gridspec_kw'")
923 gridspec_kw['width_ratios'] = width_ratios
--> 925 gs = self.add_gridspec(nrows, ncols, figure=self, **gridspec_kw)
926 axs = gs.subplots(sharex=sharex, sharey=sharey, squeeze=squeeze,
927 subplot_kw=subplot_kw)
928 return axs
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/figure.py:1608, in FigureBase.add_gridspec(self, nrows, ncols, **kwargs)
1565 """
1566 Low-level API for creating a `.GridSpec` that has this figure as a parent.
1567
(...) 1604
1605 """
1607 _ = kwargs.pop('figure', None) # pop in case user has added this...
-> 1608 gs = GridSpec(nrows=nrows, ncols=ncols, figure=self, **kwargs)
1609 return gs
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/gridspec.py:364, in GridSpec.__init__(self, nrows, ncols, figure, left, bottom, right, top, wspace, hspace, width_ratios, height_ratios)
361 self.hspace = hspace
362 self.figure = figure
--> 364 super().__init__(nrows, ncols,
365 width_ratios=width_ratios,
366 height_ratios=height_ratios)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/gridspec.py:49, in GridSpecBase.__init__(self, nrows, ncols, height_ratios, width_ratios)
34 """
35 Parameters
36 ----------
(...) 46 If not given, all rows will have the same height.
47 """
48 if not isinstance(nrows, Integral) or nrows <= 0:
---> 49 raise ValueError(
50 f"Number of rows must be a positive integer, not {nrows!r}")
51 if not isinstance(ncols, Integral) or ncols <= 0:
52 raise ValueError(
53 f"Number of columns must be a positive integer, not {ncols!r}")
ValueError: Number of rows must be a positive integer, not 0
<Figure size 1200x800 with 0 Axes>
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.
[5]:
import xrspatial.multispectral as ms
tc = ms.true_color(layers["red"], layers["green"], layers["blue"])
_show(tc)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[5], line 3
1 import xrspatial.multispectral as ms
2
----> 3 tc = ms.true_color(layers["red"], layers["green"], layers["blue"])
4 _show(tc)
KeyError: 'red'
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.
[6]:
import xrspatial.multispectral as ms
from xrspatial.multispectral import ndvi, savi
nir = layers["nir"]
red = layers["red"]
ndvi_img = ndvi(nir_agg=nir, red_agg=red)
_show_grid([('nir', nir), ('red', red), ('ndvi', ndvi_img)], ncols=2)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[6], line 4
1 import xrspatial.multispectral as ms
2 from xrspatial.multispectral import ndvi, savi
3
----> 4 nir = layers["nir"]
5 red = layers["red"]
6
7 ndvi_img = ndvi(nir_agg=nir, red_agg=red)
KeyError: 'nir'
Now, substituting the blue and green bands, we get the following image.
[7]:
_show(_shade(
ndvi(nir_agg=layers["green"], red_agg=layers["blue"]),
how='eq_hist', cmap=["purple", "black", "green"],
))
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[7], line 2
1 _show(_shade(
----> 2 ndvi(nir_agg=layers["green"], red_agg=layers["blue"]),
3 how='eq_hist', cmap=["purple", "black", "green"],
4 ))
KeyError: '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.
[8]:
_show(_shade(savi(layers["nir"], layers["red"])))
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[8], line 1
----> 1 _show(_shade(savi(layers["nir"], layers["red"])))
KeyError: 'nir'
For the next few functions, we’ll experiment with an artificial terrain. We’ll use xarray-spatial’s generate_terrain to smooth the terrain and render it
Generate Terrain#
[9]:
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"
_show(_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.
[10]:
_show(_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.
[11]:
# Composite the terrain rendering and pull the RGB channels back out as
# 2D arrays to feed ndvi (replaces shade(...).to_pil() + getchannel).
terrain_rgba = _shade(terrain, cmap=Elevation, how="linear")
r, g, b, a = [
xr.DataArray(np.flipud(terrain_rgba[..., i]) / 255.0)
for i in range(4)
]
ndvi_img = ndvi(nir_agg=g, red_agg=b)
_show(_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
heightfunction to return tree heights suitable for the given elevation rangexrspatial.bumpaccepts a function with only a single argument (locations), so we will usefunctools.partialto 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
heightfunction (from xrspatial.utils import ngjit, then put@ngjitabovedef heights(...)).
[12]:
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)
_show(_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 618 ms, sys: 4.85 ms, total: 622 ms
Wall time: 622 ms
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[12], line 42
38
39 _show(_stack(
40 _shade(terrain, cmap=["black", "white"], how="linear"),
41 _shade(hillshaded, cmap=["black", "white"], how="linear", alpha=128),
---> 42 _shade(tree_colorize, cmap="limegreen", how="linear"),
43 ))
Cell In[2], line 22, in _shade(arr, cmap, alpha, min_alpha, span, how)
18 if isinstance(cmap, (list, tuple)):
19 cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
20 elif cmap is None:
21 cmap = plt.get_cmap('terrain')
---> 22 rgba = cmap(norm, bytes=True)
23 a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
24 if min_alpha is not None:
25 a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
TypeError: 'str' object is not callable