L1C_P.py 48.8 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 ..misc.locks import MultiSlotLock
27
from ..io.input_reader import SRF
28
# from .cloud_masking import Cloud_Mask_Creator  # circular dependencies
29

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
37

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

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

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

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

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

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

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

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

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...')
90
            if self.MetaObj.ViewingAngle_arrProv:
91
                # Sentinel-2
92
93
94
95
96
97
                self._VZA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
                    {k: v.tolist() for k, v in self.MetaObj.ViewingAngle_arrProv.items()},
                    self.shape_fullArr,
                    meshwidth=10,  # for faster processing
                    subset=None,
                    bandwise=0)
98
99
            else:
                self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
100
101
                                                    float(self.MetaObj.ViewingAngle),
                                                    float(self.MetaObj.FOV),
102
                                                    self.logger,
Daniel Scheffler's avatar
Daniel Scheffler committed
103
                                                    nodata_mask=None,  # dont overwrite areas outside image with nodata
104
                                                    outFill=get_outFillZeroSaturated(np.float32)[0],
Daniel Scheffler's avatar
Daniel Scheffler committed
105
                                                    meshwidth=10)  # for faster processing
106
107
108
109
110
        return self._VZA_arr

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

112
113
114
115
    @VZA_arr.deleter
    def VZA_arr(self):
        self._VZA_arr = None

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

120
121
122
123
        :return:
        """
        if self._VAA_arr is None:
            self.logger.info('Calculating viewing azimuth array...')
124
            if self.MetaObj.IncidenceAngle_arrProv:
125
                # Sentinel-2
126
127
128
129
130
131
                self._VAA_arr = GEOP.adjust_acquisArrProv_to_shapeFullArr(
                    {k: v.tolist() for k, v in self.MetaObj.IncidenceAngle_arrProv.items()},
                    self.shape_fullArr,
                    meshwidth=10,  # for faster processing
                    subset=None,
                    bandwise=0)
132
133
134
            else:
                # only a mean VAA is available
                if self.VAA_mean is None:
135
136
                    self.VAA_mean = \
                        GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
137
138
                    assert isinstance(self.VAA_mean, float)

139
                self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
140
141
142
143
144
        return self._VAA_arr

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

146
147
148
149
    @VAA_arr.deleter
    def VAA_arr(self):
        self._VAA_arr = None

150
151
152
153
154
155
156
157
158
159
160
    @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,
161
162
                    self.MetaObj.AcqDate,
                    self.MetaObj.AcqTime,
163
164
                    self.fullSceneCornerPos,
                    self.fullSceneCornerLonLat,
165
                    self.MetaObj.overpassDurationSec,
166
167
168
169
                    self.logger,
                    meshwidth=10,
                    nodata_mask=None,  # dont overwrite areas outside the image with nodata
                    outFill=get_outFillZeroSaturated(np.float32)[0],
170
171
                    accurracy=CFG.SZA_SAA_calculation_accurracy,
                    lonlat_arr=self.lonlat_arr if CFG.SZA_SAA_calculation_accurracy == 'fine' else None)
172
173
174
175
176
177
        return self._SZA_arr

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

178
179
180
181
    @SZA_arr.deleter
    def SZA_arr(self):
        self._SZA_arr = None

182
183
184
185
186
187
188
    @property
    def SAA_arr(self):
        """Get solar azimuth angle.

        :return:
        """
        if self._SAA_arr is None:
189
190
            # noinspection PyStatementEffect
            self.SZA_arr  # getter also sets self._SAA_arr
191
192
193
194
195
196
        return self._SAA_arr

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

197
198
199
200
    @SAA_arr.deleter
    def SAA_arr(self):
        self._SAA_arr = None

201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
    @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
216

217
218
219
220
    @RAA_arr.deleter
    def RAA_arr(self):
        self._RAA_arr = None

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

        # 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
234
235
236


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

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

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

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

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

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

272
        if not re.search('Sentinel-2', self.inObjs[0].satellite, re.I):
273
274
            self.logger.debug('Calculation of acquisition geometry arrays is currently only validated for Sentinel-2!')
            # validation possible by comparing S2 angles provided by ESA with own angles  # TODO
275

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

                    logger_atmCorr.addHandler(fileHandler)

298
                    inObj.close_loggers()
Daniel Scheffler's avatar
Daniel Scheffler committed
299

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

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

319
        [inObj.close_loggers() for inObj in self.inObjs]
Daniel Scheffler's avatar
Daniel Scheffler committed
320

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

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

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

            self._data = data_dict

378
379
380
381
382
383
384
385
        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

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
412
        return ''  # FIXME
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
439
440

    @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

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
469
470
    @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
471

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

Daniel Scheffler's avatar
Daniel Scheffler committed
475
            metadata = dict(
476
                U=self.inObjs[0].MetaObj.EarthSunDist,
Daniel Scheffler's avatar
Daniel Scheffler committed
477
478
479
480
481
                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(),
482
483
                sun_mean_azimuth=self.inObjs[0].MetaObj.SunAzimuth,
                sun_mean_zenith=90 - self.inObjs[0].MetaObj.SunElevation,
Daniel Scheffler's avatar
Daniel Scheffler committed
484
485
486
487
                solar_irradiance=self._meta_get_solar_irradiance(),
                aux_data=self._meta_get_aux_data(),
                spatial_samplings=self._meta_get_spatial_samplings()
            )
488
489

            self._metadata = metadata
490
491
492

        return self._metadata

493
494
    @property
    def options(self):
495
        # type: () -> dict
496
497
498
499
500
501
        """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
502
            self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
503
            self._options["report"]["reporting"] = self.reporting
504
505
            return self._options

506
    def _meta_get_spatial_samplings(self):
507
508
509
        """

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

        for inObj in self.inObjs:

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

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

549
550
551
        return spatial_samplings

    def _meta_get_solar_irradiance(self):
552
553
554
        """

        :return:
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
        {'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:
573
574
            for bandN in inObj.arr.bandnames:
                lba_key = bandN[2:] if bandN.startswith('B0') else bandN[1:]
575
                if bandN not in solar_irradiance:
576
577
                    solar_irradiance[bandN] = inObj.MetaObj.SolIrradiance[lba_key]

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

611
612
613
614
615
616
        return viewing_azimuth

    def _meta_get_relative_viewing_azimuth(self):
        """

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

619
620
        relative_viewing_azimuth = {}

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

629
        return relative_viewing_azimuth
630

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

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

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

        return dem
673
674

    def _get_srf(self):
675
        """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
676
        """
677
678
679
        # 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
Mathias Peters's avatar
Mathias Peters committed
680
681
        GMS_identifier_fullScene.subsystem = ''
        GMS_identifier_fullScene.proc_level = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]
682
683

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

685
686
687
688
689
690
    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:
        """

691
692
        tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']

693
        # check if input GMS objects provide a cloud mask
694
        avail_cloud_masks = {inObj.GMS_identifier.subsystem: inObj.mask_clouds for inObj in self.inObjs}
695
        no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
696
697

        # compute cloud mask if not already provided
698
        if no_avail_CMs:
699
            algorithm = CFG.cloud_masking_algorithm[self.inObjs[0].satellite]
700

701
702
            if algorithm == 'SICOR':
                return None
703

704
705
706
707
708
            else:
                # FMASK or Classical Bayesian
                try:
                    from .cloud_masking import Cloud_Mask_Creator

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

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

            # 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.')
740
741
742
743

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

744
745
        # get nodata value
        self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
746

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

754
755
756
        return S2Mask(mask_array=cm_array,
                      mask_legend=cm_legend,
                      geo_coding=cm_geocoding)
757

758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
    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"]

775
776
        # NOTE: use_signals must be set to True while executed as multiprocessing worker (e.g., in multiprocessing.Pool)
        @timeout_decorator.timeout(seconds=60*5, timeout_exception=TimeoutError)
777
778
        def run_request():
            try:
779
780
781
782
783
784
785
786
787
788
789
790
791
792
                with MultiSlotLock('ECMWF download lock', allowed_slots=1, logger=self.logger):
                    t0 = time()
                    # NOTE: download_variables does not accept a logger -> so the output may be invisible in WebApp
                    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)
793
794
795
796
797
798
799
800
801

            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())

802
        try:
803
804
805
            run_request()
        except TimeoutError:
            self.logger.error("ECMWF data download failed due to API request timeout after waiting 5 minutes.")
806
807
808
809

    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"]):
810
811
            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))
812
813
814
            # 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

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

820
821
        :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
822
        """
823
824

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

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

839
840
841
842
        # 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
843
        script = False
844

845
        # check if ECMWF data are available - if not, start the download
846
        if CFG.auto_download_ecmwf:
847
            self._check_or_download_ECMWF_data()
848
849

        # validate SNR
850
851
        if CFG.ac_estimate_accuracy:
            self._validate_snr_source()
852

853
854
855
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

856
        self.ac_input = dict(
857
            rs_image=rs_image,
Daniel Scheffler's avatar
Daniel Scheffler committed
858
            options=self.options,  # type: dict
859
860
            logger=repr(self.logger),  # only a string
            script=script
861
        )
862

863
864
865
866
        # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
        # with open(path_dump, 'wb') as outF:
        #     dill.dump(self.ac_input, outF)

867
        # run AC
868
        self.logger.info('Atmospheric correction started.')
869
        try:
870
            rs_image.logger = self.logger
871
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
872

873
        except Exception as e:
874
875
876
877
878
879
            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

880
            # serialialize AC input
881
882
883
884
885
886
            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
887
                                  % path_dump)
888
889

            # delete AC input arrays
Daniel Scheffler's avatar
Daniel Scheffler committed
890
            for inObj in self.inObjs:  # type: L1C_object
891
892
                inObj.delete_ac_input_arrays()

893
894
            return list(self.inObjs)

895
896
897
898
        finally:
            # rs_image.logger must be closed properly in any case
            rs_image.logger.close()

899
        # get processing infos
900
        self.proc_info = self.ac_input['options']['processing']
901

902
        # join results
Daniel Scheffler's avatar
Daniel Scheffler committed
903
        self._join_results_to_inObjs()  # sets self.outObjs
904

905
906
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
907

908
909
910
        return self.outObjs

    def _join_results_to_inObjs(self):
911
912
913
        """
        Join results of atmospheric correction to the input GMS objects.
        """
914

915
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
Daniel Scheffler's avatar
Daniel Scheffler committed
916
917
918
        # delete logger
        # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
        del self.logger
919
920

        self._join_data_ac()
921
        self._join_data_errors(bandwise=CFG.ac_bandwise_accuracy)
922
923
924
925
926
        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]
927

928
929
930
        # update AC processing info
        [inObj.ac_options['processing'].update(self.proc_info) for inObj in self.inObjs]

931
932
933
934
        self.outObjs = self.inObjs

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

939
        if self.results.data_ac is not None:
940
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
941
942
                self.logger.info('Joining image data to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
943

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

949
950
951
                # update metadata #
                ###################

952
953
                inObj.arr_desc = 'BOA_Ref'
                inObj.MetaObj.bands = len(self.results.data_ac)
954
                inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
955
                inObj.MetaObj.LayerBandsAssignment = out_LBA
956
                inObj.LayerBandsAssignment = out_LBA
957
958
                inObj.MetaObj.filter_layerdependent_metadata()

959
960
961
962
                ##################################################################################
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config #
                ##################################################################################

963
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
964
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
965
                surf_refl *= CFG.scale_factor_BOARef  # scale using scale factor (output is float16)
966
967
968
969
970
971
972
973
974

                # set AC nodata values to GMS outFill
                # NOTE: AC nodata contains a pixel mask where at least one band is no data
                #       => Setting these pixels to outZero would also reduce pixel values of surrounding pixels in
                #          spatial homogenization (because resampling only ignores -9999).
                #       It would be possible to generate a zero-data mask here for each subsystem and apply it after
                #       spatial homogenization. Alternatively zero-data pixels could be interpolated spectrally or
                #       spatially within L1A processor (also see issue #74).
                surf_refl[nodata] = oF_refl  # overwrite AC nodata values with GMS outFill
975

Daniel Scheffler's avatar
Daniel Scheffler committed
976
                # apply the original nodata mask (indicating background values)
977
                surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
978

979
980
                # set AC NaNs to GMS outFill
                # NOTE: SICOR result has NaNs at no data positions AND non-clear positions
Daniel Scheffler's avatar
Daniel Scheffler committed
981
                if self.results.bad_data_value is np.nan:
982
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
983
                else:
984
                    surf_refl[surf_refl == self.results.bad_data_value] = oF_refl
985

986
                # use inObj.arr setter to generate a GeoArray
987
                inObj.arr = surf_refl.astype(inObj.arr.dtype)  # -> int16 (also converts NaNs to 0 if needed
988

989
990
991
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
992

993
994
    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.
995

996
997
998
        :param bandwise:    if True, a 3D array with bandwise information for each pixel is generated
        :return:
        """
999
1000
        # 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?
1001
        if self.results.data_errors is not None:
1002

1003
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
1004
1005
                inObj.logger.info('Joining AC errors to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1006

1007
1008
                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()]
1009
1010
                out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
                out_nodata_val = get_outFillZeroSaturated(out_dtype)[0]
1011

1012
                # generate raw ac_errors array
1013
                ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
1014

1015
1016
1017
1018
1019
1020
1021
                # 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
1022
                if not bandwise:
1023
1024
                    # in case of only one subsystem: directly compute median errors here
                    if len(self.inObjs) == 1:
1025
                        ac_errors = np.median(ac_errors, axis=2).astype(ac_errors.dtype)
1026
1027
1028
1029
1030

                    # 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
1031
1032

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

1035
        elif not CFG.ac_estimate_accuracy:
1036
1037
            self.logger.info("Atmospheric correction did not provide a 'data_errors' array because "
                             "'ac_estimate_accuracy' was set to False in the job configuration.")
1038
1039
1040
1041
1042
1043
1044
1045