Band math#
Band math combines two or more spectral bands into a single value per pixel that
highlights a physical property of the surface, such as vegetation vigor, water
content, or burn severity. GeoWombat exposes these operations through the
.gw accessor, so any xarray.DataArray opened with geowombat.open()
can compute an index directly, returning a new lazy xarray.DataArray.
Most indices are defined on surface reflectance, so the examples below set a
scale_factor to convert the stored digital numbers to reflectance. Passing
the sensor name (here 'rgbn') lets GeoWombat map band positions to named
bands (red, nir, …) automatically. Both can be supplied per call or, more
conveniently, through the configuration manager (see Configuration manager).
Vegetation indices#
Vegetation indices contrast the strong absorption of red light by chlorophyll with the strong reflectance of near-infrared light by leaf structure. Higher values generally indicate denser, healthier vegetation.
Enhanced Vegetation Index (EVI)#
EVI improves on the classic NDVI by correcting for atmospheric and soil-background effects, which makes it more reliable over dense canopies where NDVI saturates.
In [1]: import geowombat as gw
In [2]: from geowombat.data import rgbn
Calculate a vegetation index, returning an xarray.DataArray.
In [3]: with gw.open(rgbn) as ds:
...: print(ds)
...: evi = ds.gw.evi(sensor='rgbn', scale_factor=0.0001)
...: print(evi)
...:
<xarray.DataArray (band: 4, y: 403, x: 515)> Size: 830kB
dask.array<open_rasterio-de2373fc8e4f2835425a50d14996a676<this-array>, shape=(4, 403, 515), dtype=uint8, chunksize=(4, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) int64 32B 1 2 3 4
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/14)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan, nan, nan, nan)
_FillValue: nan
... ...
filename: /home/docs/checkouts/readthedocs.org/user_builds/geo...
resampling: nearest
DataType: Generic
AREA_OR_POINT: Area
_data_are_separate: 0
_data_are_stacked: 0
<xarray.DataArray (band: 1, y: 403, x: 515)> Size: 2MB
dask.array<where, shape=(1, 403, 515), dtype=float64, chunksize=(1, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) <U3 12B 'evi'
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/18)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan,)
_FillValue: nan
... ...
_data_are_separate: 0
_data_are_stacked: 0
pre-scaling: 0.0001
sensor: rgbn
vi: evi
drange: (0, 1)
Use the configuration context to set parameters, so the sensor and
scale_factor do not need to be repeated on every call.
In [4]: with gw.config.update(sensor='rgbn', scale_factor=0.0001):
...: with gw.open(rgbn) as ds:
...: evi = ds.gw.evi()
...: print(evi)
...:
<xarray.DataArray (band: 1, y: 403, x: 515)> Size: 2MB
dask.array<where, shape=(1, 403, 515), dtype=float64, chunksize=(1, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) <U3 12B 'evi'
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/18)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan,)
_FillValue: nan
... ...
AREA_OR_POINT: Area
_data_are_separate: 0
_data_are_stacked: 0
pre-scaling: 0.0001
vi: evi
drange: (0, 1)
Two-band Enhanced Vegetation Index (EVI2)#
EVI2 is a two-band approximation of EVI that uses only the red and near-infrared bands. It is useful for sensors that lack a blue band, or when the blue band is noisy.
In [5]: with gw.config.update(sensor='rgbn', scale_factor=0.0001):
...: with gw.open(rgbn) as ds:
...: evi2 = ds.gw.evi2()
...: print(evi2)
...:
<xarray.DataArray (band: 1, y: 403, x: 515)> Size: 2MB
dask.array<where, shape=(1, 403, 515), dtype=float64, chunksize=(1, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) <U4 16B 'evi2'
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/18)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan,)
_FillValue: nan
... ...
AREA_OR_POINT: Area
_data_are_separate: 0
_data_are_stacked: 0
pre-scaling: 0.0001
vi: evi2
drange: (0, 1)
Normalized difference indices#
Many indices share the normalized-difference form (b1 - b2) / (b1 + b2). Use
the generic xarray.DataArray.gw.norm_diff() function with any two-band
combination to build one on the fly. For example, a red/near-infrared pair
reproduces NDVI.
In [6]: with gw.config.update(sensor='rgbn'):
...: with gw.open(rgbn) as ds:
...: d = ds.gw.norm_diff('red', 'nir')
...: print(d)
...:
<xarray.DataArray (band: 1, y: 403, x: 515)> Size: 2MB
dask.array<where, shape=(1, 403, 515), dtype=float64, chunksize=(1, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) <U9 36B 'norm-diff'
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/18)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan,)
_FillValue: nan
... ...
AREA_OR_POINT: Area
_data_are_separate: 0
_data_are_stacked: 0
pre-scaling: 1.0
vi: norm-diff
drange: (-1, 1)
Tasseled cap transformations#
The tasseled cap transformation rotates the spectral bands into components that
correspond to physical scene characteristics, most commonly brightness,
greenness, and wetness. Because the rotation coefficients are sensor-specific,
supply the sensor name so GeoWombat applies the correct set.
In [7]: with gw.config.update(sensor='qb', scale_factor=0.0001):
...: with gw.open(rgbn) as ds:
...: tcap = ds.gw.tasseled_cap()
...: print(tcap)
...:
<xarray.DataArray (band: 3, y: 403, x: 515)> Size: 5MB
dask.array<transpose, shape=(3, 403, 515), dtype=float64, chunksize=(3, 64, 64), chunktype=numpy.ndarray>
Coordinates:
* band (band) <U10 120B 'brightness' 'greenness' 'wetness'
* y (y) float64 3kB 2.05e+06 2.05e+06 2.05e+06 ... 2.048e+06 2.048e+06
* x (x) float64 4kB 7.93e+05 7.93e+05 7.93e+05 ... 7.956e+05 7.956e+05
Attributes: (12/15)
transform: (5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
crs: 32618
res: (5.0, 5.0)
is_tiled: 1
nodatavals: (nan, nan, nan, nan)
_FillValue: nan
... ...
filename: /home/docs/checkouts/readthedocs.org/user_builds/geo...
resampling: nearest
DataType: Generic
AREA_OR_POINT: Area
_data_are_separate: 0
_data_are_stacked: 0
Additional useful indices are available, such as the normalized burn ratio (NBR), which measures fire severity, and the woody index (WI). For a full list of indices and their expected bands, see the API docs.