1# -*- coding: utf-8 -*- 

2 

3# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data 

4# 

5# Copyright (C) 2020 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

6# 

7# This software was developed within the context of the GeoMultiSens project funded 

8# by the German Federal Ministry of Education and Research 

9# (project grant code: 01 IS 14 010 A-C). 

10# 

11# 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# Foundation, either version 3 of the License, or (at your option) any later version. 

14# Please note the following exception: `gms_preprocessing` depends on tqdm, which 

15# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files 

16# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore". 

17# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE. 

18# 

19# This program is distributed in the hope that it will be useful, but WITHOUT 

20# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

21# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

22# details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along 

25# with this program. If not, see <http://www.gnu.org/licenses/>. 

26 

27""" 

28Level 2B Processor: Spectral homogenization 

29""" 

30 

31import numpy as np 

32from typing import Union # noqa F401 # flake8 issue 

33from geoarray import GeoArray # noqa F401 # flake8 issue 

34from spechomo.prediction import SpectralHomogenizer 

35 

36from ..options.config import GMS_config as CFG 

37from ..misc.definition_dicts import datasetid_to_sat_sen, sat_sen_to_datasetid 

38from ..model.metadata import get_LayerBandsAssignment 

39from .L2A_P import L2A_object 

40from ..model.gms_object import GMS_identifier 

41 

42__author__ = 'Daniel Scheffler' 

43 

44 

45class L2B_object(L2A_object): 

46 def __init__(self, L2A_obj=None): 

47 

48 super(L2B_object, self).__init__() 

49 

50 if L2A_obj: 

51 # 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): 

58 """Apply spectral homogenization, i.e., prediction of the spectral bands of the target sensor.""" 

59 ################################################################# 

60 # collect some information specifying the needed homogenization # 

61 ################################################################# 

62 

63 method = CFG.spechomo_method 

64 src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor) 

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 tgt_gmsid_kw = dict(satellite=tgt_sat, 

70 sensor=tgt_sen, 

71 subsystem='', 

72 image_type='RSD', 

73 dataset_ID=src_dsID, 

74 logger=None) 

75 tgt_LBA = get_LayerBandsAssignment(GMS_identifier(proc_level='L2A', **tgt_gmsid_kw)) 

76 

77 if CFG.datasetid_spectral_ref is None: 

78 tgt_cwl = CFG.target_CWL 

79 tgt_fwhm = CFG.target_FWHM 

80 else: 

81 # exclude those bands from CFG.target_CWL and CFG.target_FWHM that have been removed after AC 

82 full_LBA = get_LayerBandsAssignment(GMS_identifier(proc_level='L1A', **tgt_gmsid_kw), 

83 no_thermal=True, 

84 no_pan=False, 

85 return_fullLBA=True, 

86 sort_by_cwl=True, 

87 proc_level='L1A') 

88 tgt_cwl = [dict(zip(full_LBA, CFG.target_CWL))[bN] for bN in tgt_LBA] 

89 tgt_fwhm = [dict(zip(full_LBA, CFG.target_FWHM))[bN] for bN in tgt_LBA] 

90 

91 #################################################### 

92 # special cases where homogenization is not needed # 

93 #################################################### 

94 

95 if self.dataset_ID == CFG.datasetid_spectral_ref: 

96 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): 

101 # FIXME catch the case if LayerBandsAssignments are unequal with np.take 

102 self.logger.info("Spectral homogenization has been skipped because the current spectral characteristics " 

103 "are already equal to the target sensor's.") 

104 return 

105 

106 ################################################# 

107 # perform spectral homogenization of image data # 

108 ################################################# 

109 

110 from ..processing.multiproc import is_mainprocess 

111 SpH = SpectralHomogenizer(classifier_rootDir=CFG.path_spechomo_classif, 

112 logger=self.logger, 

113 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 # linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor) 

117 # -> no classifier for that case available -> linear interpolation 

118 if self.arr_desc != 'BOA_Ref' and CFG.target_radunit_optical == 'BOA_Ref': 

119 self.logger.warning("Spectral homogenization with an '%s' classifier is not possible because the input " 

120 "image is not atmospherically corrected (BOA reflectance is needed). Falling back " 

121 "to linear spectral interpolation." % method) 

122 

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

124 

125 if CFG.spechomo_estimate_accuracy: 

126 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 else: 

132 # a known sensor has been specified as spectral reference => apply a machine learner 

133 im, errs = SpH.predict_by_machine_learner(self.arr, 

134 method=method, 

135 src_satellite=self.satellite, 

136 src_sensor=self.sensor, 

137 src_LBA=self.LayerBandsAssignment, 

138 tgt_satellite=tgt_sat, 

139 tgt_sensor=tgt_sen, 

140 tgt_LBA=tgt_LBA, 

141 n_clusters=CFG.spechomo_n_clusters, 

142 classif_alg=CFG.spechomo_classif_alg, 

143 kNN_n_neighbors=CFG.spechomo_kNN_n_neighbors, 

144 src_nodataVal=self.arr.nodata, 

145 out_nodataVal=self.arr.nodata, 

146 compute_errors=CFG.spechomo_estimate_accuracy, 

147 bandwise_errors=CFG.spechomo_bandwise_accuracy, 

148 fallback_argskwargs=dict( 

149 arrcube=self.arr, 

150 source_CWLs=src_cwls, 

151 target_CWLs=tgt_cwl, 

152 kind='linear') 

153 ) 

154 

155 ################### 

156 # update metadata # 

157 ################### 

158 

159 self.LayerBandsAssignment = tgt_LBA 

160 self.MetaObj.CWL = dict(zip(tgt_LBA, tgt_cwl)) 

161 self.MetaObj.FWHM = dict(zip(tgt_LBA, tgt_fwhm)) 

162 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 

166 

167 ######################################################################################### 

168 # perform spectral homogenization of bandwise error information from earlier processors # 

169 ######################################################################################### 

170 

171 if self.ac_errors and self.ac_errors.ndim == 3: 

172 from scipy.interpolate import interp1d 

173 

174 self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands " 

175 "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)