L2B_P.py 9.43 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2 3 4

# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data
#
5
# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6 7 8 9 10 11
#
# 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
12
# the terms of the GNU General Public License as published by the Free Software
13 14
# 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
15 16 17
# 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.
18 19 20 21 22 23 24 25 26
#
# 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/>.

Daniel Scheffler's avatar
Daniel Scheffler committed
27 28 29 30
"""
Level 2B Processor: Spectral homogenization
"""

31
import numpy as np
32
from typing import Union  # noqa F401  # flake8 issue
33
from geoarray import GeoArray  # noqa F401  # flake8 issue
34
from spechomo.prediction import SpectralHomogenizer
35

36
from ..options.config import GMS_config as CFG
37
from ..misc.definition_dicts import datasetid_to_sat_sen, sat_sen_to_datasetid
38
from ..model.metadata import get_LayerBandsAssignment
39
from .L2A_P import L2A_object
40
from ..model.gms_object import GMS_identifier
Daniel Scheffler's avatar
Daniel Scheffler committed
41

42
__author__ = 'Daniel Scheffler'
43 44 45


class L2B_object(L2A_object):
46 47
    def __init__(self, L2A_obj=None):

Daniel Scheffler's avatar
Daniel Scheffler committed
48
        super(L2B_object, self).__init__()
49 50 51

        if L2A_obj:
            # populate attributes
52
            [setattr(self, key, value) for key, value in L2A_obj.__dict__.items()]
53

54
        self.proc_level = 'L2B'
55
        self.proc_status = 'initialized'
56

57
    def spectral_homogenization(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
58
        """Apply spectral homogenization, i.e., prediction of the spectral bands of the target sensor."""
59 60 61
        #################################################################
        # collect some information specifying the needed homogenization #
        #################################################################
62

63 64
        method = CFG.spechomo_method
        src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
Daniel Scheffler's avatar
Bugfix.  
Daniel Scheffler committed
65
        src_cwls = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment]
66
        # FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment
67
        tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref)
68
        # NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC
69 70 71 72 73 74 75
        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))
76

77 78 79 80 81
        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
82 83 84 85 86 87
            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')
88 89 90
            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]

91 92 93
        ####################################################
        # special cases where homogenization is not needed #
        ####################################################
94

95 96
        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 "
97
                             "the spectral reference sensor.")
98
            return
99

100
        if src_cwls == CFG.target_CWL or (self.satellite == tgt_sat and self.sensor == tgt_sen):
Daniel Scheffler's avatar
Bugfix  
Daniel Scheffler committed
101
            # FIXME catch the case if LayerBandsAssignments are unequal with np.take
102 103
            self.logger.info("Spectral homogenization has been skipped because the current spectral characteristics "
                             "are already equal to the target sensor's.")
104 105
            return

106 107 108 109
        #################################################
        # perform spectral homogenization of image data #
        #################################################

110 111 112 113
        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)
114

115
        if method == 'LI' or CFG.datasetid_spectral_ref is None or self.arr_desc != 'BOA_Ref':
116 117
            # linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
            # -> no classifier for that case available -> linear interpolation
118 119 120 121 122
            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)

123
            im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwl, kind='linear')
124 125 126

            if CFG.spechomo_estimate_accuracy:
                self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm "
127
                                    "is set to 'LI' (Linear Interpolation).")
128 129

            errs = None
130 131 132

        else:
            # a known sensor has been specified as spectral reference => apply a machine learner
133 134 135 136 137 138 139 140
            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,
141 142 143
                                                      n_clusters=CFG.spechomo_n_clusters,
                                                      classif_alg=CFG.spechomo_classif_alg,
                                                      kNN_n_neighbors=CFG.spechomo_kNN_n_neighbors,
144 145
                                                      src_nodataVal=self.arr.nodata,
                                                      out_nodataVal=self.arr.nodata,
146
                                                      compute_errors=CFG.spechomo_estimate_accuracy,
147
                                                      bandwise_errors=CFG.spechomo_bandwise_accuracy,
148
                                                      fallback_argskwargs=dict(
149 150 151 152 153
                                                          arrcube=self.arr,
                                                          source_CWLs=src_cwls,
                                                          target_CWLs=tgt_cwl,
                                                          kind='linear')
                                                      )
154

155 156 157
        ###################
        # update metadata #
        ###################
158

159
        self.LayerBandsAssignment = tgt_LBA
160 161 162
        self.MetaObj.CWL = dict(zip(tgt_LBA, tgt_cwl))
        self.MetaObj.FWHM = dict(zip(tgt_LBA, tgt_fwhm))
        self.MetaObj.bands = len(tgt_LBA)
163

164
        self.arr = im  # type: GeoArray
165
        self.spec_homo_errors = errs  # type: Union[np.ndarray, None]  # int16, None if ~CFG.spechomo_estimate_accuracy
Daniel Scheffler's avatar
Daniel Scheffler committed
166

167 168 169 170 171
        #########################################################################################
        # perform spectral homogenization of bandwise error information from earlier processors #
        #########################################################################################

        if self.ac_errors and self.ac_errors.ndim == 3:
172 173
            from scipy.interpolate import interp1d

174 175
            self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands "
                             "number..")
176
            outarr = interp1d(np.array(src_cwls), self.ac_errors,
177
                              axis=2, kind='linear', fill_value='extrapolate')(tgt_cwl)
178
            self.ac_errors = outarr.astype(self.ac_errors.dtype)