Source code for eolearn.features.radiometric_normalization

"""
Module for radiometric normalization

Credits:
Copyright (c) 2018-2019 Johannes Schmid (GeoVille)
Copyright (c) 2017-2019 Matej Aleksandrov, Matic Lubej, Devis Peresutti (Sinergise)

This source code is licensed under the MIT license found in the LICENSE
file in the root directory of this source tree.
"""

import numpy as np

from eolearn.core import EOTask, FeatureType
from eolearn.core.utilities import renamed_and_deprecated


[docs]class ReferenceScenesTask(EOTask): """ Creates a layer of reference scenes which have the highest fraction of valid pixels. The number of reference scenes is limited to a definable number. Contributor: Johannes Schmid, GeoVille Information Systems GmbH, 2018 :param feature: Name of the eopatch data layer. Needs to be of the FeatureType "DATA". :type feature: (FeatureType, str) or (FeatureType, str, str) :param valid_fraction_feature: Name of the layer containing the valid fraction obtained with the EOTask 'AddValidDataFraction'. Needs to be of the FeatureType "SCALAR". :type valid_fraction_feature: (FeatureType, str) :param max_scene_number: Maximum number of reference scenes taken for the creation of the composite. By default, the maximum number of scenes equals the number of time frames :type max_scene_number: int """ def __init__(self, feature, valid_fraction_feature, max_scene_number=None): self.feature = self._parse_features(feature, new_names=True, default_feature_type=FeatureType.DATA, rename_function='{}_REFERENCE'.format) self.valid_fraction_feature = self._parse_features(valid_fraction_feature, default_feature_type=FeatureType.SCALAR) self.number = max_scene_number
[docs] def execute(self, eopatch): feature_type, feature_name, new_feature_name = next(self.feature(eopatch)) valid_fraction_feature_type, valid_fraction_feature_name = next(self.valid_fraction_feature(eopatch)) valid_frac = list(eopatch[valid_fraction_feature_type][valid_fraction_feature_name].flatten()) data = eopatch[feature_type][feature_name] number = data.shape[0] if self.number is None else self.number eopatch[feature_type][new_feature_name] = np.array([data[x] for _, x in sorted(zip(valid_frac, range(data.shape[0])), reverse=True) if x <= number-1]) return eopatch
[docs]class BaseCompositingTask(EOTask): """ Base class to create a composite of reference scenes Contributor: Johannes Schmid, GeoVille Information Systems GmbH, 2018 :param feature: Feature holding the input time-series. Default type is FeatureType.DATA :type feature: (FeatureType, str) :param feature_composite: Type and name of output composite image. Default type is FeatureType.DATA_TIMELESS :type feature_composite: (FeatureType, str) :param percentile: Percentile along the time dimension used for compositing. Methods use different percentiles :type percentile: int or list :param max_index: Value used to flag indices with NaNs. Could be integer or NaN. Default is 255 :type max_index: int or NaN :param interpolation: Method used to compute percentile. Allowed values are {'geoville', 'linear', 'lower', 'higher', 'midpoint', 'nearest'}. 'geoville' interpolation performs a custom implementation, while the other methods use the numpy `percentile` function. Default is 'lower' :type interpolation: str :param no_data_value: Value in the composite assigned to non valid data points. Default is NaN :type no_data_value: float or NaN """ def __init__(self, feature, feature_composite, percentile=None, max_index=255, interpolation='lower', no_data_value=np.nan): self.feature = self._parse_features(feature, default_feature_type=FeatureType.DATA, rename_function='{}_COMPOSITE'.format) self.composite_type, self.composite_name = next( self._parse_features(feature_composite, default_feature_type=FeatureType.DATA_TIMELESS)()) self.percentile = percentile self.max_index = max_index self.interpolation = interpolation self._index_by_percentile = self._geoville_index_by_percentile \ if self.interpolation.lower() == 'geoville' else self._numpy_index_by_percentile self.no_data_value = no_data_value def _numpy_index_by_percentile(self, data, percentile): """ Calculate percentile of numpy stack and return the index of the chosen pixel. numpy percentile function is used with one of the following interpolations {'linear', 'lower', 'higher', 'midpoint', 'nearest'} """ data_perc_low = np.nanpercentile(data, percentile, axis=0, interpolation=self.interpolation) indices = np.empty(data_perc_low.shape, dtype=np.uint8) indices[:] = np.nan abs_diff = np.where(np.isnan(data_perc_low), np.inf, abs(data - data_perc_low)) indices = np.where(np.isnan(data_perc_low), self.max_index, np.nanargmin(abs_diff, axis=0)) return indices def _geoville_index_by_percentile(self, data, percentile): """ Calculate percentile of numpy stack and return the index of the chosen pixel. """ # no_obs = bn.allnan(arr_tmp["data"], axis=0) data_tmp = np.array(data, copy=True) valid_obs = np.sum(np.isfinite(data_tmp), axis=0) # replace NaN with maximum max_val = np.nanmax(data_tmp) + 1 data_tmp[np.isnan(data_tmp)] = max_val # sort - former NaNs will move to the end ind_tmp = np.argsort(data_tmp, kind="mergesort", axis=0) # desired position as well as floor and ceiling of it k_arr = (valid_obs - 1) * (percentile / 100.0) k_arr = np.where(k_arr < 0, 0, k_arr) f_arr = np.floor(k_arr + 0.5) f_arr = f_arr.astype(int) # get floor value of reference band and index band ind = f_arr.astype("int16") y_val, x_val = ind_tmp.shape[1], ind_tmp.shape[2] y_val, x_val = np.ogrid[0:y_val, 0:x_val] idx = np.where(valid_obs == 0, self.max_index, ind_tmp[ind, y_val, x_val]) return idx def _get_reference_band(self, data): """ Extract reference band from input 4D data according to compositing method :param data: 4D array from which to extract reference band (e.g. blue, maxNDVI, ..) :type data: numpy array :return: 3D array containing reference band according to compositing method """ raise NotImplementedError def _get_indices(self, data): """ Compute indices along temporal dimension corresponding to the sought percentile :param data: Input 3D array holding the reference band :type data: numpy array :return: 2D array holding the temporal index corresponding to percentile """ indices = self._index_by_percentile(data, self.percentile) return indices
[docs] def execute(self, eopatch): """ Compute composite array merging temporal frames according to the compositing method :param eopatch: eopatch holding time-series :return: eopatch with composite image of time-series """ feature_type, feature_name = next(self.feature(eopatch)) data = eopatch[feature_type][feature_name].copy() # compute band according to compositing method (e.g. blue, maxNDVI, maxNDWI) reference_bands = self._get_reference_band(data) # find temporal indices corresponding to pre-defined percentile indices = self._get_indices(reference_bands) # compute composite image selecting values along temporal dimension corresponding to percentile indices composite_image = np.empty((data.shape[1:]), np.float32) composite_image[:] = self.no_data_value for scene_id, scene in enumerate(data): composite_image = np.where(np.dstack([indices]) == scene_id, scene, composite_image) eopatch[self.composite_type][self.composite_name] = composite_image return eopatch
[docs]class BlueCompositingTask(BaseCompositingTask): """ Blue band compositing method - blue (25th percentile of the blue band) :param blue_idx: Index of blue band in `feature` array :type blue_idx: int """ def __init__(self, feature, feature_composite, blue_idx, interpolation='lower'): super().__init__(feature, feature_composite, percentile=25, interpolation=interpolation) self.blue_idx = blue_idx if not isinstance(blue_idx, int): raise ValueError('Incorrect value of blue band index specified') def _get_reference_band(self, data): """ Extract the blue band from time-series :param data: 4D array from which to extract the blue reference band :type data: numpy array :return: 3D array containing the blue reference band """ return data[..., self.blue_idx].astype("float32")
[docs]class HOTCompositingTask(BaseCompositingTask): """ HOT compositing method - HOT (Index using bands blue and red) The HOT index is defined as per Zhu, Z., & Woodcock, C. E. (2012). "Object-based cloud and cloud shadow detection in Landsat imagery." Remote Sensing of Environment, 118, 83-94. :param blue_idx: Index of blue band in `feature` array :type blue_idx: int :param red_idx: Index of red band in `feature` array :type red_idx: int """ def __init__(self, feature, feature_composite, blue_idx, red_idx, interpolation='lower'): super().__init__(feature, feature_composite, percentile=25, interpolation=interpolation) self.blue_idx = blue_idx self.red_idx = red_idx if not isinstance(blue_idx, int) or not isinstance(red_idx, int): raise ValueError('Incorrect values of blue and red band indices specified') def _get_reference_band(self, data): """ Extract the HOT band from time-series :param data: 4D array from which to extract the HOT reference band :type data: numpy array :return: 3D array containing the HOT reference band """ return data[..., self.blue_idx] - 0.5 * data[..., self.red_idx] - 0.08
[docs]class MaxNDVICompositingTask(BaseCompositingTask): """ maxNDVI compositing method - maxNDVI (temporal maximum of NDVI) :param red_idx: Index of red band in `feature` array :type red_idx: int :param nir_idx: Index of NIR band in `feature` array :type nir_idx: int """ def __init__(self, feature, feature_composite, red_idx, nir_idx, interpolation='lower'): super().__init__(feature, feature_composite, percentile=[0, 100], interpolation=interpolation) self.red_idx = red_idx self.nir_idx = nir_idx if not isinstance(nir_idx, int) or not isinstance(red_idx, int): raise ValueError('Incorrect values of red and NIR band indices specified') def _get_reference_band(self, data): """ Extract the NDVI band from time-series :param data: 4D array from which to compute the NDVI reference band :type data: numpy array :return: 3D array containing the NDVI reference band """ nir = data[..., self.nir_idx].astype("float32") red = data[..., self.red_idx].astype("float32") return (nir - red) / (nir + red) def _get_indices(self, data): median = np.nanmedian(data, axis=0) indices_min = self._index_by_percentile(data, self.percentile[0]) indices_max = self._index_by_percentile(data, self.percentile[1]) indices = np.where(median < -0.05, indices_min, indices_max) return indices
[docs]class MaxNDWICompositingTask(BaseCompositingTask): """ maxNDWI compositing method - maxNDWI (temporal maximum of NDWI) :param nir_idx: Index of NIR band in `feature` array :type nir_idx: int :param swir1_idx: Index of SWIR1 band in `feature` array :type swir1_idx: int """ def __init__(self, feature, feature_composite, nir_idx, swir1_idx, interpolation='lower'): super().__init__(feature, feature_composite, percentile=100, interpolation=interpolation) self.nir_idx = nir_idx self.swir1_idx = swir1_idx if not isinstance(nir_idx, int) or not isinstance(swir1_idx, int): raise ValueError('Incorrect values of NIR and SWIR1 band indices specified') def _get_reference_band(self, data): """ Extract the NDWI band from time-series :param data: 4D array from which to compute the NDWI reference band :type data: numpy array :return: 3D array containing the NDWI reference band """ nir = data[..., self.nir_idx].astype("float32") swir1 = data[..., self.swir1_idx].astype("float32") return (nir - swir1) / (nir + swir1)
[docs]class MaxRatioCompositingTask(BaseCompositingTask): """ maxRatio compositing method - maxRatio (temporal maximum of a ratio using bands blue, NIR and SWIR) :param blue_idx: Index of blue band in `feature` array :type blue_idx: int :param nir_idx: Index of NIR band in `feature` array :type nir_idx: int :param swir1_idx: Index of SWIR1 band in `feature` array :type swir1_idx: int """ def __init__(self, feature, feature_composite, blue_idx, nir_idx, swir1_idx, interpolation='lower'): super().__init__(feature, feature_composite, percentile=100, interpolation=interpolation) self.blue_idx = blue_idx self.nir_idx = nir_idx self.swir1_idx = swir1_idx if not isinstance(blue_idx, int) or not isinstance(nir_idx, int) or not isinstance(swir1_idx, int): raise ValueError('Incorrect values for either blue, NIR or SWIR1 band indices specified') def _get_reference_band(self, data): """ Extract the max-ratio band from time-series The max-ratio is defined as max(NIR,SWIR1)/BLUE :param data: 4D array from which to compute the max-ratio reference band :type data: numpy array :return: 3D array containing the max-ratio reference band """ blue = data[..., self.blue_idx].astype("float32") nir = data[..., self.nir_idx].astype("float32") swir1 = data[..., self.swir1_idx].astype("float32") return np.nanmax(np.array([nir, swir1]), axis=0) / blue
[docs]class HistogramMatchingTask(EOTask): """ Histogram match of each band of each scene within a time-series with respect to the corresponding band of a reference composite. Contributor: Johannes Schmid, GeoVille Information Systems GmbH, 2018 :param feature: Name of the eopatch data layer that will undergo a histogram match. Should be of the FeatureType "DATA". :type feature: (FeatureType, str) or (FeatureType, str, str) :param reference: Name of the eopatch data layer that represents the reference for the histogram match. Should be of the FeatureType "DATA_TIMELESS". :type reference: (FeatureType, str) """ def __init__(self, feature, reference): self.feature = self._parse_features(feature, new_names=True, default_feature_type=FeatureType.DATA, rename_function='{}_NORMALISED'.format) self.reference = self._parse_features(reference, default_feature_type=FeatureType.DATA_TIMELESS)
[docs] def execute(self, eopatch): """ Perform histogram matching of the time-series with respect to a reference scene :param eopatch: eopatch holding the time-series and reference data :type eopatch: EOPatch :return: The same eopatch instance with the normalised time-series """ feature_type, feature_name, new_feature_name = next(self.feature(eopatch)) reference_type, reference_name = next(self.reference(eopatch)) reference_scene = eopatch[reference_type][reference_name] # check if band dimension matches if reference_scene.shape[-1] != eopatch[feature_type][feature_name].shape[-1]: raise ValueError('Time-series and reference scene must have corresponding bands') eopatch[feature_type][new_feature_name] = np.zeros_like(eopatch[feature_type][feature_name]) for source_id, source in enumerate(eopatch[feature_type][feature_name]): # mask-out same invalid pixels src_masked = np.where(np.isnan(reference_scene), np.nan, source) ref_masked = np.where(np.isnan(source), np.nan, reference_scene) # compute statistics std_ref = np.nanstd(ref_masked, axis=(0, 1), dtype=np.float64) std_src = np.nanstd(src_masked, axis=(0, 1), dtype=np.float64) mean_ref = np.nanmean(ref_masked, axis=(0, 1), dtype=np.float64) mean_src = np.nanmean(src_masked, axis=(0, 1), dtype=np.float64) # normalise values eopatch[feature_type][new_feature_name][source_id] = \ source * (std_ref / std_src) + (mean_ref - (mean_src * (std_ref / std_src))) return eopatch
[docs]@renamed_and_deprecated class ReferenceScenes(ReferenceScenesTask): """ A deprecated version of ReferenceScenesTask """
[docs]@renamed_and_deprecated class BaseCompositing(BaseCompositingTask): """ A deprecated version of BaseCompositingTask """
[docs]@renamed_and_deprecated class BlueCompositing(BlueCompositingTask): """ A deprecated version of BlueCompositingTask """
[docs]@renamed_and_deprecated class MaxNDVICompositing(MaxNDVICompositingTask): """ A deprecated version of MaxNDVICompositingTask """
[docs]@renamed_and_deprecated class MaxNDWICompositing(MaxNDWICompositingTask): """ A deprecated version of MaxNDWICompositingTask """
[docs]@renamed_and_deprecated class HOTCompositing(HOTCompositingTask): """ A deprecated version of HOTCompositingTask """
[docs]@renamed_and_deprecated class MaxRatioCompositing(MaxRatioCompositingTask): """ A deprecated version of MaxRatioCompositingTask """
[docs]@renamed_and_deprecated class HistogramMatching(HistogramMatchingTask): """ A deprecated version of HistogramMatchingTask """