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: .. csv-table:: :header: "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 `_: .. code-block:: text 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: .. code-block:: python 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. .. code-block:: python 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* 1. Spectral information can be found in different locations of the GDAL Metadata system: as band metadata, as dataset metadata, in different metadata domains. 2. Since GDAL 3.10, a standardized way exists to access wavelength and fwhm, but not the bad band multiplier. 3. Depending on data providers, spectral information may still be defined in different places of the GDAL metadata model. 4. 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: .. code-block:: text ... { 0.005800, 0.005800, 0.005800, ..., , 0.009100, 0.009100, 0.009100} { 0.46, 0.465, 0.47, ..., 2.393, 2.401, 2.409} { 1, 1, 1, ..., 1, 1, 1} Micrometers 5. 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 .. 1. Look into the QgsRasterLayer custom properties and search for keys with *enmapbox/* prefix. This is mainly relevant for GUI applications, where we need to modify spectral properties using QgsRasterLayer objects, without being able to modify the underlying raster source: .. code-block:: python layer = QgsRasterLayer(path) pam: dict = layer.customProperty('enmapbox/pam') layer.setCustomProperty('enmapbox/wavelength', [400, 420, 300]) wavelength = pam['dataset'] == ds.GetMetadata_Dict() wavelength_band1 = pam['bands'][0]['']['wavelength'] == ds.GetRasterBand(1).GetMetadata_Dict() fwhm = layer.customProperty('enmapbox/fwhm') bbl = layer.customProperty('enmapbox/bbl') 1. Search in GDAL band metadata. Search first within the IMAGERY, standard and ENVI domains: .. code-block:: python 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) 2. If still undefined, search GDAL dataset metadata. Search first within the standard and ENVI domain: .. code-block:: python 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: .. code-block:: python 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')