RFC 2: Spectral Properties

Authors:
Andreas Janz andreas.janz@geo.hu-berlin.de
Benjamin Jakimow, benjamin.jakimow@geo.hu-berlin.de

Started: 2021-12-23

Last modified: 2025-10-23

Status: Proposed

Summary

Based on the ENVI header format and the STAC Electro-Optical Extension Specification, it is proposed to implement functions or methods for reading and writing band-specific spectral properties from/to QGIS Raster Layers using the custom properties. Addressed spectral properties are:

Band Property	Examples	ENVI Header Field	STAC EO
(center) wavelength	0.38638 (float)	wavelength	eo:center_wavelength (in micrometers)
wavelength unit	micrometers, μm (string)	wavelength units	N/A
full with at half maximum (fwhm)	0.00545 (float)	fwhm	eo:full_width_half_max (in micrometers)
bad band multiplier	1 or 0 (int)	bbl	N/A

QGIS doesn’t provides a concept for spectral property management, and GDAL just recently implemented (version 3.10) support for wavelengths and fwhm. The here proposed approach integrates spectral property handling into QGIS custom properties. This allows to set/update spectral properties for each band of a QgsRasterLayer and to set spectral information for layers, whose data sources are immutable / read-only.

Exemplary API support is implemented in the enmapboxprocessing.rasterreader.RasterReader and qps.qgsrasterlayerproperties.QgsRasterLayerSpectralProperties classes

Motivation

The ENVI format is very common in the hyperspectral community. It describes spectral properties of raster images in the fields of ENVI header files:

ENVI
...
wavelength units = Micrometers
wavelength = {0.46, 0.465, 0.47, ..., 2.393, 2.401, 2.409}
fwhm = {0.0058, 0.0058, 0.0058, ..., 0.0091, 0.0091, 0.0091}
bbl = {1, 1, 1, ..., 1, 1, 1}

In GDAL, these values can be read from the ENVI metadata domain:

ds: gdal.Dataset = gdal.Open(path_envifile)

wl: str = ds.GetMetadataItem('wavelength', 'ENVI')
# note the required underscore in wavelength_units
wlu: str = ds.GetMetadataItem('wavelength_units', 'ENVI')
fwhm: str = ds.GetMetadataItem('fwhm', 'ENVI')

All values are returned as strings, e.g. wavelength as ‘{0.46, 0.465, 0.47, …, 2.393, 2.401, 2.409}’.

Since version 3.10, GDAL supports reading and writing of wavelength and fwhm at band level.

ds: gdal.Dataset = gdal.Open(path_envifile)
band: gdal.Band = ds.GetRasterBand(1)
wl = float(band.GetMetadataItem('CENTRAL_WAVELENGTH_UM', 'IMAGERY')))
fwhm: float(band.GetMetadataItem('FWHM_UM', 'IMAGERY')))

If implemented by the GDAL driver, the returned strings can be converted directly into float values. As now, this is supported by the GTiff, SENTINEL2, ENVI and VRT driver.

Problems

Spectral information can be found in different locations of the GDAL Metadata system: as band metadata, as dataset metadata, in different metadata domains.
Since GDAL 3.10, a standardized way exists to access wavelength and fwhm, but not the bad band multiplier.
Depending on data providers, spectral information may still be defined in different places of the GDAL metadata model.

GDAL is flexible to write spectral information into different locations. If not handled by the driver, GDAL will write it into a separated aux.xml file:

<Metadata domain="ENVI">
  ...
  <MDI key="fwhm">{ 0.005800, 0.005800, 0.005800, ..., , 0.009100, 0.009100, 0.009100}</MDI>
  <MDI key="wavelength">{ 0.46, 0.465, 0.47, ..., 2.393, 2.401, 2.409}</MDI>
  <MDI key="bbl">{ 1, 1, 1, ..., 1, 1, 1}</MDI>
  <MDI key="wavelength_units">Micrometers</MDI>
</Metadata>

Correcting or supplementing missing spectral information requires write permissions, which may not be available. Furthermore, while being opened in QGIS as QgsRasterLayer, modifications of the datasource’s GDAL metadata may get overridden by QGIS. This can lead to the problem that the EnMAP-Box cannot modify metadata from rasters when they are opened in QGIS / EnMAP Box.

Approach

Band-wise spectral properties like wavelength, wavelength units, full width at half maximum (fwhm), and bad band multiplier (bbl) values, are fetched with the following prioritization:

Todo

Add support for QgsRasterLayer::customProperties to overwrite GDAL properties

Search in GDAL band metadata. Search first within the IMAGERY, standard and ENVI domains:

ds = gdal.Dataset(path)
band = ds.GetRasterBand(42)

# test IMAGERY domain
fwhm = band.GetMetadataItem('FWHM_UM', 'IMAGERY')
wl =   band.GetMetadataItem('CENTRAL_WAVELENGTH_UM', 'IMAGERY')

# bad band multiplier
bbl = band.GetMetadataItem('bbl')

# test other domains
if wl is None:
    for domain in ['', 'ENVI']:
        wlu = band.GetMetadataItem('wavelength_units', domain)
        fwhm = band.GetMetadataItem('fwhm', domain)
        wl = band.GetMetadataItem('wavelength', domain)

If still undefined, search GDAL dataset metadata. Search first within the standard and ENVI domain:

ds = gdal.Dataset(path)

for domain in ['', 'ENVI']:
    # bad band multiplier
    bbl = band.GetMetadataItem('bbl')
    wlu = band.GetMetadataItem('wavelength_units', domain)
    wl = band.GetMetadataItem('wavelength', domain)
    fwhm = band.GetMetadataItem('fwhm', domain)

Note that the parseEnviString function helps to parse returned ENVI list string like ‘{1, 2, 3}’ into Python lists like [1, 2, 3].

Guide line 1:: If you need to set band-wise spectral properties for GDAL processing output: set it to the gdal.Band IMAGERY or default (‘’) domain. This way, i) the information is accessible with the GDAL API, and ii) consecutive band subsetting via gdal.Translate and gdal.BuildVrt can easily copy the band domains to the destination dataset.
Guide line 2:: If you need to set band-wise metadata for a QgsRasterLayer in an EnMAP-Box GUI application: set it to layer’s custom properties This is most flexible and secure. The spectral properties are i) available from the layer instance, ii) can be inspected in the layer property dialog, and iii) can be saved to and reloaded from the QGIS project and QML layer style files.
Guide line 3:: Do not update the GDAL Metadata model (*.aux.xml files), while the corresponding source is opened as a QgsRasterLayer in QGIS. QGIS may overwrite theses changes, e.g. when closing the layer or calculating new band statistics.

Implementation

Example implementations are given by the QgsRasterLayerSpectralProperties and RasterReader classes:

from enmapbox.qgispluginsupport.qps.qgsrasterlayerproperties import QgsRasterLayerSpectralProperties
from enmapboxprocessing.rasterreader import RasterReader

properties = QgsRasterLayerSpectralProperties(layer)
wavelengths: List[float] = properties.wavelength()
wavelength_subset = properties.bandValues([1,2,3], 'wavelength')

reader = RasterReader(layer)
wavelength = reader.wavelength(42)
fwhm = reader.fwhm(42)
badBandMultiplier = reader.badBandMultiplier(42)

The destination wavelength units (default is nanometers) can also be stated:

wavelengthInNanometers = reader.wavelength(42, 'nanometers')
wavelengthInMicrometers = reader.wavelength(42, 'micrometers')
wavelengthInMillimeters = reader.wavelength(42, 'millimeters')
wavelengthInMeter = reader.wavelength(42, 'meters')

To find the band closest to 850 nanometers use:

bandNo = reader.findWavelength(850)
bandNo = reader.findWavelength(0.850, 'micrometers')