# -*- coding: utf-8 -*-
# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data
#
# Copyright (C) 2020 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later version.
# Please note the following exception: `gms_preprocessing` depends on tqdm, which
# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
"""
Level 2B Processor: Spectral homogenization
"""
import numpy as np
from typing import Union # noqa F401 # flake8 issue
from geoarray import GeoArray # noqa F401 # flake8 issue
from spechomo.prediction import SpectralHomogenizer
from ..options.config import GMS_config as CFG
from ..misc.definition_dicts import datasetid_to_sat_sen, sat_sen_to_datasetid
from ..model.metadata import get_LayerBandsAssignment
from .L2A_P import L2A_object
from ..model.gms_object import GMS_identifier
__author__ = 'Daniel Scheffler'
[docs]class L2B_object(L2A_object):
def __init__(self, L2A_obj=None):
super(L2B_object, self).__init__()
if L2A_obj:
# populate attributes
[setattr(self, key, value) for key, value in L2A_obj.__dict__.items()]
self.proc_level = 'L2B'
self.proc_status = 'initialized'
[docs] def spectral_homogenization(self):
"""Apply spectral homogenization, i.e., prediction of the spectral bands of the target sensor."""
#################################################################
# collect some information specifying the needed homogenization #
#################################################################
method = CFG.spechomo_method
src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
src_cwls = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment]
# FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment
tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref)
# NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC
tgt_gmsid_kw = dict(satellite=tgt_sat,
sensor=tgt_sen,
subsystem='',
image_type='RSD',
dataset_ID=src_dsID,
logger=None)
tgt_LBA = get_LayerBandsAssignment(GMS_identifier(proc_level='L2A', **tgt_gmsid_kw))
if CFG.datasetid_spectral_ref is None:
tgt_cwl = CFG.target_CWL
tgt_fwhm = CFG.target_FWHM
else:
# exclude those bands from CFG.target_CWL and CFG.target_FWHM that have been removed after AC
full_LBA = get_LayerBandsAssignment(GMS_identifier(proc_level='L1A', **tgt_gmsid_kw),
no_thermal=True,
no_pan=False,
return_fullLBA=True,
sort_by_cwl=True,
proc_level='L1A')
tgt_cwl = [dict(zip(full_LBA, CFG.target_CWL))[bN] for bN in tgt_LBA]
tgt_fwhm = [dict(zip(full_LBA, CFG.target_FWHM))[bN] for bN in tgt_LBA]
####################################################
# special cases where homogenization is not needed #
####################################################
if self.dataset_ID == CFG.datasetid_spectral_ref:
self.logger.info("Spectral homogenization has been skipped because the dataset id equals the dataset id of "
"the spectral reference sensor.")
return
if src_cwls == CFG.target_CWL or (self.satellite == tgt_sat and self.sensor == tgt_sen):
# FIXME catch the case if LayerBandsAssignments are unequal with np.take
self.logger.info("Spectral homogenization has been skipped because the current spectral characteristics "
"are already equal to the target sensor's.")
return
#################################################
# perform spectral homogenization of image data #
#################################################
from ..processing.multiproc import is_mainprocess
SpH = SpectralHomogenizer(classifier_rootDir=CFG.path_spechomo_classif,
logger=self.logger,
CPUs=CFG.CPUs if is_mainprocess() else 1)
if method == 'LI' or CFG.datasetid_spectral_ref is None or self.arr_desc != 'BOA_Ref':
# linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
# -> no classifier for that case available -> linear interpolation
if self.arr_desc != 'BOA_Ref' and CFG.target_radunit_optical == 'BOA_Ref':
self.logger.warning("Spectral homogenization with an '%s' classifier is not possible because the input "
"image is not atmospherically corrected (BOA reflectance is needed). Falling back "
"to linear spectral interpolation." % method)
im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwl, kind='linear')
if CFG.spechomo_estimate_accuracy:
self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm "
"is set to 'LI' (Linear Interpolation).")
errs = None
else:
# a known sensor has been specified as spectral reference => apply a machine learner
im, errs = SpH.predict_by_machine_learner(self.arr,
method=method,
src_satellite=self.satellite,
src_sensor=self.sensor,
src_LBA=self.LayerBandsAssignment,
tgt_satellite=tgt_sat,
tgt_sensor=tgt_sen,
tgt_LBA=tgt_LBA,
n_clusters=CFG.spechomo_n_clusters,
classif_alg=CFG.spechomo_classif_alg,
kNN_n_neighbors=CFG.spechomo_kNN_n_neighbors,
src_nodataVal=self.arr.nodata,
out_nodataVal=self.arr.nodata,
compute_errors=CFG.spechomo_estimate_accuracy,
bandwise_errors=CFG.spechomo_bandwise_accuracy,
fallback_argskwargs=dict(
arrcube=self.arr,
source_CWLs=src_cwls,
target_CWLs=tgt_cwl,
kind='linear')
)
###################
# update metadata #
###################
self.LayerBandsAssignment = tgt_LBA
self.MetaObj.CWL = dict(zip(tgt_LBA, tgt_cwl))
self.MetaObj.FWHM = dict(zip(tgt_LBA, tgt_fwhm))
self.MetaObj.bands = len(tgt_LBA)
self.arr = im # type: GeoArray
self.spec_homo_errors = errs # type: Union[np.ndarray, None] # int16, None if ~CFG.spechomo_estimate_accuracy
#########################################################################################
# perform spectral homogenization of bandwise error information from earlier processors #
#########################################################################################
if self.ac_errors and self.ac_errors.ndim == 3:
from scipy.interpolate import interp1d
self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands "
"number..")
outarr = interp1d(np.array(src_cwls), self.ac_errors,
axis=2, kind='linear', fill_value='extrapolate')(tgt_cwl)
self.ac_errors = outarr.astype(self.ac_errors.dtype)