L1C_P.py 48 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
import timeout_decorator
13

14
import numpy as np
15
16
from ecmwfapi.api import APIKeyFetchError
from ecmwfapi.api import APIException as ECMWFAPIException
Daniel Scheffler's avatar
Daniel Scheffler committed
17

18
from geoarray import GeoArray
19
from py_tools_ds.geo.map_info import mapinfo2geotransform
20

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

29
from sicor.sicor_ac import ac_gms
30
from sicor.sensors import RSImage
31
from sicor.Mask import S2Mask
32
from sicor.ECMWF import download_variables
33

Daniel Scheffler's avatar
Daniel Scheffler committed
34
35
__author__ = 'Daniel Scheffler'

Daniel Scheffler's avatar
Daniel Scheffler committed
36

37
class L1C_object(L1B_object):
38
    def __init__(self, L1B_obj=None):
39
        super(L1C_object, self).__init__()
40
41
42

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

45
46
47
48
49
50
51
52
        # 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

53
        self.proc_level = 'L1C'
54
        self.proc_status = 'initialized'
55

56
57
58
    @property
    def lonlat_arr(self):
        """Calculates pixelwise 2D-array with longitude and latitude coordinates.
59

60
61
62
63
64
65
66
67
        :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
68
69
                                            meshwidth=10,  # for faster processing
                                            nodata_mask=None,  # dont overwrite areas outside the image with nodata
70
71
                                            outFill=get_outFillZeroSaturated(np.float32)[0])[0]
        return self._lonlat_arr
72

73
74
75
    @lonlat_arr.setter
    def lonlat_arr(self, lonlat_arr):
        self._lonlat_arr = lonlat_arr
76

77
78
79
80
    @lonlat_arr.deleter
    def lonlat_arr(self):
        self._lonlat_arr = None

81
82
83
84
85
86
87
88
89
90
91
92
    @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
93
                                                                          meshwidth=10,  # for faster processing
94
95
96
97
                                                                          subset=None,
                                                                          bandwise=0)
            else:
                self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
98
99
100
                                                    float(self.meta_odict['ViewingAngle']),
                                                    float(self.meta_odict['FieldOfView']),
                                                    self.logger,
Daniel Scheffler's avatar
Daniel Scheffler committed
101
                                                    nodata_mask=None,  # dont overwrite areas outside image with nodata
102
                                                    outFill=get_outFillZeroSaturated(np.float32)[0],
Daniel Scheffler's avatar
Daniel Scheffler committed
103
                                                    meshwidth=10)  # for faster processing
104
105
106
107
108
        return self._VZA_arr

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

110
111
112
113
    @VZA_arr.deleter
    def VZA_arr(self):
        self._VZA_arr = None

114
115
116
    @property
    def VAA_arr(self):
        """Get viewing azimuth angle.
117

118
119
120
121
122
123
124
125
        :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
126
                                                                          meshwidth=10,  # for faster processing
127
128
129
130
131
                                                                          subset=None,
                                                                          bandwise=0)
            else:
                # only a mean VAA is available
                if self.VAA_mean is None:
132
133
                    self.VAA_mean = \
                        GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
134
135
                    assert isinstance(self.VAA_mean, float)

136
                self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
137
138
139
140
141
        return self._VAA_arr

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

143
144
145
146
    @VAA_arr.deleter
    def VAA_arr(self):
        self._VAA_arr = None

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    @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],
167
168
                    accurracy=CFG.SZA_SAA_calculation_accurracy,
                    lonlat_arr=self.lonlat_arr if CFG.SZA_SAA_calculation_accurracy == 'fine' else None)
169
170
171
172
173
174
        return self._SZA_arr

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

175
176
177
178
    @SZA_arr.deleter
    def SZA_arr(self):
        self._SZA_arr = None

179
180
181
182
183
184
185
    @property
    def SAA_arr(self):
        """Get solar azimuth angle.

        :return:
        """
        if self._SAA_arr is None:
186
187
            # noinspection PyStatementEffect
            self.SZA_arr  # getter also sets self._SAA_arr
188
189
190
191
192
193
        return self._SAA_arr

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

194
195
196
197
    @SAA_arr.deleter
    def SAA_arr(self):
        self._SAA_arr = None

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    @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
213

214
215
216
217
    @RAA_arr.deleter
    def RAA_arr(self):
        self._RAA_arr = None

218
    def delete_ac_input_arrays(self):
219
220
221
222
223
224
225
        """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
226
227
228
229
230

        # 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
231
232
233


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

        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
        """
242
        # FIXME not yet usable for data < 2012 due to missing ECMWF archive
243
244
245
        L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)

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

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

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

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

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

274
275
276
277
278
    @property
    def logger(self):
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
279
            if len(self.inObjs) == 1:
280
281
282
283
284
285
286
287
288
289
                # 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
290
                    fileHandler = logging.FileHandler(path_logfile, mode='a')
291
                    fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
292
                    fileHandler.setLevel(CFG.log_level)
293
294
295

                    logger_atmCorr.addHandler(fileHandler)

Daniel Scheffler's avatar
Daniel Scheffler committed
296
297
                    inObj.close_GMS_loggers()

298
299
300
301
302
303
            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
304
            "AtmCorr.logger can not be set to %s." % logger
305
306
307
308
309
310
311
312
        if logger in ['not set', None]:
            self._logger.close()
            self._logger = logger
        else:
            self._logger = logger

    @logger.deleter
    def logger(self):
313
314
315
        if self._logger not in [None, 'not set']:
            self._logger.close()
            self._logger = None
316

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

319
320
321
322
323
324
325
326
    @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
327
                              (' (%s)' % obj.subsystem if obj.subsystem else '') +
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
                              '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:
354
355
            data_dict = {}

356
            for inObj in self.inObjs:
357
                for bandN, bandIdx in inObj.arr.bandnames.items():
358
                    if bandN not in data_dict:
Daniel Scheffler's avatar
Daniel Scheffler committed
359
360
361
362
                        # 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)
363
                    else:
364
                        inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
Daniel Scheffler's avatar
Daniel Scheffler committed
365
                                             "exists multiple times." % bandN)
366

367
            # validate: data must have all bands needed for AC
Daniel Scheffler's avatar
Daniel Scheffler committed
368
369
            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]
370
371
            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
372
                                   'Received: %s' % (str(all_bNs_AC), str(list(sorted(data_dict.keys())))))
373
374
375

            self._data = data_dict

376
377
378
379
380
381
382
383
        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

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
    @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):
403
        """Returns S2A tile name.
404
        NOTE: this is only needed if no DEM is passed to ac_gms
405
406
407
408
409

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

Daniel Scheffler's avatar
Daniel Scheffler committed
410
        return ''  # FIXME
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
436
437
438

    @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

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
466
467
468
    @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
469

470
        if not self._metadata:
Daniel Scheffler's avatar
Daniel Scheffler committed
471
            del self.logger  # otherwise each input object would have multiple fileHandlers
472

Daniel Scheffler's avatar
Daniel Scheffler committed
473
474
475
476
477
478
479
480
481
482
483
484
485
            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()
            )
486
487

            self._metadata = metadata
488
489
490

        return self._metadata

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

504
    def _meta_get_spatial_samplings(self):
505
506
507
        """

        :return:
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
         {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}}
526
        """
527
528
        # set corner coordinates and dims
        spatial_samplings = {}
529
530
531

        for inObj in self.inObjs:

532
533
534
535
536
            # 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.')
537

538
539
            # set spatial information
            spatial_samplings[inObj.arr.xgsd] = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
540
541
542
543
544
545
                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)
546

547
548
549
        return spatial_samplings

    def _meta_get_solar_irradiance(self):
550
551
552
        """

        :return:
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
581
582
583
        {'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
584
        for inObj in self.inObjs:  # type: L1C_object
585
            for bandN, bandIdx in inObj.arr.bandnames.items():
586
                if bandN not in viewing_zenith:
587
588
                    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
589
                    # viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
590
591
592
593
594
595
596
597
598
599
        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
600
        for inObj in self.inObjs:  # type: L1C_object
601
            for bandN, bandIdx in inObj.arr.bandnames.items():
602
                if bandN not in viewing_azimuth:
Daniel Scheffler's avatar
Daniel Scheffler committed
603
                    arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim == 3 else inObj.VAA_arr
604
                    viewing_azimuth[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
605
                    # viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
606

607
608
609
610
611
612
        return viewing_azimuth

    def _meta_get_relative_viewing_azimuth(self):
        """

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

615
616
        relative_viewing_azimuth = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
617
        for inObj in self.inObjs:  # type: L1C_object
618
            for bandN, bandIdx in inObj.arr.bandnames.items():
619
                if bandN not in relative_viewing_azimuth:
620
621
                    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
622
623
                    # relative_viewing_azimuth[bandN] = \
                    #     inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr
624

625
        return relative_viewing_azimuth
626

627
628
629
630
631
632
633
634
    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
635
636
            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°
637
            # FIXME correct to reduce resolution here by factor 10?
638
639
640
641
642
643
644
645
646
647
648
649
        )

        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
650
            inObj4dem = [obj for obj in self.inObjs if obj.arr.xgsd == 20]
651
652
653
            if not inObj4dem:
                self.logger.warning('Sentinel-2 20m subsystem could not be found. DEM passed to '
                                    'atmospheric correction might have wrong resolution.')
654
655
656
657
            inObj4dem = inObj4dem[0]
        else:
            inObj4dem = self.inObjs[0]

658
659
660
661
        try:
            dem = inObj4dem.dem[:].astype(np.float32)
        except Exception as e:
            dem = None
Daniel Scheffler's avatar
Daniel Scheffler committed
662
            self.logger.warning('A static elevation is assumed during atmospheric correction due to an error during '
663
664
665
                                '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:")
666
            self.logger.warning(traceback.format_exc())
667
668

        return dem
669
670

    def _get_srf(self):
671
        """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
672
        """
673
674
675
        # 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
676
        GMS_identifier_fullScene['Subsystem'] = ''
677
678
679
        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)
680

681
682
683
684
685
686
    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:
        """

687
688
        tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']

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

        # compute cloud mask if not already provided
694
        if no_avail_CMs:
695
            algorithm = CFG.cloud_masking_algorithm[self.inObjs[0].satellite]
696

697
698
            if algorithm == 'SICOR':
                return None
699

700
701
702
703
704
            else:
                # FMASK or Classical Bayesian
                try:
                    from .cloud_masking import Cloud_Mask_Creator

705
                    CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm, tempdir_root=CFG.path_tempdir)
706
707
708
709
                    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
710
                except Exception:
711
712
                    self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
                    self.logger.error(traceback.format_exc())
713
                    return None
714

715
716
        else:
            # check if there is a cloud mask with suitable GSD
Daniel Scheffler's avatar
Daniel Scheffler committed
717
            inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd == tgt_res]
718
719
            if not inObjs2use:
                raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
Daniel Scheffler's avatar
Daniel Scheffler committed
720
                                 'GMS object provides a cloud mask with spatial resolution of %s.' % tgt_res)
721
722
723
724
725
726
727
728
            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)
729
            #    {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40}  # FIXME hardcoded
730
731
732
733
734
735

            # 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.')
736
737
738
739

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

740
741
        # get nodata value
        self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
742

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

750
751
752
        return S2Mask(mask_array=cm_array,
                      mask_legend=cm_legend,
                      geo_coding=cm_geocoding)
753

754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
    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"]

771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
        # NOTE: use_signals must be set to false if run_request runs within a multiprocessing worker!
        @timeout_decorator.timeout(seconds=60*5, timeout_exception=TimeoutError, use_signals=False)
        def run_request():
            try:
                t0 = time()
                # TODO  Implement Lock to avoid too many parallel requests to ECMWF API. This might be an issue in
                # TODO  multiprocessing.
                results = download_variables(date_from=self.inObjs[0].acq_datetime,
                                             date_to=self.inObjs[0].acq_datetime,
                                             db_path=CFG.path_ECMWF_db,
                                             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 APIKeyFetchError:
                self.logger.error("ECMWF data download failed due to missing API credentials.")

            except (ECMWFAPIException, Exception):
                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())

798
        try:
799
800
801
            run_request()
        except TimeoutError:
            self.logger.error("ECMWF data download failed due to API request timeout after waiting 5 minutes.")
802
803
804
805

    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"]):
806
807
            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))
808
809
810
            # 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

811
812
    def run_atmospheric_correction(self, dump_ac_input=False):
        # type: (bool) -> list
813
814
815
        """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
        containing surface reflectance.

816
817
        :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
818
        """
819
820

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

823
        rs_data = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
824
825
826
827
828
829
830
831
832
            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
833
        )  # NOTE: all keys of this dict are later converted to attributes of RSImage
834

835
836
837
838
        # 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
839
        script = False
840

841
        # check if ECMWF data are available - if not, start the download
842
        if CFG.auto_download_ecmwf:
843
            self._check_or_download_ECMWF_data()
844
845

        # validate SNR
846
847
        if CFG.ac_estimate_accuracy:
            self._validate_snr_source()
848

849
850
851
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

852
        self.ac_input = dict(
853
            rs_image=rs_image,
Daniel Scheffler's avatar
Daniel Scheffler committed
854
            options=self.options,  # type: dict
855
856
            logger=repr(self.logger),  # only a string
            script=script
857
        )
858

859
860
861
862
        # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
        # with open(path_dump, 'wb') as outF:
        #     dill.dump(self.ac_input, outF)

863
        # run AC
864
        self.logger.info('Atmospheric correction started.')
865
        try:
866
            rs_image.logger = self.logger
867
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
868

869
        except Exception as e:
870
871
872
873
874
875
            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'
                              % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
            self.logger.error(traceback.format_exc())
            # TODO include that in the job summary

876
            # serialialize AC input
877
878
879
880
881
882
            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
883
                                  % path_dump)
884
885

            # delete AC input arrays
Daniel Scheffler's avatar
Daniel Scheffler committed
886
            for inObj in self.inObjs:  # type: L1C_object
887
888
                inObj.delete_ac_input_arrays()

889
890
            return list(self.inObjs)

891
        # get processing infos
892
        self.proc_info = self.ac_input['options']['processing']
893

894
        # join results
Daniel Scheffler's avatar
Daniel Scheffler committed
895
        self._join_results_to_inObjs()  # sets self.outObjs
896

897
898
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
899

900
901
902
        return self.outObjs

    def _join_results_to_inObjs(self):
903
904
905
        """
        Join results of atmospheric correction to the input GMS objects.
        """
906

907
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
Daniel Scheffler's avatar
Daniel Scheffler committed
908
909
910
        # delete logger
        # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
        del self.logger
911
912

        self._join_data_ac()
913
        self._join_data_errors(bandwise=CFG.ac_bandwise_accuracy)
914
915
916
917
918
        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]
919

920
921
922
        # update AC processing info
        [inObj.ac_options['processing'].update(self.proc_info) for inObj in self.inObjs]

923
924
925
926
        self.outObjs = self.inObjs

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

931
        if self.results.data_ac is not None:
932
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
933
934
                self.logger.info('Joining image data to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
935

936
                assert isinstance(inObj, L1B_object)
937
                nodata = self.results.nodata[inObj.arr.xgsd]  # 2D mask with True outside of image coverage
Daniel Scheffler's avatar
Daniel Scheffler committed
938
                ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
939
                out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
940

941
942
943
                # update metadata
                inObj.arr_desc = 'BOA_Ref'
                inObj.MetaObj.bands = len(self.results.data_ac)
944
                inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
945
946
                inObj.MetaObj.LayerBandsAssignment = out_LBA
                inObj.MetaObj.filter_layerdependent_metadata()
947
                inObj.meta_odict = inObj.MetaObj.to_odict()  # actually auto-updated by getter
948

949
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config
950
951
                # FIXME AC output nodata values = 0 -> new nodata areas but mask not updated
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
952
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
953
                surf_refl *= CFG.scale_factor_BOARef  # scale using scale factor (output is float16)
954
955
                # 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
956
                # apply the original nodata mask (indicating background values)
957
                surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
958

Daniel Scheffler's avatar
Daniel Scheffler committed
959
                if self.results.bad_data_value is np.nan:
960
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
961
                else:
Daniel Scheffler's avatar
Daniel Scheffler committed
962
963
                    surf_refl[
                        surf_refl == self.results.bad_data_value] = oF_refl  # FIXME meaningful to set AC nans to -9999?
964
965
966

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

969
970
971
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
972

973
974
    def _join_data_errors(self, bandwise=False):
        """Join ERRORS ARRAY as 3D or 2D int8 or int16 BOA reflectance array, scaled to scale factor from config.
975

976
977
978
        :param bandwise:    if True, a 3D array with bandwise information for each pixel is generated
        :return:
        """
979
980
        # TODO ac_error values are very close to 0 -> a scale factor of 255 yields int8 values below 10
        # TODO => better to stretch the whole array to values between 0 and 100 and save original min/max?
981
        if self.results.data_errors is not None:
982

983
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
984
985
                inObj.logger.info('Joining AC errors to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
986

987
988
                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()]
989
990
                out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
                out_nodata_val = get_outFillZeroSaturated(out_dtype)[0]
991

992
                # generate raw ac_errors array
993
                ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
994

995
996
997
998
999
1000
1001
                # apply scale factor from config to data pixels and overwrite nodata area with nodata value
                ac_errors *= CFG.ac_scale_factor_errors  # scale using scale factor (output is float16)
                ac_errors[np.isnan(ac_errors)] = out_nodata_val  # replace NaNs with outnodata
                ac_errors[nodata] = out_nodata_val  # set all positions where SICOR nodata mask is 1 to outnodata
                ac_errors = np.around(ac_errors).astype(out_dtype)  # round floats to next even int8/int16 value

                # average array over bands if no bandwise information is requested
1002
                if not bandwise:
1003
1004
                    # in case of only one subsystem: directly compute median errors here
                    if len(self.inObjs) == 1:
1005
                        ac_errors = np.median(ac_errors, axis=2).astype(ac_errors.dtype)
1006
1007
1008
1009
1010

                    # in case of multiple subsystems: dont compute median here but first apply geometric homogenization
                    # -> median could only be computed for each subsystem separately
                    else:
                        pass
1011
1012

                # set inObj.ac_errors
1013
                inObj.ac_errors = ac_errors  # setter generates a GeoArray with the same bandnames like inObj.arr
1014

1015
        elif not CFG.ac_estimate_accuracy:
1016
1017
            self.logger.info("Atmospheric correction did not provide a 'data_errors' array because "
                             "'ac_estimate_accuracy' was set to False in the job configuration.")
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
        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
1030

1031
1032
            joined = False
            for inObj in self.inObjs:
1033

1034
1035
                # delete all previous cloud masks
                del inObj.mask_clouds
1036
1037
1038

                # append mask_clouds only to the input GMS object with the same dimensions
                if inObj.arr.shape[:2] == mask_clouds_ac.shape:
1039
                    inObj.logger.info('Joining mask_clouds to %s.' % (inObj.entity_ID if not inObj.subsystem else
Daniel Scheffler's avatar
Daniel Scheffler committed
1040
                                                                      '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1041

1042
                    inObj.mask_clouds = mask_clouds_ac
1043
1044
                    inObj.mask_clouds.legend = self.results.mask_clouds.mask_legend  # dict(value=string, string=value))
                    # FIXME legend is not used later
1045
1046

                    # set cloud mask nodata value
1047
                    tgt_nodata = get_outFillZeroSaturated(mask_clouds_ac.dtype)[0]
1048
1049
                    ac_out_nodata = self.ac_input['options']['cld_mask']['nodata_value_mask']
                    if tgt_nodata not in self.results.mask_clouds.mask_legend.keys():
1050
                        inObj.mask_clouds[inObj.mask_clouds[:] == ac_out_nodata] = tgt_nodata
1051
1052
1053
1054
                        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.'
1055
                                      % (tgt_nodata, self.results.mask_clouds.mask_legend[tgt_nodata], ac_out_nodata))
1056
1057
1058
1059
                        mask_clouds_nodata = ac_out_nodata

                    inObj.mask_clouds.nodata = mask_clouds_nodata

1060
                    joined = True
1061

1062
1063
1064
1065
            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.')

1066
        else:
1067
1068
            self.logger.warning("Atmospheric correction did not provide a 'mask_clouds.mask_array' array. "
                                "GMS_object.mask_clouds kept None.")