L2B_P.py 8.56 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
32
import numpy as np
from scipy.interpolate import interp1d
33
from typing import Union  # noqa F401  # flake8 issue
34
from geoarray import GeoArray  # noqa F401  # flake8 issue
35
from spechomo.prediction import SpectralHomogenizer
36

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

43
__author__ = 'Daniel Scheffler'
44
45
46


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

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

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

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

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

64
65
        method = CFG.spechomo_method
        src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
Daniel Scheffler's avatar
Bugfix.    
Daniel Scheffler committed
66
        src_cwls = [float(self.MetaObj.CWL[bN]) for bN in self.MetaObj.LayerBandsAssignment]
67
        # FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment
68
        tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref)
69
        # NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC
70
        tgt_LBA = get_LayerBandsAssignment(
71
            GMS_identifier(satellite=tgt_sat, sensor=tgt_sen, subsystem='',
72
                           image_type='RSD', proc_level='L2A', dataset_ID=src_dsID, logger=None))
73

74
75
76
77
78
79
80
81
82
83
84
85
        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(satellite=tgt_sat, sensor=tgt_sen, subsystem='',
                               image_type='RSD', proc_level='L1A', dataset_ID=src_dsID, logger=None),
                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]

86
87
88
        ####################################################
        # special cases where homogenization is not needed #
        ####################################################
89

90
91
92
93
        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 refernce sensor.")
            return
94

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

101
102
103
104
        #################################################
        # perform spectral homogenization of image data #
        #################################################

105
106
107
108
        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)
109
110

        if method == 'LI' or CFG.datasetid_spectral_ref is None:
111
112
            # linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
            # -> no classifier for that case available -> linear interpolation
113
            im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwl, kind='linear')
114
115
116
117
118
119

            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
120
121
122

        else:
            # a known sensor has been specified as spectral reference => apply a machine learner
123
124
125
126
127
128
129
130
            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,
131
132
133
                                                      n_clusters=CFG.spechomo_n_clusters,
                                                      classif_alg=CFG.spechomo_classif_alg,
                                                      kNN_n_neighbors=CFG.spechomo_kNN_n_neighbors,
134
                                                      nodataVal=self.arr.nodata,
135
                                                      compute_errors=CFG.spechomo_estimate_accuracy,
136
                                                      bandwise_errors=CFG.spechomo_bandwise_accuracy,
137
                                                      fallback_argskwargs=dict(
138
                                                          args=dict(source_CWLs=src_cwls, target_CWLs=tgt_cwl,),
139
140
141
                                                          kwargs=dict(kind='linear')
                                                      ))

142
143
144
        ###################
        # update metadata #
        ###################
145

146
        self.LayerBandsAssignment = tgt_LBA
147
148
149
        self.MetaObj.CWL = dict(zip(tgt_LBA, tgt_cwl))
        self.MetaObj.FWHM = dict(zip(tgt_LBA, tgt_fwhm))
        self.MetaObj.bands = len(tgt_LBA)
150

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

154
155
156
157
158
159
160
        #########################################################################################
        # perform spectral homogenization of bandwise error information from earlier processors #
        #########################################################################################

        if self.ac_errors and self.ac_errors.ndim == 3:
            self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands "
                             "number..")
161
            outarr = interp1d(np.array(src_cwls), self.ac_errors,
162
                              axis=2, kind='linear', fill_value='extrapolate')(tgt_cwl)
163
            self.ac_errors = outarr.astype(self.ac_errors.dtype)