L1C_P.py 50.4 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2
3
4

# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data
#
5
# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6
7
8
9
10
11
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
12
13
14
15
16
17
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version. Please note the following exception: `spechomo` depends on tqdm, which
# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
18
19
20
21
22
23
24
25
26
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

Daniel Scheffler's avatar
Daniel Scheffler committed
27
"""Level 1C Processor:   Atmospheric correction of TOA-reflectance data."""
Daniel Scheffler's avatar
Daniel Scheffler committed
28

29
import warnings
30
31
import re
import logging
32
import dill
33
import traceback
Daniel Scheffler's avatar
Daniel Scheffler committed
34
from typing import List  # noqa F401  # flake8 issue
35
36
from time import time
import os
37
import timeout_decorator
38

39
import numpy as np
40
41
from ecmwfapi.api import APIKeyFetchError
from ecmwfapi.api import APIException as ECMWFAPIException
Daniel Scheffler's avatar
Daniel Scheffler committed
42

43
from geoarray import GeoArray
44
from py_tools_ds.geo.map_info import mapinfo2geotransform
45
from pyrsr import RSR
46

47
from ..options.config import GMS_config as CFG
48
from . import geoprocessing as GEOP
Daniel Scheffler's avatar
Daniel Scheffler committed
49
from .L1B_P import L1B_object
50
from ..model.metadata import get_LayerBandsAssignment
51
from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
52
from ..misc.locks import MultiSlotLock
53
# from .cloud_masking import Cloud_Mask_Creator  # circular dependencies
54

55
from sicor.sicor_ac import ac_gms
56
from sicor.sensors import RSImage
57
from sicor.Mask import S2Mask
58
from sicor.ECMWF import download_variables
59

Daniel Scheffler's avatar
Daniel Scheffler committed
60
61
__author__ = 'Daniel Scheffler'

Daniel Scheffler's avatar
Daniel Scheffler committed
62

63
class L1C_object(L1B_object):
64
    def __init__(self, L1B_obj=None):
65
        super(L1C_object, self).__init__()
66
67
68

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

71
72
73
74
75
76
77
78
        # 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

79
        self.proc_level = 'L1C'
80
        self.proc_status = 'initialized'
81

82
83
84
    @property
    def lonlat_arr(self):
        """Calculates pixelwise 2D-array with longitude and latitude coordinates.
85

86
87
88
89
90
91
        :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,
92
93
                                            mapinfo2geotransform(self.MetaObj.map_info),
                                            self.MetaObj.projection,
Daniel Scheffler's avatar
Daniel Scheffler committed
94
95
                                            meshwidth=10,  # for faster processing
                                            nodata_mask=None,  # dont overwrite areas outside the image with nodata
96
97
                                            outFill=get_outFillZeroSaturated(np.float32)[0])[0]
        return self._lonlat_arr
98

99
100
101
    @lonlat_arr.setter
    def lonlat_arr(self, lonlat_arr):
        self._lonlat_arr = lonlat_arr
102

103
104
105
106
    @lonlat_arr.deleter
    def lonlat_arr(self):
        self._lonlat_arr = None

107
108
109
110
111
112
113
114
    @property
    def VZA_arr(self):
        """Get viewing zenith angle.

        :return:
        """
        if self._VZA_arr is None:
            self.logger.info('Calculating viewing zenith array...')
115
            if self.MetaObj.ViewingAngle_arrProv:
116
                # Sentinel-2
117
118
119
120
121
122
                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)
123
124
            else:
                self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
125
126
                                                    float(self.MetaObj.ViewingAngle),
                                                    float(self.MetaObj.FOV),
127
                                                    self.logger,
Daniel Scheffler's avatar
Daniel Scheffler committed
128
                                                    nodata_mask=None,  # dont overwrite areas outside image with nodata
129
                                                    outFill=get_outFillZeroSaturated(np.float32)[0],
Daniel Scheffler's avatar
Daniel Scheffler committed
130
                                                    meshwidth=10)  # for faster processing
131
132
133
134
135
        return self._VZA_arr

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

137
138
139
140
    @VZA_arr.deleter
    def VZA_arr(self):
        self._VZA_arr = None

141
142
143
    @property
    def VAA_arr(self):
        """Get viewing azimuth angle.
144

145
146
147
148
        :return:
        """
        if self._VAA_arr is None:
            self.logger.info('Calculating viewing azimuth array...')
149
            if self.MetaObj.IncidenceAngle_arrProv:
150
                # Sentinel-2
151
152
153
154
155
156
                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)
157
158
159
            else:
                # only a mean VAA is available
                if self.VAA_mean is None:
160
161
                    self.VAA_mean = \
                        GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
162
163
                    assert isinstance(self.VAA_mean, float)

164
                self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
165
166
167
168
169
        return self._VAA_arr

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

171
172
173
174
    @VAA_arr.deleter
    def VAA_arr(self):
        self._VAA_arr = None

175
176
177
178
179
180
181
182
183
184
185
    @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,
186
187
                    self.MetaObj.AcqDate,
                    self.MetaObj.AcqTime,
188
189
                    self.fullSceneCornerPos,
                    self.fullSceneCornerLonLat,
190
                    self.MetaObj.overpassDurationSec,
191
192
193
194
                    self.logger,
                    meshwidth=10,
                    nodata_mask=None,  # dont overwrite areas outside the image with nodata
                    outFill=get_outFillZeroSaturated(np.float32)[0],
195
196
                    accurracy=CFG.SZA_SAA_calculation_accurracy,
                    lonlat_arr=self.lonlat_arr if CFG.SZA_SAA_calculation_accurracy == 'fine' else None)
197
198
199
200
201
202
        return self._SZA_arr

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

203
204
205
206
    @SZA_arr.deleter
    def SZA_arr(self):
        self._SZA_arr = None

207
208
209
210
211
212
213
    @property
    def SAA_arr(self):
        """Get solar azimuth angle.

        :return:
        """
        if self._SAA_arr is None:
214
215
            # noinspection PyStatementEffect
            self.SZA_arr  # getter also sets self._SAA_arr
216
217
218
219
220
221
        return self._SAA_arr

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

222
223
224
225
    @SAA_arr.deleter
    def SAA_arr(self):
        self._SAA_arr = None

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
    @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
241

242
243
244
245
    @RAA_arr.deleter
    def RAA_arr(self):
        self._RAA_arr = None

246
    def delete_ac_input_arrays(self):
247
248
249
250
251
252
253
        """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
254
255
256
257
258

        # 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
259
260
261


class AtmCorr(object):
262
    def __init__(self, *L1C_objs, reporting=False):
263
        """Wrapper around atmospheric correction by Andre Hollstein, GFZ Potsdam
264
265
266
267
268
269

        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
        """
270
        # FIXME not yet usable for data < 2012 due to missing ECMWF archive
271
272
273
        L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)

        # hidden attributes
Daniel Scheffler's avatar
Daniel Scheffler committed
274
275
276
        self._logger = None
        self._GSDs = []
        self._data = {}
277
        self._metadata = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
278
        self._nodata = {}
279
        self._band_spatial_sampling = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
280
        self._options = {}
281
282
283

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

Daniel Scheffler's avatar
Daniel Scheffler committed
287
        self.inObjs = L1C_objs  # type: List[L1C_object]
288
        self.reporting = reporting
Daniel Scheffler's avatar
Daniel Scheffler committed
289
290
        self.ac_input = {}  # set by self.run_atmospheric_correction()
        self.results = None  # direct output of external atmCorr module (set by run_atmospheric_correction)
291
        self.proc_info = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
292
        self.outObjs = []  # atmospherically corrected L1C objects
293
294

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

297
        if not re.search(r'Sentinel-2', self.inObjs[0].satellite, re.I):
298
299
            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
300

301
302
303
304
305
    @property
    def logger(self):
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
306
            if len(self.inObjs) == 1:
307
308
309
310
311
312
313
314
315
316
                # 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
317
                    fileHandler = logging.FileHandler(path_logfile, mode='a')
318
                    fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
319
                    fileHandler.setLevel(CFG.log_level)
320
321
322

                    logger_atmCorr.addHandler(fileHandler)

323
                    inObj.close_loggers()
Daniel Scheffler's avatar
Daniel Scheffler committed
324

325
326
327
328
329
330
            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
331
            "AtmCorr.logger can not be set to %s." % logger
332
333
334
335
336
337
338
339
        if logger in ['not set', None]:
            self._logger.close()
            self._logger = logger
        else:
            self._logger = logger

    @logger.deleter
    def logger(self):
340
341
342
        if self._logger not in [None, 'not set']:
            self._logger.close()
            self._logger = None
343

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

346
347
348
349
350
351
352
353
    @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
354
                              (' (%s)' % obj.subsystem if obj.subsystem else '') +
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
                              '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:
381
382
            data_dict = {}

383
            for inObj in self.inObjs:
384
                for bandN, bandIdx in inObj.arr.bandnames.items():
385
                    if bandN not in data_dict:
Daniel Scheffler's avatar
Daniel Scheffler committed
386
387
388
                        # 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
389
                        data_dict[bandN] = (arr2pass / inObj.MetaObj.ScaleFactor).astype(np.float16)
390
                    else:
391
                        inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
Daniel Scheffler's avatar
Daniel Scheffler committed
392
                                             "exists multiple times." % bandN)
393

394
            # validate: data must have all bands needed for AC
Daniel Scheffler's avatar
Daniel Scheffler committed
395
396
            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]
397
398
            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
399
                                   'Received: %s' % (str(all_bNs_AC), str(list(sorted(data_dict.keys())))))
400
401
402

            self._data = data_dict

403
404
405
406
407
408
409
410
        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

411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
    @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):
430
        """Returns S2A tile name.
431
        NOTE: this is only needed if no DEM is passed to ac_gms
432
433
434
435
436

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

Daniel Scheffler's avatar
Daniel Scheffler committed
437
        return ''  # FIXME
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465

    @property
    def 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

466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
    @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
496

497
        if not self._metadata:
Daniel Scheffler's avatar
Daniel Scheffler committed
498
            del self.logger  # otherwise each input object would have multiple fileHandlers
499

Daniel Scheffler's avatar
Daniel Scheffler committed
500
            metadata = dict(
501
                U=self.inObjs[0].MetaObj.EarthSunDist,
Daniel Scheffler's avatar
Daniel Scheffler committed
502
503
504
505
506
                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(),
507
508
                sun_mean_azimuth=self.inObjs[0].MetaObj.SunAzimuth,
                sun_mean_zenith=90 - self.inObjs[0].MetaObj.SunElevation,
Daniel Scheffler's avatar
Daniel Scheffler committed
509
510
511
512
                solar_irradiance=self._meta_get_solar_irradiance(),
                aux_data=self._meta_get_aux_data(),
                spatial_samplings=self._meta_get_spatial_samplings()
            )
513
514

            self._metadata = metadata
515
516
517

        return self._metadata

518
519
    @property
    def options(self):
520
        # type: () -> dict
521
522
523
524
525
526
        """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
527
            self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
528
            self._options["report"]["reporting"] = self.reporting
529
530
            return self._options

531
    def _meta_get_spatial_samplings(self):
532
533
534
        """

        :return:
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
         {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}}
553
        """
554
555
        # set corner coordinates and dims
        spatial_samplings = {}
556
557
558

        for inObj in self.inObjs:

559
560
561
562
563
            # 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.')
564

565
566
            # set spatial information
            spatial_samplings[inObj.arr.xgsd] = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
567
568
569
570
571
572
                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)
573

574
575
576
        return spatial_samplings

    def _meta_get_solar_irradiance(self):
577
578
579
        """

        :return:
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
        {'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:
598
599
            for bandN in inObj.arr.bandnames:
                lba_key = bandN[2:] if bandN.startswith('B0') else bandN[1:]
600
                if bandN not in solar_irradiance:
601
602
                    solar_irradiance[bandN] = inObj.MetaObj.SolIrradiance[lba_key]

603
604
605
606
607
608
609
610
611
612
        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
613
        for inObj in self.inObjs:  # type: L1C_object
614
            for bandN, bandIdx in inObj.arr.bandnames.items():
615
                if bandN not in viewing_zenith:
616
617
                    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
618
                    # viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
619
620
621
622
623
624
625
626
627
628
        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
629
        for inObj in self.inObjs:  # type: L1C_object
630
            for bandN, bandIdx in inObj.arr.bandnames.items():
631
                if bandN not in viewing_azimuth:
Daniel Scheffler's avatar
Daniel Scheffler committed
632
                    arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim == 3 else inObj.VAA_arr
633
                    viewing_azimuth[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
634
                    # viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
635

636
637
638
639
640
641
        return viewing_azimuth

    def _meta_get_relative_viewing_azimuth(self):
        """

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

644
645
        relative_viewing_azimuth = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
646
        for inObj in self.inObjs:  # type: L1C_object
647
            for bandN, bandIdx in inObj.arr.bandnames.items():
648
                if bandN not in relative_viewing_azimuth:
649
650
                    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
651
652
                    # relative_viewing_azimuth[bandN] = \
                    #     inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr
653

654
        return relative_viewing_azimuth
655

656
657
658
659
660
661
662
663
    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
664
665
            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°
666
            # FIXME correct to reduce resolution here by factor 10?
667
668
669
670
671
672
673
674
675
676
        )

        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
677
        if re.search(r'Sentinel-2', self.inObjs[0].satellite):
678
            # in case of Sentinel-2 the 20m DEM must be passed
Daniel Scheffler's avatar
Daniel Scheffler committed
679
            inObj4dem = [obj for obj in self.inObjs if obj.arr.xgsd == 20]
680
681
682
            if not inObj4dem:
                self.logger.warning('Sentinel-2 20m subsystem could not be found. DEM passed to '
                                    'atmospheric correction might have wrong resolution.')
683
684
685
686
            inObj4dem = inObj4dem[0]
        else:
            inObj4dem = self.inObjs[0]

687
688
689
690
        try:
            dem = inObj4dem.dem[:].astype(np.float32)
        except Exception as e:
            dem = None
Daniel Scheffler's avatar
Daniel Scheffler committed
691
            self.logger.warning('A static elevation is assumed during atmospheric correction due to an error during '
692
693
694
                                '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:")
695
            self.logger.warning(traceback.format_exc())
696
697

        return dem
698
699

    def _get_srf(self):
700
        """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
701
        """
702
        # FIXME calculation of center wavelengths within RSR() used not the GMS algorithm
703
        # SRF instance must be created for all bands and the previous proc level
704
705
706
707
708
709
710
711
712
713
714
715
716
717
        GMSid_fullScene = self.inObjs[0].GMS_identifier
        GMSid_fullScene.subsystem = ''
        GMSid_fullScene.proc_level = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]

        return RSR(satellite=GMSid_fullScene.satellite,
                   sensor=GMSid_fullScene.sensor,
                   subsystem=GMSid_fullScene.subsystem,
                   wvl_unit='nanometers',
                   format_bandnames=True,
                   no_pan=CFG.skip_pan,
                   no_thermal=CFG.skip_thermal,
                   after_ac=GMSid_fullScene.proc_level in ['L1C', 'L2A', 'L2B', 'L2C'],
                   sort_by_cwl=CFG.sort_bands_by_cwl
                   )
718

719
720
721
722
723
724
    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:
        """

725
726
        tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']

727
        # check if input GMS objects provide a cloud mask
728
        avail_cloud_masks = {inObj.GMS_identifier.subsystem: inObj.mask_clouds for inObj in self.inObjs}
729
        no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
730
731

        # compute cloud mask if not already provided
732
        if no_avail_CMs:
733
            algorithm = CFG.cloud_masking_algorithm[self.inObjs[0].satellite]
734

735
736
            if algorithm == 'SICOR':
                return None
737

738
739
740
741
742
            else:
                # FMASK or Classical Bayesian
                try:
                    from .cloud_masking import Cloud_Mask_Creator

743
                    CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm, tempdir_root=CFG.path_tempdir)
744
745
746
747
                    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
748
                except Exception:
749
750
                    self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
                    self.logger.error(traceback.format_exc())
751
                    return None
752

753
754
        else:
            # check if there is a cloud mask with suitable GSD
Daniel Scheffler's avatar
Daniel Scheffler committed
755
            inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd == tgt_res]
756
757
            if not inObjs2use:
                raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
Daniel Scheffler's avatar
Daniel Scheffler committed
758
                                 'GMS object provides a cloud mask with spatial resolution of %s.' % tgt_res)
759
760
761
762
763
764
765
766
            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)
767
            #    {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40}  # FIXME hardcoded
768
769
770
771
772
773

            # 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.')
774
775
776
777

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

778
779
        # get nodata value
        self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
780

781
        # append cloud mask to input object with the same spatial resolution if there was no mask before
782
        for inObj in self.inObjs:
783
            if inObj.arr.xgsd == cm_geoarray.xgsd:
784
785
                inObj.mask_clouds = cm_geoarray
                inObj.build_combined_masks_array()
Daniel Scheffler's avatar
Daniel Scheffler committed
786
                break  # appending it to one inObj is enough
787

788
789
790
        return S2Mask(mask_array=cm_array,
                      mask_legend=cm_legend,
                      geo_coding=cm_geocoding)
791

792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
    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"]

809
810
        # 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)
811
812
        def run_request():
            try:
813
814
815
816
817
818
819
820
821
822
823
824
825
826
                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)
827
828
829
830
831
832
833
834
835

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

836
        try:
837
838
839
            run_request()
        except TimeoutError:
            self.logger.error("ECMWF data download failed due to API request timeout after waiting 5 minutes.")
840
841
842
843

    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"]):
844
845
            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))
846
847
848
            # 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

849
850
    def run_atmospheric_correction(self, dump_ac_input=False):
        # type: (bool) -> list
851
852
853
        """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
        containing surface reflectance.

854
855
        :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
856
        """
857
858

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

861
        rs_data = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
862
863
864
865
866
867
868
869
870
            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
871
        )  # NOTE: all keys of this dict are later converted to attributes of RSImage
872

873
874
875
876
        # 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
877
        script = False
878

879
        # check if ECMWF data are available - if not, start the download
880
        if CFG.auto_download_ecmwf:
881
            self._check_or_download_ECMWF_data()
882
883

        # validate SNR
884
885
        if CFG.ac_estimate_accuracy:
            self._validate_snr_source()
886

887
888
889
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

890
        self.ac_input = dict(
891
            rs_image=rs_image,
Daniel Scheffler's avatar
Daniel Scheffler committed
892
            options=self.options,  # type: dict
893
894
            logger=repr(self.logger),  # only a string
            script=script
895
        )
896

897
898
899
900
        # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
        # with open(path_dump, 'wb') as outF:
        #     dill.dump(self.ac_input, outF)

901
        # run AC
902
        self.logger.info('Atmospheric correction started.')
903
        try:
904
            rs_image.logger = self.logger
905
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
906

907
        except Exception as e:
908
909
910
911
912
913
            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

914
            # serialialize AC input
915
916
917
918
919
920
            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
921
                                  % path_dump)
922
923

            # delete AC input arrays
Daniel Scheffler's avatar
Daniel Scheffler committed
924
            for inObj in self.inObjs:  # type: L1C_object
925
926
                inObj.delete_ac_input_arrays()

927
928
            return list(self.inObjs)

929
930
931
932
        finally:
            # rs_image.logger must be closed properly in any case
            rs_image.logger.close()

933
        # get processing infos
934
        self.proc_info = self.ac_input['options']['processing']
935

936
        # join results
Daniel Scheffler's avatar
Daniel Scheffler committed
937
        self._join_results_to_inObjs()  # sets self.outObjs
938

939
940
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
941

942
943
944
        return self.outObjs

    def _join_results_to_inObjs(self):
945
946
947
        """
        Join results of atmospheric correction to the input GMS objects.
        """
948

949
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
Daniel Scheffler's avatar
Daniel Scheffler committed
950
951
952
        # delete logger
        # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
        del self.logger
953
954

        self._join_data_ac()
955
        self._join_data_errors(bandwise=CFG.ac_bandwise_accuracy)
956
957
958
959
960
        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]
961

962
963
964
        # update AC processing info
        [inObj.ac_options['processing'].update(self.proc_info) for inObj in self.inObjs]

965
966
967
968
        self.outObjs = self.inObjs

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

973
        if self.results.data_ac is not None:
974
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
975
976
                self.logger.info('Joining image data to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
977

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

983
984
985
                # update metadata #
                ###################

986
987
                inObj.arr_desc = 'BOA_Ref'
                inObj.MetaObj.bands = len(self.results.data_ac)
988
                inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
989
                inObj.MetaObj.LayerBandsAssignment = out_LBA
990
                inObj.LayerBandsAssignment = out_LBA
991
992
                inObj.MetaObj.filter_layerdependent_metadata()

993
994
995
996
                ##################################################################################
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config #
                ##################################################################################

997
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
998
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
999
                surf_refl *= CFG.scale_factor_BOARef  # scale using scale factor (output is float16)
1000
1001
1002
1003
1004
1005
1006
1007
1008

                # 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
1009

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

1013
1014
                # 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
1015
                if self.results.bad_data_value is np.nan:
1016
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
1017
                else:
1018
                    surf_refl[surf_refl == self.results.bad_data_value] = oF_refl
1019

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

1023
1024
1025
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
1026

1027
1028
    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.
1029

1030
1031
1032
        :param bandwise:    if True, a 3D array with bandwise information for each pixel is generated
        :return:
        """
1033
1034
        # 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?
1035
        if self.results.data_errors is not None:
1036

1037
            for inObj in self.inObjs:
Daniel Scheffler's avatar
Daniel Scheffler committed
1038
1039
                inObj.logger.info('Joining AC errors to %s.' % (inObj.entity_ID if not inObj.subsystem else
                                                                '%s %s' % (inObj.entity_ID, inObj.subsystem)))
1040

1041
1042
                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()]
1043
1044
                out_dtype = np.int8 if CFG.ac_scale_factor_errors <= 255 else np.int16
                out_nodata_val = get_outFillZeroSaturated(out_dtype)[0]
1045

1046
                # generate raw ac_errors array
1047
                ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
1048

1049
1050
1051
1052
1053
1054
1055
                # 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 outnodat