L1C_P.py 44.2 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
Daniel Scheffler's avatar
Daniel Scheffler committed
2
"""Level 1C Processor:   Atmospheric correction of TOA-reflectance data."""
Daniel Scheffler's avatar
Daniel Scheffler committed
3

4
import warnings
5
6
import re
import logging
7
import dill
8
import traceback
Daniel Scheffler's avatar
Daniel Scheffler committed
9
from typing import List  # noqa F401  # flake8 issue
10
11
from time import time
import os
12

13
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
14

15
from geoarray import GeoArray
16
from py_tools_ds.geo.map_info import mapinfo2geotransform
17

18
from ..options.config import GMS_config as CFG
19
from . import geoprocessing as GEOP
Daniel Scheffler's avatar
Daniel Scheffler committed
20
from .L1B_P import L1B_object
21
from ..model.metadata import get_LayerBandsAssignment
22
from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
23
from ..io.input_reader import SRF
24
# from .cloud_masking import Cloud_Mask_Creator  # circular dependencies
25

26
from sicor.sicor_ac import ac_gms
27
from sicor.sensors import RSImage
28
from sicor.Mask import S2Mask
29
from sicor.ECMWF import download_variables
30

Daniel Scheffler's avatar
Daniel Scheffler committed
31
32
__author__ = 'Daniel Scheffler'

Daniel Scheffler's avatar
Daniel Scheffler committed
33

34
class L1C_object(L1B_object):
35
    def __init__(self, L1B_obj=None):
36
        super(L1C_object, self).__init__()
37
38
39

        if L1B_obj:
            # populate attributes
Daniel Scheffler's avatar
Daniel Scheffler committed
40
            [setattr(self, key, value) for key, value in L1B_obj.__dict__.items()]
41

42
43
44
45
46
47
48
49
        # private attributes
        self._VZA_arr = None
        self._VAA_arr = None
        self._SZA_arr = None
        self._SAA_arr = None
        self._RAA_arr = None
        self._lonlat_arr = None

50
        self.proc_level = 'L1C'
51
        self.proc_status = 'initialized'
52

53
54
55
    @property
    def lonlat_arr(self):
        """Calculates pixelwise 2D-array with longitude and latitude coordinates.
56

57
58
59
60
61
62
63
64
        :return:
        """
        if self._lonlat_arr is None:
            self.logger.info('Calculating LonLat array...')
            self._lonlat_arr = \
                GEOP.get_lonlat_coord_array(self.shape_fullArr, self.arr_pos,
                                            mapinfo2geotransform(self.meta_odict['map info']),
                                            self.meta_odict['coordinate system string'],
Daniel Scheffler's avatar
Daniel Scheffler committed
65
66
                                            meshwidth=10,  # for faster processing
                                            nodata_mask=None,  # dont overwrite areas outside the image with nodata
67
68
                                            outFill=get_outFillZeroSaturated(np.float32)[0])[0]
        return self._lonlat_arr
69

70
71
72
    @lonlat_arr.setter
    def lonlat_arr(self, lonlat_arr):
        self._lonlat_arr = lonlat_arr
73

74
75
76
77
    @lonlat_arr.deleter
    def lonlat_arr(self):
        self._lonlat_arr = None

78
79
80
81
82
83
84
85
86
87
88
89
    @property
    def VZA_arr(self):
        """Get viewing zenith angle.

        :return:
        """
        if self._VZA_arr is None:
            self.logger.info('Calculating viewing zenith array...')
            if 'ViewingAngle_arrProv' in self.meta_odict and self.meta_odict['ViewingAngle_arrProv']:
                # Sentinel-2
                self._VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['ViewingAngle_arrProv'],
                                                                          self.shape_fullArr,
Daniel Scheffler's avatar
Daniel Scheffler committed
90
                                                                          meshwidth=10,  # for faster processing
91
92
93
94
                                                                          subset=None,
                                                                          bandwise=0)
            else:
                self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
95
96
97
                                                    float(self.meta_odict['ViewingAngle']),
                                                    float(self.meta_odict['FieldOfView']),
                                                    self.logger,
Daniel Scheffler's avatar
Daniel Scheffler committed
98
                                                    nodata_mask=None,  # dont overwrite areas outside image with nodata
99
                                                    outFill=get_outFillZeroSaturated(np.float32)[0],
Daniel Scheffler's avatar
Daniel Scheffler committed
100
                                                    meshwidth=10)  # for faster processing
101
102
103
104
105
        return self._VZA_arr

    @VZA_arr.setter
    def VZA_arr(self, VZA_arr):
        self._VZA_arr = VZA_arr
106

107
108
109
110
    @VZA_arr.deleter
    def VZA_arr(self):
        self._VZA_arr = None

111
112
113
    @property
    def VAA_arr(self):
        """Get viewing azimuth angle.
114

115
116
117
118
119
120
121
122
        :return:
        """
        if self._VAA_arr is None:
            self.logger.info('Calculating viewing azimuth array...')
            if 'IncidenceAngle_arrProv' in self.meta_odict and self.meta_odict['IncidenceAngle_arrProv']:
                # Sentinel-2
                self._VAA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(self.meta_odict['IncidenceAngle_arrProv'],
                                                                          self.shape_fullArr,
Daniel Scheffler's avatar
Daniel Scheffler committed
123
                                                                          meshwidth=10,  # for faster processing
124
125
126
127
128
                                                                          subset=None,
                                                                          bandwise=0)
            else:
                # only a mean VAA is available
                if self.VAA_mean is None:
129
130
                    self.VAA_mean = \
                        GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
131
132
                    assert isinstance(self.VAA_mean, float)

133
                self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
134
135
136
137
138
        return self._VAA_arr

    @VAA_arr.setter
    def VAA_arr(self, VAA_arr):
        self._VAA_arr = VAA_arr
139

140
141
142
143
    @VAA_arr.deleter
    def VAA_arr(self):
        self._VAA_arr = None

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    @property
    def SZA_arr(self):
        """Get solar zenith angle.

        :return:
        """
        if self._SZA_arr is None:
            self.logger.info('Calculating solar zenith and azimuth arrays...')
            self._SZA_arr, self._SAA_arr = \
                GEOP.calc_SZA_SAA_array(
                    self.shape_fullArr, self.arr_pos,
                    self.meta_odict['AcqDate'],
                    self.meta_odict['AcqTime'],
                    self.fullSceneCornerPos,
                    self.fullSceneCornerLonLat,
                    self.meta_odict['overpass duraction sec'],
                    self.logger,
                    meshwidth=10,
                    nodata_mask=None,  # dont overwrite areas outside the image with nodata
                    outFill=get_outFillZeroSaturated(np.float32)[0],
164
165
                    accurracy=CFG.SZA_SAA_calculation_accurracy,
                    lonlat_arr=self.lonlat_arr if CFG.SZA_SAA_calculation_accurracy == 'fine' else None)
166
167
168
169
170
171
        return self._SZA_arr

    @SZA_arr.setter
    def SZA_arr(self, SZA_arr):
        self._SZA_arr = SZA_arr

172
173
174
175
    @SZA_arr.deleter
    def SZA_arr(self):
        self._SZA_arr = None

176
177
178
179
180
181
182
    @property
    def SAA_arr(self):
        """Get solar azimuth angle.

        :return:
        """
        if self._SAA_arr is None:
183
184
            # noinspection PyStatementEffect
            self.SZA_arr  # getter also sets self._SAA_arr
185
186
187
188
189
190
        return self._SAA_arr

    @SAA_arr.setter
    def SAA_arr(self, SAA_arr):
        self._SAA_arr = SAA_arr

191
192
193
194
    @SAA_arr.deleter
    def SAA_arr(self):
        self._SAA_arr = None

195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
    @property
    def RAA_arr(self):
        """Get relative azimuth angle.

        :return:
        """
        if self._RAA_arr is None:
            self.logger.info('Calculating relative azimuth array...')
            self._RAA_arr = GEOP.calc_RAA_array(self.SAA_arr, self.VAA_mean,
                                                nodata_mask=None, outFill=get_outFillZeroSaturated(np.float32)[0])
        return self._RAA_arr

    @RAA_arr.setter
    def RAA_arr(self, RAA_arr):
        self._RAA_arr = RAA_arr
210

211
212
213
214
    @RAA_arr.deleter
    def RAA_arr(self):
        self._RAA_arr = None

215
    def delete_ac_input_arrays(self):
216
217
218
219
220
221
222
        """Delete AC input arrays if they are not needed anymore."""
        self.logger.info('Deleting input arrays for atmospheric correction...')
        del self.VZA_arr
        del self.SZA_arr
        del self.SAA_arr
        del self.RAA_arr
        del self.lonlat_arr
Daniel Scheffler's avatar
Daniel Scheffler committed
223
224
225
226
227

        # use self.dem deleter
        # would have to be resampled when writing MGRS tiles
        # -> better to directly warp it to the output dims and projection
        del self.dem
228
229
230


class AtmCorr(object):
231
    def __init__(self, *L1C_objs, reporting=False):
232
        """Wrapper around atmospheric correction by Andre Hollstein, GFZ Potsdam
233
234
235
236
237
238

        Creates the input arguments for atmospheric correction from one or multiple L1C_object instance(s) belonging to
        the same scene ID, performs the atmospheric correction and returns the atmospherically corrected L1C object(s).

        :param L1C_objs: one or more instances of L1C_object belonging to the same scene ID
        """
239
        # FIXME not yet usable for data < 2012 due to missing ECMWF archive
240
241
242
        L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)

        # hidden attributes
Daniel Scheffler's avatar
Daniel Scheffler committed
243
244
245
        self._logger = None
        self._GSDs = []
        self._data = {}
246
        self._metadata = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
247
        self._nodata = {}
248
        self._band_spatial_sampling = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
249
        self._options = {}
250
251
252

        # assertions
        scene_IDs = [obj.scene_ID for obj in L1C_objs]
Daniel Scheffler's avatar
Daniel Scheffler committed
253
        assert len(list(set(scene_IDs))) == 1, \
Daniel Scheffler's avatar
Daniel Scheffler committed
254
            "Input GMS objects for 'AtmCorr' must all belong to the same scene ID!. Received %s." % scene_IDs
255

Daniel Scheffler's avatar
Daniel Scheffler committed
256
        self.inObjs = L1C_objs  # type: List[L1C_object]
257
        self.reporting = reporting
Daniel Scheffler's avatar
Daniel Scheffler committed
258
259
        self.ac_input = {}  # set by self.run_atmospheric_correction()
        self.results = None  # direct output of external atmCorr module (set by run_atmospheric_correction)
260
        self.proc_info = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
261
        self.outObjs = []  # atmospherically corrected L1C objects
262
263

        # append AtmCorr object to input L1C objects
Daniel Scheffler's avatar
Daniel Scheffler committed
264
        # [setattr(L1C_obj, 'AtmCorr', self) for L1C_obj in self.inObjs] # too big for serialization
265

266
        if not re.search('Sentinel-2', self.inObjs[0].satellite, re.I):
Daniel Scheffler's avatar
Daniel Scheffler committed
267
268
            self.logger.warning('Calculation of acquisition geometry arrays is currently only validated for '
                                'Sentinel-2!')
269
270
            # validation possible by comparing S2 angles provided by ESA with own angles

271
272
273
274
275
    @property
    def logger(self):
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
276
            if len(self.inObjs) == 1:
277
278
279
280
281
282
283
284
285
286
                # just use the logger of the inObj
                logger_atmCorr = self.inObjs[0].logger
            else:
                # in case of multiple GMS objects to be processed at once:
                # get the logger of the first inObj
                logger_atmCorr = self.inObjs[0].logger

                # add additional file handlers for the remaining inObj (that belong to the same scene_ID)
                for inObj in self.inObjs[1:]:
                    path_logfile = inObj.pathGen.get_path_logfile()
Daniel Scheffler's avatar
Daniel Scheffler committed
287
                    fileHandler = logging.FileHandler(path_logfile, mode='a')
288
                    fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
289
                    fileHandler.setLevel(CFG.log_level)
290
291
292

                    logger_atmCorr.addHandler(fileHandler)

Daniel Scheffler's avatar
Daniel Scheffler committed
293
294
                    inObj.close_GMS_loggers()

295
296
297
298
299
300
            self._logger = logger_atmCorr
            return self._logger

    @logger.setter
    def logger(self, logger):
        assert isinstance(logger, logging.Logger) or logger in ['not set', None], \
Daniel Scheffler's avatar
Daniel Scheffler committed
301
            "AtmCorr.logger can not be set to %s." % logger
302
303
304
305
306
307
308
309
        if logger in ['not set', None]:
            self._logger.close()
            self._logger = logger
        else:
            self._logger = logger

    @logger.deleter
    def logger(self):
310
311
312
        if self._logger not in [None, 'not set']:
            self._logger.close()
            self._logger = None
313

Daniel Scheffler's avatar
Daniel Scheffler committed
314
315
        [inObj.close_GMS_loggers() for inObj in self.inObjs]

316
317
318
319
320
321
322
323
    @property
    def GSDs(self):
        """
        Returns a list of spatial samplings within the input GMS objects, e.g. [10,20,60].
        """
        for obj in self.inObjs:
            if obj.arr.xgsd != obj.arr.ygsd:
                warnings.warn("X/Y GSD is not equal for entity ID %s" % obj.entity_ID +
Daniel Scheffler's avatar
Daniel Scheffler committed
324
                              (' (%s)' % obj.subsystem if obj.subsystem else '') +
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
                              'Using X-GSD as key for spatial sampling dictionary.')
                self._GSDs.append(obj.arr.xgsd)

        return self._GSDs

    @property
    def data(self):
        """

        :return:
            ___ attribute: data, type:<class 'dict'>
            ______ key:B05, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 085998540.0803833 ]]
            ______ key:B01, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 131225590.13208008]]
            ______ key:B06, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .14965820.13977051]]
            ______ key:B11, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .11492920.10192871]]
            ______ key:B02, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 104187010.10308838]]
            ______ key:B10, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 013099670.01300049]]
            ______ key:B08, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .16857910.15783691]]
            ______ key:B04, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 065490720.06228638]]
            ______ key:B03, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 082702640.08148193]]
            ______ key:B12, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 068420410.06060791]]
            ______ key:B8A, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 192138670.17553711]]
            ______ key:B09, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] .09600830.09887695]]
            ______ key:B07, value_type:<class 'numpy.ndarray'>, repr: [[nan nan nan ...,0. [..] 173339840.15600586]]
        """
        if not self._data:
351
352
            data_dict = {}

353
            for inObj in self.inObjs:
354
                for bandN, bandIdx in inObj.arr.bandnames.items():
355
                    if bandN not in data_dict:
Daniel Scheffler's avatar
Daniel Scheffler committed
356
357
358
359
                        # float32! -> conversion to np.float16 will convert -9999 to -10000
                        arr2pass = inObj.arr[:, :, bandIdx].astype(np.float32)
                        arr2pass[arr2pass == inObj.arr.nodata] = np.nan  # set nodata values to np.nan
                        data_dict[bandN] = (arr2pass / inObj.meta_odict['ScaleFactor']).astype(np.float16)
360
                    else:
361
                        inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
Daniel Scheffler's avatar
Daniel Scheffler committed
362
                                             "exists multiple times." % bandN)
363

364
            # validate: data must have all bands needed for AC
Daniel Scheffler's avatar
Daniel Scheffler committed
365
366
            full_LBA = get_LayerBandsAssignment(self.inObjs[0].GMS_identifier, return_fullLBA=True)
            all_bNs_AC = ['B%s' % i if len(i) == 2 else 'B0%s' % i for i in full_LBA]
367
368
            if not all([bN in list(data_dict.keys()) for bN in all_bNs_AC]):
                raise RuntimeError('Atmospheric correction did not receive all the needed bands. \n\tExpected: %s;\n\t'
Daniel Scheffler's avatar
Daniel Scheffler committed
369
                                   'Received: %s' % (str(all_bNs_AC), str(list(sorted(data_dict.keys())))))
370
371
372

            self._data = data_dict

373
374
375
376
377
378
379
380
        return self._data

    @data.setter
    def data(self, data_dict):
        assert isinstance(data_dict, dict), \
            "'data' can only be set to a dictionary with band names as keys and numpy arrays as values."
        self._data = data_dict

381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
    @property
    def nodata(self):
        """

        :return:
            ___ attribute: nodata, type:<class 'dict'>
            ______ key:60.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..]  False False False]]
            ______ key:10.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..]  False False False]]
            ______ key:20.0, value_type:<class 'numpy.ndarray'>, repr: [[ TrueTrueTrue ..., [..]  False False False]]
        """

        if not self._nodata:
            for inObj in self.inObjs:
                self._nodata[inObj.arr.xgsd] = ~inObj.arr.mask_nodata[:]

        return self._nodata

    @property
    def tile_name(self):
400
        """Returns S2A tile name.
401
        NOTE: this is only needed if no DEM is passed to ac_gms
402
403
404
405
406

        :return: e.g.
            '32UMA'
        """

Daniel Scheffler's avatar
Daniel Scheffler committed
407
        return ''  # FIXME
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435

    @property
    def band_spatial_sampling(self):
        """

        :return: e.g.
            {'B01': 60.0,
             'B02': 10.0,
             'B03': 10.0,
             'B04': 10.0,
             'B05': 20.0,
             'B06': 20.0,
             'B07': 20.0,
             'B08': 10.0,
             'B09': 60.0,
             'B10': 60.0,
             'B11': 20.0,
             'B12': 20.0,
             'B8A': 20.0}
        """

        if not self._band_spatial_sampling:
            for inObj in self.inObjs:
                for bandN in inObj.arr.bandnames:
                    if bandN not in self._band_spatial_sampling:
                        self._band_spatial_sampling[bandN] = inObj.arr.xgsd
        return self._band_spatial_sampling

436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
    @property
    def metadata(self):
        """

        :return:
            ___ attribute: metadata, type:<class 'dict'>
            ______ key:spatial_samplings
            _________ key:60.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 1830
            ____________ key:XDIM, value_type:<class 'int'>, repr: 60
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 1830
            ____________ key:YDIM, value_type:<class 'int'>, repr: -60
            _________ key:10.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 10980
            ____________ key:XDIM, value_type:<class 'int'>, repr: 10
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 10980
            ____________ key:YDIM, value_type:<class 'int'>, repr: -10
            _________ key:20.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 5490
            ____________ key:XDIM, value_type:<class 'int'>, repr: 20
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 5490
            ____________ key:YDIM, value_type:<class 'int'>, repr: -20
            ______ key:SENSING_TIME, value_type:<class 'datetime.datetime'>, repr: 2016-03-26 10:34:06.538000+00:00
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
466

467
        if not self._metadata:
Daniel Scheffler's avatar
Daniel Scheffler committed
468
            del self.logger  # otherwise each input object would have multiple fileHandlers
469

Daniel Scheffler's avatar
Daniel Scheffler committed
470
471
472
473
474
475
476
477
478
479
480
481
482
            metadata = dict(
                U=self.inObjs[0].meta_odict['EarthSunDist'],
                SENSING_TIME=self.inObjs[0].acq_datetime,
                # SENSING_TIME=datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z'),
                viewing_zenith=self._meta_get_viewing_zenith(),
                viewing_azimuth=self._meta_get_viewing_azimuth(),
                relative_viewing_azimuth=self._meta_get_relative_viewing_azimuth(),
                sun_mean_azimuth=self.inObjs[0].meta_odict['SunAzimuth'],
                sun_mean_zenith=90 - self.inObjs[0].meta_odict['SunElevation'],
                solar_irradiance=self._meta_get_solar_irradiance(),
                aux_data=self._meta_get_aux_data(),
                spatial_samplings=self._meta_get_spatial_samplings()
            )
483
484

            self._metadata = metadata
485
486
487

        return self._metadata

488
489
490
491
    @property
    def options(self):
        """Returns a dictionary containing AC options.
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
492
        # type: -> dict
493
494
495
496
        if self._options:
            return self._options
        else:
            self._options = self.inObjs[0].ac_options
Daniel Scheffler's avatar
Daniel Scheffler committed
497
            self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
498
            self._options["report"]["reporting"] = self.reporting
499
500
            return self._options

501
    def _meta_get_spatial_samplings(self):
502
503
504
        """

        :return:
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
         {10.0: {'NCOLS': 10980,
           'NROWS': 10980,
           'ULX': 499980.0,
           'ULY': 5800020.0,
           'XDIM': 10.0,
           'YDIM': -10.0},
          20.0: {'NCOLS': 5490,
           'NROWS': 5490,
           'ULX': 499980.0,
           'ULY': 5800020.0,
           'XDIM': 20.0,
           'YDIM': -20.0},
          60.0: {'NCOLS': 1830,
           'NROWS': 1830,
           'ULX': 499980.0,
           'ULY': 5800020.0,
           'XDIM': 60.0,
           'YDIM': -60.0}}
523
        """
524
525
        # set corner coordinates and dims
        spatial_samplings = {}
526
527
528

        for inObj in self.inObjs:

529
530
531
532
533
            # validate GSD
            if inObj.arr.xgsd != inObj.arr.ygsd:
                warnings.warn("X/Y GSD is not equal for entity ID %s" % inObj.entity_ID +
                              (' (%s)' % inObj.subsystem if inObj.subsystem else '') +
                              'Using X-GSD as key for spatial sampling dictionary.')
534

535
536
            # set spatial information
            spatial_samplings[inObj.arr.xgsd] = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
537
538
539
540
541
542
                ULX=inObj.arr.box.boxMapYX[0][1],
                ULY=inObj.arr.box.boxMapYX[0][0],
                XDIM=inObj.arr.xgsd,
                YDIM=-inObj.arr.ygsd,
                NROWS=inObj.arr.rows,
                NCOLS=inObj.arr.cols)
543

544
545
546
        return spatial_samplings

    def _meta_get_solar_irradiance(self):
547
548
549
        """

        :return:
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
        {'B01': 1913.57,
         'B02': 1941.63,
         'B03': 1822.61,
         'B04': 1512.79,
         'B05': 1425.56,
         'B06': 1288.32,
         'B07': 1163.19,
         'B08': 1036.39,
         'B09': 813.04,
         'B10': 367.15,
         'B11': 245.59,
         'B12': 85.25,
         'B8A': 955.19}
        """

        solar_irradiance = {}

        for inObj in self.inObjs:
            for bandN, bandIdx in inObj.arr.bandnames.items():
                if bandN not in solar_irradiance:
                    solar_irradiance[bandN] = inObj.meta_odict['SolIrradiance'][bandIdx]
        return solar_irradiance

    def _meta_get_viewing_zenith(self):
        """

        :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
        """

        viewing_zenith = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
581
        for inObj in self.inObjs:  # type: L1C_object
582
            for bandN, bandIdx in inObj.arr.bandnames.items():
583
                if bandN not in viewing_zenith:
584
585
                    arr2pass = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim == 3 else inObj.VZA_arr
                    viewing_zenith[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
586
                    # viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
587
588
589
590
591
592
593
594
595
596
        return viewing_zenith

    def _meta_get_viewing_azimuth(self):
        """

        :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
        """

        viewing_azimuth = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
597
        for inObj in self.inObjs:  # type: L1C_object
598
            for bandN, bandIdx in inObj.arr.bandnames.items():
599
                if bandN not in viewing_azimuth:
Daniel Scheffler's avatar
Daniel Scheffler committed
600
                    arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim == 3 else inObj.VAA_arr
601
                    viewing_azimuth[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
602
                    # viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
603

604
605
606
607
608
609
        return viewing_azimuth

    def _meta_get_relative_viewing_azimuth(self):
        """

        :return: {B10:ndarray(dtype=float16),[...],B09:ndarray(dtype=float16)}
610
611
        """

612
613
        relative_viewing_azimuth = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
614
        for inObj in self.inObjs:  # type: L1C_object
615
            for bandN, bandIdx in inObj.arr.bandnames.items():
616
                if bandN not in relative_viewing_azimuth:
617
618
                    arr2pass = inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim == 3 else inObj.RAA_arr
                    relative_viewing_azimuth[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
619
620
                    # relative_viewing_azimuth[bandN] = \
                    #     inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr
621

622
        return relative_viewing_azimuth
623

624
625
626
627
628
629
630
631
    def _meta_get_aux_data(self):
        """

        :return:  {lons:ndarray(dtype=float16),,lats:ndarray(dtype=float16)}
        """

        aux_data = dict(
            # set lons and lats (a 2D array for all bands is enough (different band resolutions dont matter))
Daniel Scheffler's avatar
Daniel Scheffler committed
632
633
            lons=self.inObjs[0].lonlat_arr[::10, ::10, 0].astype(np.float16),  # 2D array of lon values: 0° - 360°
            lats=self.inObjs[0].lonlat_arr[::10, ::10, 1].astype(np.float16)  # 2D array of lat values: -90° - 90°
634
            # FIXME correct to reduce resolution here by factor 10?
635
636
637
638
639
640
641
642
643
644
645
646
        )

        return aux_data

    def _get_dem(self):
        """Get a DEM to be used in atmospheric correction.

        :return: <np.ndarray> 2D array (with 20m resolution in case of Sentinel-2)
        """
        # determine which input GMS object is used to generate DEM
        if re.search('Sentinel-2', self.inObjs[0].satellite):
            # in case of Sentinel-2 the 20m DEM must be passed
Daniel Scheffler's avatar
Daniel Scheffler committed
647
            inObj4dem = [obj for obj in self.inObjs if obj.arr.xgsd == 20]
648
649
650
            if not inObj4dem:
                self.logger.warning('Sentinel-2 20m subsystem could not be found. DEM passed to '
                                    'atmospheric correction might have wrong resolution.')
651
652
653
654
            inObj4dem = inObj4dem[0]
        else:
            inObj4dem = self.inObjs[0]

655
656
657
658
        try:
            dem = inObj4dem.dem[:].astype(np.float32)
        except Exception as e:
            dem = None
Daniel Scheffler's avatar
Daniel Scheffler committed
659
            self.logger.warning('A static elevation is assumed during atmospheric correction due to an error during '
660
661
662
                                'creation of the DEM corresponding to scene %s (entity ID: %s). Error message was: '
                                '\n%s\n' % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
            self.logger.info("Print traceback in case you care:")
663
            self.logger.warning(traceback.format_exc())
664
665

        return dem
666
667

    def _get_srf(self):
668
        """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
669
        """
670
671
672
        # FIXME calculation of center wavelengths within SRF() used not the GMS algorithm
        # SRF instance must be created for all bands and the previous proc level
        GMS_identifier_fullScene = self.inObjs[0].GMS_identifier
Daniel Scheffler's avatar
Daniel Scheffler committed
673
        GMS_identifier_fullScene['Subsystem'] = ''
674
675
676
        GMS_identifier_fullScene['proc_level'] = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]

        return SRF(GMS_identifier_fullScene, wvl_unit='nanometers', format_bandnames=True)
677

678
679
680
681
682
683
    def _get_mask_clouds(self):
        """Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned.

        :return:
        """

684
685
        tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']

686
687
        # check if input GMS objects provide a cloud mask
        avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs}
688
        no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
689
690

        # compute cloud mask if not already provided
691
        if no_avail_CMs:
692
            algorithm = CFG.cloud_masking_algorithm[self.inObjs[0].satellite]
693

694
695
            if algorithm == 'SICOR':
                return None
696

697
698
699
700
701
            else:
                # FMASK or Classical Bayesian
                try:
                    from .cloud_masking import Cloud_Mask_Creator

702
                    CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm, tempdir_root=CFG.path_tempdir)
703
704
705
706
                    CMC.calc_cloud_mask()
                    cm_geoarray = CMC.cloud_mask_geoarray
                    cm_array = CMC.cloud_mask_array
                    cm_legend = CMC.cloud_mask_legend
Daniel Scheffler's avatar
Daniel Scheffler committed
707
                except Exception:
708
709
                    self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
                    self.logger.error(traceback.format_exc())
710
                    return None
711

712
713
        else:
            # check if there is a cloud mask with suitable GSD
Daniel Scheffler's avatar
Daniel Scheffler committed
714
            inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd == tgt_res]
715
716
            if not inObjs2use:
                raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
Daniel Scheffler's avatar
Daniel Scheffler committed
717
                                 'GMS object provides a cloud mask with spatial resolution of %s.' % tgt_res)
718
719
720
721
722
723
724
725
            inObj2use = inObjs2use[0]

            # get mask (geo)array
            cm_geoarray = inObj2use.mask_clouds
            cm_array = inObj2use.mask_clouds[:]

            # get legend
            cm_legend = get_mask_classdefinition('mask_clouds', inObj2use.satellite)
726
            #    {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40}  # FIXME hardcoded
727
728
729
730
731
732

            # validate that xGSD equals yGSD
            if cm_geoarray.xgsd != cm_geoarray.ygsd:
                warnings.warn("Cloud mask X/Y GSD is not equal for entity ID %s" % inObj2use.entity_ID +
                              (' (%s)' % inObj2use.subsystem if inObj2use.subsystem else '') +
                              'Using X-GSD as key for cloud mask geocoding.')
733
734
735
736

        # get geocoding
        cm_geocoding = self.metadata["spatial_samplings"][tgt_res]

737
738
        # get nodata value
        self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
739

740
        # append cloud mask to input object with the same spatial resolution if there was no mask before
741
        for inObj in self.inObjs:
742
            if inObj.arr.xgsd == cm_geoarray.xgsd:
743
744
                inObj.mask_clouds = cm_geoarray
                inObj.build_combined_masks_array()
Daniel Scheffler's avatar
Daniel Scheffler committed
745
                break  # appending it to one inObj is enough
746

747
748
749
        return S2Mask(mask_array=cm_array,
                      mask_legend=cm_legend,
                      geo_coding=cm_geocoding)
750

751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
    def _check_or_download_ECMWF_data(self):
        """Check if ECMWF files are already downloaded. If not, start the downloader."""
        self.logger.info('Checking if ECMWF data are available... (if not, run download!)')

        default_products = [
            "fc_T2M",
            "fc_O3",
            "fc_SLP",
            "fc_TCWV",
            "fc_GMES_ozone",
            "fc_total_AOT_550nm",
            "fc_sulphate_AOT_550nm",
            "fc_black_carbon_AOT_550nm",
            "fc_dust_AOT_550nm",
            "fc_organic_matter_AOT_550nm",
            "fc_sea_salt_AOT_550nm"]

        try:
            t0 = time()
            results = download_variables(date_from=self.inObjs[0].acq_datetime,
                                         date_to=self.inObjs[0].acq_datetime,
772
                                         db_path=CFG.path_ECMWF_db,
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
                                         max_step=120,  # default
                                         ecmwf_variables=default_products,
                                         processes=0,  # singleprocessing
                                         force=False)  # dont force download if files already exist
            t1 = time()
            self.logger.info("Runtime: %.2f" % (t1 - t0))
            for result in results:
                self.logger.info(result)

        except Exception as err:
            self.logger.error("ECMWF data download failed for scene %s (entity ID: %s). Traceback: "
                              % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID))
            self.logger.error(traceback.format_exc())

    def _validate_snr_source(self):
        """Check if the given file path for the SNR model exists - if not, use a constant SNR of 500."""
        if not os.path.isfile(self.options["uncertainties"]["snr_model"]):
790
791
            self.logger.warning('No valid SNR model found for %s %s. Using constant SNR to compute uncertainties of '
                                'atmospheric correction.' % (self.inObjs[0].satellite, self.inObjs[0].sensor))
792
793
794
            # self.options["uncertainties"]["snr_model"] = np.nan  # causes the computed uncertainties to be np.nan
            self.options["uncertainties"]["snr_model"] = 500  # use a constant SNR of 500 to compute uncertainties

795
796
    def run_atmospheric_correction(self, dump_ac_input=False):
        # type: (bool) -> list
797
798
799
        """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
        containing surface reflectance.

800
801
        :param dump_ac_input:   allows to dump the inputs of AC to the scene's processing folder in case AC fails
        :return:                list of L1C_object instances containing atmospherically corrected data
802
        """
803
804

        # collect input args/kwargs for AC
805
806
        self.logger.info('Calculating input data for atmospheric correction...')

807
        rs_data = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
808
809
810
811
812
813
814
815
816
            data=self.data,
            metadata=self.metadata,
            nodata=self.nodata,
            band_spatial_sampling=self.band_spatial_sampling,
            tile_name=self.tile_name,
            dem=self._get_dem(),
            srf=self._get_srf(),
            mask_clouds=self._get_mask_clouds()
            # returns an instance of S2Mask or None if cloud mask is not given by input GMS objects
Daniel Scheffler's avatar
Daniel Scheffler committed
817
        )  # NOTE: all keys of this dict are later converted to attributes of RSImage
818

819
820
821
822
        # remove empty values from RSImage kwargs because SICOR treats any kind of RSImage attributes as given
        # => 'None'-attributes may cause issues
        rs_data = {k: v for k, v in rs_data.items() if v is not None}

Daniel Scheffler's avatar
Daniel Scheffler committed
823
        script = False
824

825
        # check if ECMWF data are available - if not, start the download
826
        if CFG.auto_download_ecmwf:
827
            self._check_or_download_ECMWF_data()
828
829
830
831

        # validate SNR
        self._validate_snr_source()

832
833
834
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

835
        self.ac_input = dict(
836
            rs_image=rs_image,
Daniel Scheffler's avatar
Daniel Scheffler committed
837
            options=self.options,  # type: dict
838
839
            logger=repr(self.logger),  # only a string
            script=script
840
        )
841

842
843
844
845
        # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
        # with open(path_dump, 'wb') as outF:
        #     dill.dump(self.ac_input, outF)

846
        # run AC
847
        self.logger.info('Atmospheric correction started.')
848
        try:
849
            rs_image.logger = self.logger
850
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
851

852
        except Exception as e:
853
            # serialialize AC input
854
855
856
857
858
859
            if dump_ac_input:
                path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
                with open(path_dump, 'wb') as outF:
                    dill.dump(self.ac_input, outF)

                self.logger.error('An error occurred during atmospheric correction. Inputs have been dumped to %s.'
Daniel Scheffler's avatar
Daniel Scheffler committed
860
                                  % path_dump)
861
862

            # delete AC input arrays
Daniel Scheffler's avatar
Daniel Scheffler committed
863
            for inObj in self.inObjs:  # type: L1C_object
864
865
                inObj.delete_ac_input_arrays()

866
867
            self.logger.error('\nAn error occurred during atmospheric correction. BE AWARE THAT THE SCENE %s '
                              '(ENTITY ID %s) HAS NOT BEEN ATMOSPHERICALLY CORRECTED! Error message was: \n%s\n'
868
                              % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
869
            self.logger.error(traceback.format_exc())
870
            # TODO include that in the job summary
871

872
873
            return list(self.inObjs)

874
        # get processing infos
Daniel Scheffler's avatar
Daniel Scheffler committed
875
        self.proc_info = self.ac_input['options']['processing']  # FIXME this is not appended to GMS objects
876

877
        # join results
Daniel Scheffler's avatar
Daniel Scheffler committed
878
        self._join_results_to_inObjs()  # sets self.outObjs
879

880
881
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
882

883
884
885
        return self.outObjs

    def _join_results_to_inObjs(self):
886
887
888
        """
        Join results of atmospheric correction to the input GMS objects.
        """
889

890
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
Daniel Scheffler's avatar
Daniel Scheffler committed
891
892
893
        # delete logger
        # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
        del self.logger
894
895
896
897
898
899
900
901

        self._join_data_ac()
        self._join_data_errors()
        self._join_mask_clouds()
        self._join_mask_confidence_array()

        # update masks (always do that because masks can also only contain one layer)
        [inObj.build_combined_masks_array() for inObj in self.inObjs]
902

903
904
905
906
        self.outObjs = self.inObjs

    def _join_data_ac(self):
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
907
908
        Join ATMOSPHERICALLY CORRECTED ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from
        config.
909
        """
910

911
        if self.results.data_ac is not None:
912
            for inObj in self.inObjs:
913
                assert isinstance(inObj, L1B_object)
914
                nodata = self.results.nodata[inObj.arr.xgsd]  # 2D mask with True outside of image coverage
Daniel Scheffler's avatar
Daniel Scheffler committed
915
                ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
916
                out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
917

918
919
920
                # update metadata
                inObj.arr_desc = 'BOA_Ref'
                inObj.MetaObj.bands = len(self.results.data_ac)
921
                inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
922
923
                inObj.MetaObj.LayerBandsAssignment = out_LBA
                inObj.MetaObj.filter_layerdependent_metadata()
924
                inObj.meta_odict = inObj.MetaObj.to_odict()  # actually auto-updated by getter
925

926
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config
927
928
                # FIXME AC output nodata values = 0 -> new nodata areas but mask not updated
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
929
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
930
                surf_refl *= CFG.scale_factor_BOARef  # scale using scale factor (output is float16)
931
932
                # FIXME really set AC nodata values to GMS outZero?
                surf_refl[nodata] = oZ_refl  # overwrite AC nodata values with GMS outZero
Daniel Scheffler's avatar
Daniel Scheffler committed
933
                # apply the original nodata mask (indicating background values)
934
                surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
935

Daniel Scheffler's avatar
Daniel Scheffler committed
936
                if self.results.bad_data_value is np.nan:
937
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
938
                else:
Daniel Scheffler's avatar
Daniel Scheffler committed
939
940
                    surf_refl[
                        surf_refl == self.results.bad_data_value] = oF_refl  # FIXME meaningful to set AC nans to -9999?
941
942
943

                # overwrite LayerBandsAssignment and use inObj.arr setter to generate a GeoArray
                inObj.LayerBandsAssignment = out_LBA
944
                inObj.arr = surf_refl.astype(inObj.arr.dtype)  # -> int16 (also converts NaNs to 0 if needed
945

946
947
948
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
949

950
951
952
953
954
955
956
957
958
959
960
    def _join_data_errors(self):
        """
        Join ERRORS ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config.
        """

        if self.results.data_errors is not None:
            for inObj in self.inObjs:
                nodata = self.results.nodata[inObj.arr.xgsd]  # 2D mask with True outside of image coverage
                ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]

                ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
961
962
                ac_errors *= CFG.ac_scale_factor_errors  # scale using scale factor (output is float16)
                out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
                ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
                ac_errors = ac_errors.astype(out_dtype)
                inObj.ac_errors = ac_errors  # setter generates a GeoArray with the same bandnames like inObj.arr
                # TODO how to handle nans?
        else:
            self.logger.warning("Atmospheric correction did not provide a 'data_errors' array. Maybe due to "
                                "missing SNR model? GMS_object.ac_errors kept None.")

    def _join_mask_clouds(self):
        """
        Join CLOUD MASK as 2D uint8 array.
        NOTE: mask_clouds has also methods 'export_mask_rgb()', 'export_confidence_to_jpeg2000()', ...
        """

        if self.results.mask_clouds.mask_array is not None:
            mask_clouds_ac = self.results.mask_clouds.mask_array  # uint8 2D array
979

980
981
            joined = False
            for inObj in self.inObjs:
982
983
                # delete all previous cloud masks
                del inObj.mask_clouds
984
985
986
987

                # append mask_clouds only to the input GMS object with the same dimensions
                if inObj.arr.shape[:2] == mask_clouds_ac.shape:
                    inObj.mask_clouds = mask_clouds_ac
988
989
                    inObj.mask_clouds.legend = self.results.mask_clouds.mask_legend  # dict(value=string, string=value))
                    # FIXME legend is not used later
990
991

                    # set cloud mask nodata value
992
                    tgt_nodata = get_outFillZeroSaturated(mask_clouds_ac.dtype)[0]
993
994
                    ac_out_nodata = self.ac_input['options']['cld_mask']['nodata_value_mask']
                    if tgt_nodata not in self.results.mask_clouds.mask_legend.keys():
995
                        inObj.mask_clouds[inObj.mask_clouds[:] == ac_out_nodata] = tgt_nodata
996
997
998
999
                        mask_clouds_nodata = tgt_nodata
                    else:
                        warnings.warn('The cloud mask from AC output already uses the desired nodata value %s for the '
                                      'class %s. Using AC output nodata value %s.'
1000
                                      % (tgt_nodata, self.results.mask_clouds.mask_legend[tgt_nodata], ac_out_nodata))
1001
1002
1003
1004
                        mask_clouds_nodata = ac_out_nodata

                    inObj.mask_clouds.nodata = mask_clouds_nodata

1005
                    joined = True
1006

1007
1008
1009
1010
            if not joined:
                self.logger.warning('Cloud mask has not been appended to one of the AC inputs because there was no'
                                    'input GMS object with the same dimensions.')

1011
        else:
1012
1013
            self.logger.warning("Atmospheric correction did not provide a 'mask_clouds.mask_array' array. "
                                "GMS_object.mask_clouds kept None.")
1014

1015
1016
1017
1018
1019
1020
1021
1022
    def _join_mask_confidence_array(self):
        """
        Join confidence array for mask_clouds.
        """

        if self.results.mask_clouds.mask_confidence_array is not None:
            cfd_arr = self.results.mask_clouds.mask_confidence_array  # float32 2D array, scaled [0-1, nodata 255]
            cfd_arr[cfd_arr == self.ac_input['options']['cld_mask']['nodata_value_mask']] = -1
1023
1024
            cfd_arr = (cfd_arr * CFG.scale_factor_BOARef).astype(np.int16)
            cfd_arr[cfd_arr == -CFG.scale_factor_BOARef] = get_outFillZeroSaturated(cfd_arr.dtype)[0]
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034

            joined = False
            for inObj in self.inObjs:

                # append mask_clouds only to the input GMS object with the same dimensions
                if inObj.arr.shape[:2] == cfd_arr.shape:
                    # set cloud mask confidence array
                    inObj.mask_clouds_confidence = GeoArray(cfd_arr, inObj.arr.gt, inObj.arr.prj,
                                                            nodata=get_outFillZeroSaturated(cfd_arr.dtype)[0])
                    joined = True
1035

1036
1037
1038
            if not joined:
                self.logger.warning('Cloud mask confidence array has not been appended to one of the AC inputs because '
                                    'there was no input GMS object with the same dimensions.')
1039

1040
1041
        else:
            self.logger.warning("Atmospheric correction did not provide a 'mask_confidence_array' array for "
1042
                                "attribute 'mask_clouds. GMS_object.mask_clouds_confidence kept None.")