L1C_P.py 40.9 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
10

11
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
12

13
14
15
16
try:
    from osgeo import osr
except ImportError:
    import osr
Daniel Scheffler's avatar
Daniel Scheffler committed
17

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

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

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
35

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

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

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

52
        self.proc_level = 'L1C'
53

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

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

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

75
76
77
78
79
80
81
82
83
84
85
86
    @property
    def VZA_arr(self):
        """Get viewing zenith angle.

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

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

104
105
106
    @property
    def VAA_arr(self):
        """Get viewing azimuth angle.
107

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

126
                self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
127
128
129
130
131
        return self._VAA_arr

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

133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    @property
    def SZA_arr(self):
        """Get solar zenith angle.

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

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

    @property
    def SAA_arr(self):
        """Get solar azimuth angle.

        :return:
        """
        if self._SAA_arr is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
168
            _ = self.SZA_arr  # getter also sets self._SAA_arr
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
        return self._SAA_arr

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

    @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
190

191
    def delete_ac_input_arrays(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
192
193
194
195
196
        self.VZA_arr = None  # not needed anymore
        self.SZA_arr = None  # not needed anymore
        self.SAA_arr = None  # not needed anymore
        self.RAA_arr = None  # not needed anymore
        self.lonlat_arr = None  # not needed anymore
Daniel Scheffler's avatar
Daniel Scheffler committed
197
198
199
200
201

        # 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
202
203
204


class AtmCorr(object):
205
    def __init__(self, *L1C_objs, reporting=False):
206
        """Wrapper around atmospheric correction by Andre Hollstein, GFZ Potsdam
207
208
209
210
211
212

        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
        """
213
        # FIXME not yet usable for data < 2012 due to missing ECMWF archive
214
215
216
        L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)

        # hidden attributes
Daniel Scheffler's avatar
Daniel Scheffler committed
217
218
219
        self._logger = None
        self._GSDs = []
        self._data = {}
220
        self._metadata = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
221
        self._nodata = {}
222
        self._band_spatial_sampling = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
223
        self._options = {}
224
225
226

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

Daniel Scheffler's avatar
Daniel Scheffler committed
230
        self.inObjs = L1C_objs  # type: List[L1C_object]
231
        self.reporting = reporting
Daniel Scheffler's avatar
Daniel Scheffler committed
232
233
        self.ac_input = {}  # set by self.run_atmospheric_correction()
        self.results = None  # direct output of external atmCorr module (set by run_atmospheric_correction)
234
        self.proc_info = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
235
        self.outObjs = []  # atmospherically corrected L1C objects
236
237

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

240
241
242
243
        if not re.search('Sentinel-2', self.inObjs[0].satellite, re.I):
            warnings.warn('Calculation of acquisition geometry arrays is currently only validated for Sentinel-2!')
            # validation possible by comparing S2 angles provided by ESA with own angles

244
245
246
247
248
    @property
    def logger(self):
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
249
            if len(self.inObjs) == 1:
250
251
252
253
254
255
256
257
258
259
                # 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
260
                    fileHandler = logging.FileHandler(path_logfile, mode='a')
261
262
263
264
265
                    fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
                    fileHandler.setLevel(logging.DEBUG)

                    logger_atmCorr.addHandler(fileHandler)

Daniel Scheffler's avatar
Daniel Scheffler committed
266
267
                    inObj.close_GMS_loggers()

268
269
270
271
272
273
            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
274
            "AtmCorr.logger can not be set to %s." % logger
275
276
277
278
279
280
281
282
        if logger in ['not set', None]:
            self._logger.close()
            self._logger = logger
        else:
            self._logger = logger

    @logger.deleter
    def logger(self):
283
284
285
        if self._logger not in [None, 'not set']:
            self._logger.close()
            self._logger = None
286

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

289
290
291
292
293
294
295
296
    @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
297
                              (' (%s)' % obj.subsystem if obj.subsystem else '') +
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
                              '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:
324
325
            data_dict = {}

326
            for inObj in self.inObjs:
327
                for bandN, bandIdx in inObj.arr.bandnames.items():
328
                    if bandN not in data_dict:
Daniel Scheffler's avatar
Daniel Scheffler committed
329
330
331
332
                        # float32! -> conversion to np.float16 will convert -9999 to -10000
                        arr2pass = inObj.arr[:, :, bandIdx].astype(np.float32)
                        arr2pass[arr2pass == inObj.arr.nodata] = np.nan  # set nodata values to np.nan
                        data_dict[bandN] = (arr2pass / inObj.meta_odict['ScaleFactor']).astype(np.float16)
333
                    else:
334
                        inObj.logger.warning("Band '%s' cannot be included into atmospheric correction because it "
Daniel Scheffler's avatar
Daniel Scheffler committed
335
                                             "exists multiple times." % bandN)
336

337
            # validate: data must have all bands needed for AC
Daniel Scheffler's avatar
Daniel Scheffler committed
338
339
            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]
340
341
            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
342
                                   'Received: %s' % (str(all_bNs_AC), str(list(sorted(data_dict.keys())))))
343
344
345

            self._data = data_dict

346
347
348
349
350
351
352
353
        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

354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    @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):
373
        """Returns S2A tile name.
374
        NOTE: this is only needed if no DEM is passed to ac_gms
375
376
377
378
379

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

Daniel Scheffler's avatar
Daniel Scheffler committed
380
        return ''  # FIXME
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408

    @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

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
    @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
        """
439
        # TODO add SRF object
Daniel Scheffler's avatar
Daniel Scheffler committed
440

441
        if not self._metadata:
Daniel Scheffler's avatar
Daniel Scheffler committed
442
            del self.logger  # otherwise each input object would have multiple fileHandlers
443

Daniel Scheffler's avatar
Daniel Scheffler committed
444
445
446
447
448
449
450
451
452
453
454
455
456
            metadata = dict(
                U=self.inObjs[0].meta_odict['EarthSunDist'],
                SENSING_TIME=self.inObjs[0].acq_datetime,
                # SENSING_TIME=datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z'),
                viewing_zenith=self._meta_get_viewing_zenith(),
                viewing_azimuth=self._meta_get_viewing_azimuth(),
                relative_viewing_azimuth=self._meta_get_relative_viewing_azimuth(),
                sun_mean_azimuth=self.inObjs[0].meta_odict['SunAzimuth'],
                sun_mean_zenith=90 - self.inObjs[0].meta_odict['SunElevation'],
                solar_irradiance=self._meta_get_solar_irradiance(),
                aux_data=self._meta_get_aux_data(),
                spatial_samplings=self._meta_get_spatial_samplings()
            )
457
458

            self._metadata = metadata
459
460
461

        return self._metadata

462
463
464
465
    @property
    def options(self):
        """Returns a dictionary containing AC options.
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
466
        # type: -> dict
467
468
469
470
        if self._options:
            return self._options
        else:
            self._options = self.inObjs[0].ac_options
Daniel Scheffler's avatar
Daniel Scheffler committed
471
            self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
472
            self._options["report"]["reporting"] = self.reporting
473
474
            return self._options

475
    def _meta_get_spatial_samplings(self):
476
477
478
        """

        :return:
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
         {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}}
497
        """
498
499
        # set corner coordinates and dims
        spatial_samplings = {}
500
501
502

        for inObj in self.inObjs:

503
504
505
506
507
            # 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.')
508

509
510
            # set spatial information
            spatial_samplings[inObj.arr.xgsd] = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
511
512
513
514
515
516
                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)
517

518
519
520
        return spatial_samplings

    def _meta_get_solar_irradiance(self):
521
522
523
        """

        :return:
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
        {'B01': 1913.57,
         'B02': 1941.63,
         'B03': 1822.61,
         'B04': 1512.79,
         'B05': 1425.56,
         'B06': 1288.32,
         'B07': 1163.19,
         'B08': 1036.39,
         'B09': 813.04,
         'B10': 367.15,
         'B11': 245.59,
         'B12': 85.25,
         'B8A': 955.19}
        """

        solar_irradiance = {}

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

    def _meta_get_viewing_zenith(self):
        """

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

        viewing_zenith = {}

Daniel Scheffler's avatar
Daniel Scheffler committed
555
        for inObj in self.inObjs:  # type: L1C_object
556
            for bandN, bandIdx in inObj.arr.bandnames.items():
557
                if bandN not in viewing_zenith:
558
559
                    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
560
                    # viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
561
562
563
564
565
566
567
568
569
570
        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
571
        for inObj in self.inObjs:  # type: L1C_object
572
            for bandN, bandIdx in inObj.arr.bandnames.items():
573
                if bandN not in viewing_azimuth:
Daniel Scheffler's avatar
Daniel Scheffler committed
574
                    arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim == 3 else inObj.VAA_arr
575
                    viewing_azimuth[bandN] = arr2pass.astype(np.float16)
Daniel Scheffler's avatar
Daniel Scheffler committed
576
                    # viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
577

578
579
580
581
582
583
        return viewing_azimuth

    def _meta_get_relative_viewing_azimuth(self):
        """

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

586
587
        relative_viewing_azimuth = {}

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 relative_viewing_azimuth:
591
592
                    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
593
594
                    # relative_viewing_azimuth[bandN] = \
                    #     inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr
595

596
        return relative_viewing_azimuth
597

598
599
600
601
602
603
604
605
    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
606
607
            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°
608
            # FIXME correct to reduce resolution here by factor 10?
609
610
611
612
613
614
615
616
617
618
619
620
        )

        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
621
            inObj4dem = [obj for obj in self.inObjs if obj.arr.xgsd == 20]
622
623
624
            if not inObj4dem:
                self.logger.warning('Sentinel-2 20m subsystem could not be found. DEM passed to '
                                    'atmospheric correction might have wrong resolution.')
625
626
627
628
            inObj4dem = inObj4dem[0]
        else:
            inObj4dem = self.inObjs[0]

629
630
631
632
        try:
            dem = inObj4dem.dem[:].astype(np.float32)
        except Exception as e:
            dem = None
Daniel Scheffler's avatar
Daniel Scheffler committed
633
            self.logger.warning('A static elevation is assumed during atmospheric correction due to an error during '
634
635
636
                                '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:")
637
            self.logger.warning(traceback.format_exc())
638
639

        return dem
640
641

    def _get_srf(self):
642
        """Returns an instance of SRF in the same structure like sicor.sensors.SRF.SensorSRF
643
        """
644
645
646
        # FIXME calculation of center wavelengths within SRF() used not the GMS algorithm
        # SRF instance must be created for all bands and the previous proc level
        GMS_identifier_fullScene = self.inObjs[0].GMS_identifier
Daniel Scheffler's avatar
Daniel Scheffler committed
647
        GMS_identifier_fullScene['Subsystem'] = ''
648
649
650
        GMS_identifier_fullScene['proc_level'] = proc_chain[proc_chain.index(self.inObjs[0].proc_level) - 1]

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

652
653
654
655
656
657
    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:
        """

658
659
        tgt_res = self.inObjs[0].ac_options['cld_mask']['target_resolution']

660
661
        # check if input GMS objects provide a cloud mask
        avail_cloud_masks = {inObj.GMS_identifier['Subsystem']: inObj.mask_clouds for inObj in self.inObjs}
662
        no_avail_CMs = list(set(avail_cloud_masks.values())) == [None]
663
664

        # compute cloud mask if not already provided
665
666
        if no_avail_CMs:
            algorithm = CFG.job.cloud_masking_algorithm[self.inObjs[0].satellite]
667

668
669
            if algorithm == 'SICOR':
                return None
670

671
672
673
674
675
676
677
678
679
680
681
            else:
                # FMASK or Classical Bayesian
                try:
                    from .cloud_masking import Cloud_Mask_Creator

                    CMC = Cloud_Mask_Creator(self.inObjs[0], algorithm=algorithm)
                    CMC.calc_cloud_mask()
                    cm_geoarray = CMC.cloud_mask_geoarray
                    cm_array = CMC.cloud_mask_array
                    cm_legend = CMC.cloud_mask_legend
                except Exception as err:
682
683
                    self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
                    self.logger.error(traceback.format_exc())
684
                    return None
685

686
687
        else:
            # check if there is a cloud mask with suitable GSD
Daniel Scheffler's avatar
Daniel Scheffler committed
688
            inObjs2use = [obj for obj in self.inObjs if obj.mask_clouds is not None and obj.mask_clouds.xgsd == tgt_res]
689
690
            if not inObjs2use:
                raise ValueError('Error appending cloud mask to input arguments of atmospheric correction. No input '
Daniel Scheffler's avatar
Daniel Scheffler committed
691
                                 'GMS object provides a cloud mask with spatial resolution of %s.' % tgt_res)
692
693
694
695
696
697
698
699
            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)
700
            #    {'Clear': 10, 'Thick Clouds': 20, 'Thin Clouds': 30, 'Snow': 40}  # FIXME hardcoded
701
702
703
704
705
706

            # 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.')
707
708
709
710

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

711
712
        # get nodata value
        self.options['cld_mask']['nodata_value_mask'] = cm_geoarray.nodata
713

714
        # append cloud mask to input object with the same spatial resolution if there was no mask before
715
        for inObj in self.inObjs:
716
            if inObj.arr.xgsd == cm_geoarray.xgsd:
717
718
                inObj.mask_clouds = cm_geoarray
                inObj.build_combined_masks_array()
Daniel Scheffler's avatar
Daniel Scheffler committed
719
                break  # appending it to one inObj is enough
720

721
722
723
        return S2Mask(mask_array=cm_array,
                      mask_legend=cm_legend,
                      geo_coding=cm_geocoding)
724

725
726
    def run_atmospheric_correction(self, dump_ac_input=False):
        # type: (bool) -> list
727
728
729
        """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
        containing surface reflectance.

730
731
        :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
732
        """
733
734

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

737
        rs_data = dict(
Daniel Scheffler's avatar
Daniel Scheffler committed
738
739
740
741
742
743
744
745
746
            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
747
        )  # NOTE: all keys of this dict are later converted to attributes of RSImage
748

Daniel Scheffler's avatar
Daniel Scheffler committed
749
        script = False
750

751
752
753
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

754
        self.ac_input = dict(
755
            rs_image=rs_image,
Daniel Scheffler's avatar
Daniel Scheffler committed
756
            options=self.options,  # type: dict
757
758
            logger=repr(self.logger),  # only a string
            script=script
759
        )
760

761
762
763
764
        # path_dump = self.inObjs[0].pathGen.get_path_ac_input_dump()
        # with open(path_dump, 'wb') as outF:
        #     dill.dump(self.ac_input, outF)

765
        # run AC
766
        self.logger.info('Atmospheric correction started.')
767
        try:
768
            rs_image.logger = self.logger
769
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
770

771
        except Exception as e:
772
            # serialialize AC input
773
774
775
776
777
778
            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
779
                                  % path_dump)
780
781

            # delete AC input arrays
Daniel Scheffler's avatar
Daniel Scheffler committed
782
            for inObj in self.inObjs:  # type: L1C_object
783
784
                inObj.delete_ac_input_arrays()

785
786
            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'
787
                              % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
788
            self.logger.error(traceback.format_exc())
789
            # TODO include that in the job summary
790

791
792
            return list(self.inObjs)

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

796
        # join results
Daniel Scheffler's avatar
Daniel Scheffler committed
797
        self._join_results_to_inObjs()  # sets self.outObjs
798

799
800
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
801

802
803
804
        return self.outObjs

    def _join_results_to_inObjs(self):
805
806
807
        """
        Join results of atmospheric correction to the input GMS objects.
        """
808

809
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
Daniel Scheffler's avatar
Daniel Scheffler committed
810
811
812
        # delete logger
        # -> otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)
        del self.logger
813
814
815
816
817
818
819
820

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

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

822
823
824
825
        self.outObjs = self.inObjs

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

830
        if self.results.data_ac is not None:
831
            for inObj in self.inObjs:
832
                assert isinstance(inObj, L1B_object)
833
                nodata = self.results.nodata[inObj.arr.xgsd]  # 2D mask with True outside of image coverage
Daniel Scheffler's avatar
Daniel Scheffler committed
834
                ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
835
                out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
836

837
838
839
840
841
842
                # update metadata
                inObj.arr_desc = 'BOA_Ref'
                inObj.MetaObj.bands = len(self.results.data_ac)
                inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.usecase.scale_factor_BOARef
                inObj.MetaObj.LayerBandsAssignment = out_LBA
                inObj.MetaObj.filter_layerdependent_metadata()
843
                inObj.meta_odict = inObj.MetaObj.to_odict()  # actually auto-updated by getter
844

845
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config
846
847
                # FIXME AC output nodata values = 0 -> new nodata areas but mask not updated
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
848
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
849
850
851
                surf_refl *= CFG.usecase.scale_factor_BOARef  # scale using scale factor (output is float16)
                # FIXME really set AC nodata values to GMS outZero?
                surf_refl[nodata] = oZ_refl  # overwrite AC nodata values with GMS outZero
Daniel Scheffler's avatar
Daniel Scheffler committed
852
853
                # apply the original nodata mask (indicating background values)
                surf_refl[np.array(inObj.mask_nodata) == False] = oF_refl
854

Daniel Scheffler's avatar
Daniel Scheffler committed
855
                if self.results.bad_data_value is np.nan:
856
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
857
                else:
Daniel Scheffler's avatar
Daniel Scheffler committed
858
859
                    surf_refl[
                        surf_refl == self.results.bad_data_value] = oF_refl  # FIXME meaningful to set AC nans to -9999?
860
861
862

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

865
866
867
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
868

869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
    def _join_data_errors(self):
        """
        Join ERRORS ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config.
        """

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

                ac_errors = np.dstack((self.results.data_errors[bandN] for bandN in ac_bandNs))
                ac_errors *= CFG.usecase.scale_factor_errors_ac  # scale using scale factor (output is float16)
                out_dtype = np.int8 if CFG.usecase.scale_factor_errors_ac <= 255 else np.int16
                ac_errors[nodata] = get_outFillZeroSaturated(out_dtype)[0]
                ac_errors = ac_errors.astype(out_dtype)
                inObj.ac_errors = ac_errors  # setter generates a GeoArray with the same bandnames like inObj.arr
                # TODO how to handle nans?
        else:
            self.logger.warning("Atmospheric correction did not provide a 'data_errors' array. Maybe due to "
                                "missing SNR model? GMS_object.ac_errors kept None.")

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

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

899
900
            joined = False
            for inObj in self.inObjs:
901
902
                # delete all previous cloud masks
                del inObj.mask_clouds
903
904
905
906

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

                    # set cloud mask nodata value
911
                    tgt_nodata = get_outFillZeroSaturated(mask_clouds_ac.dtype)[0]
912
913
                    ac_out_nodata = self.ac_input['options']['cld_mask']['nodata_value_mask']
                    if tgt_nodata not in self.results.mask_clouds.mask_legend.keys():
914
                        inObj.mask_clouds[inObj.mask_clouds[:] == ac_out_nodata] = tgt_nodata
915
916
917
918
                        mask_clouds_nodata = tgt_nodata
                    else:
                        warnings.warn('The cloud mask from AC output already uses the desired nodata value %s for the '
                                      'class %s. Using AC output nodata value %s.'
919
                                      % (tgt_nodata, self.results.mask_clouds.mask_legend[tgt_nodata], ac_out_nodata))
920
921
922
923
                        mask_clouds_nodata = ac_out_nodata

                    inObj.mask_clouds.nodata = mask_clouds_nodata

924
                    joined = True
925

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

930
        else:
931
932
            self.logger.warning("Atmospheric correction did not provide a 'mask_clouds.mask_array' array. "
                                "GMS_object.mask_clouds kept None.")
933

934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
    def _join_mask_confidence_array(self):
        """
        Join confidence array for mask_clouds.
        """

        if self.results.mask_clouds.mask_confidence_array is not None:
            cfd_arr = self.results.mask_clouds.mask_confidence_array  # float32 2D array, scaled [0-1, nodata 255]
            cfd_arr[cfd_arr == self.ac_input['options']['cld_mask']['nodata_value_mask']] = -1
            cfd_arr = (cfd_arr * CFG.usecase.scale_factor_BOARef).astype(np.int16)
            cfd_arr[cfd_arr == -CFG.usecase.scale_factor_BOARef] = get_outFillZeroSaturated(cfd_arr.dtype)[0]

            joined = False
            for inObj in self.inObjs:

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

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

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