L1C_P.py 41.2 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
# -*- coding: utf-8 -*-
###############################################################################
#
#   Level 1C Processor:
5
#
Daniel Scheffler's avatar
Daniel Scheffler committed
6
#   Performed operations:
7
#   Atmospheric correction of TOA-reflectance data:
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
12
13
#
#   Written by Daniel Scheffler
#
###############################################################################

########################### Library import ####################################
14
import warnings
15
16
import re
import logging
17
import dill
18
import traceback
19

20
import numpy as np
21
22
23
24
try:
    from osgeo import osr
except ImportError:
    import osr
Daniel Scheffler's avatar
Daniel Scheffler committed
25

26
from geoarray import GeoArray
27
from py_tools_ds.geo.map_info import mapinfo2geotransform
28
29

from ..config import GMS_config       as CFG
30
31
from .        import GEOPROCESSING    as GEOP
from .L1B_P   import L1B_object
32
from ..model.METADATA import get_LayerBandsAssignment
33
from ..misc.definition_dicts import get_outFillZeroSaturated, proc_chain, get_mask_classdefinition
34
from ..io.Input_reader import SRF
35
# from .cloud_masking import Cloud_Mask_Creator  # circular dependencies
36

37
from sicor.sicor_ac import ac_gms
38
from sicor.sensors import RSImage
39
from sicor.Mask import S2Mask
40

Daniel Scheffler's avatar
Daniel Scheffler committed
41
42

########################### core functions ####################################
43
class L1C_object(L1B_object):
44
45
46
47
48
49
50
    def __init__(self, L1B_obj=None):
        super().__init__()

        if L1B_obj:
            # populate attributes
            [setattr(self, key, value) for key,value in L1B_obj.__dict__.items()]

51
52
53
54
55
56
57
58
        # 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

59
        self.proc_level = 'L1C'
60

61

62
63
64
    @property
    def lonlat_arr(self):
        """Calculates pixelwise 2D-array with longitude and latitude coordinates.
65

66
67
68
69
70
71
72
73
        :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'],
74
                                            meshwidth=10, # for faster processing
75
76
77
                                            nodata_mask=None, # dont overwrite areas outside the image with nodata
                                            outFill=get_outFillZeroSaturated(np.float32)[0])[0]
        return self._lonlat_arr
78
79


80
81
82
    @lonlat_arr.setter
    def lonlat_arr(self, lonlat_arr):
        self._lonlat_arr = lonlat_arr
83

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

    @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,
                                                                          meshwidth=10, # for faster processing
                                                                          subset=None,
                                                                          bandwise=0)
            else:
                self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
102
103
104
105
106
107
                                                    float(self.meta_odict['ViewingAngle']),
                                                    float(self.meta_odict['FieldOfView']),
                                                    self.logger,
                                                    nodata_mask=None,  # dont overwrite areas outside the image with nodata
                                                    outFill=get_outFillZeroSaturated(np.float32)[0],
                                                    meshwidth=10) # for faster processing
108
109
110
111
112
113
        return self._VZA_arr


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

115

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

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

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


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

146

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
    @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:
            _ = self.SZA_arr # getter also sets self._SAA_arr
        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
209

210

211
212
213
214
215
216
217
    def delete_ac_input_arrays(self):
        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
        del self.dem # uses deleter # would have to be resampled when writing MGRS tiles -> better to directly warp it to the output dims and projection
218

219
220
221


class AtmCorr(object):
222
    def __init__(self, *L1C_objs, reporting=False):
223
        """Wrapper around atmospheric correction by Andre Hollstein, GFZ Potsdam
224
225
226
227
228
229

        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
        """
230
        # FIXME not yet usable for data < 2012 due to missing ECMWF archive
231
232
233
        L1C_objs = L1C_objs if isinstance(L1C_objs, tuple) else (L1C_objs,)

        # hidden attributes
234
        self._logger   = None
235
236
237
238
        self._GSDs     = []
        self._data     = {}
        self._metadata = {}
        self._nodata   = {}
239
        self._band_spatial_sampling = {}
240
        self._options  = {}
241
242
243
244
245
246

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

247
        self.inObjs    = L1C_objs
248
        self.reporting = reporting
249
250
251
252
253
254
255
256
        self.ac_input  = {} # set by self.run_atmospheric_correction()
        self.results   = None # direct output of external atmCorr module (set by run_atmospheric_correction)
        self.proc_info = {}
        self.outObjs   = [] # atmospherically corrected L1C objects

        # append AtmCorr object to input L1C objects
        #[setattr(L1C_obj, 'AtmCorr', self) for L1C_obj in self.inObjs] # too big for serialization

257
258
259
260
        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

261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283

    @property
    def logger(self):
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
            if len(self.inObjs)==1:
                # 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()
                    fileHandler  = logging.FileHandler(path_logfile, mode='a')
                    fileHandler.setFormatter(logger_atmCorr.formatter_fileH)
                    fileHandler.setLevel(logging.DEBUG)

                    logger_atmCorr.addHandler(fileHandler)

Daniel Scheffler's avatar
Daniel Scheffler committed
284
285
                    inObj.close_GMS_loggers()

286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
            self._logger = logger_atmCorr
            return self._logger


    @logger.setter
    def logger(self, logger):
        assert isinstance(logger, logging.Logger) or logger in ['not set', None], \
            "AtmCorr.logger can not be set to %s." %logger
        if logger in ['not set', None]:
            self._logger.close()
            self._logger = logger
        else:
            self._logger = logger


    @logger.deleter
    def logger(self):
303
304
305
        if self._logger not in [None, 'not set']:
            self._logger.close()
            self._logger = None
306

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

309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346

    @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 +
                              (' (%s)'%obj.subsystem if obj.subsystem else '') +
                              '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:
347
348
            data_dict = {}

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

359
360
361
362
363
364
365
366
367
            # validate: data must have all bands needed for AC
            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]
            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'
                                   'Received: %s' %(str(all_bNs_AC), str(list(sorted(data_dict.keys())))))

            self._data = data_dict

368
369
370
371
372
373
374
375
376
377
        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


378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
    @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):
398
        """Returns S2A tile name.
399
        NOTE: this is only needed if no DEM is passed to ac_gms
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435

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

        return '' # FIXME


    @property
    def band_spatial_sampling(self):
        """

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

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


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

        :return:
            ___ attribute: metadata, type:<class 'dict'>
            ______ key:spatial_samplings
            _________ key:60.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 1830
            ____________ key:XDIM, value_type:<class 'int'>, repr: 60
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 1830
            ____________ key:YDIM, value_type:<class 'int'>, repr: -60
            _________ key:10.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 10980
            ____________ key:XDIM, value_type:<class 'int'>, repr: 10
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 10980
            ____________ key:YDIM, value_type:<class 'int'>, repr: -10
            _________ key:20.0
            ____________ key:ULY, value_type:<class 'int'>, repr: 4900020
            ____________ key:NCOLS, value_type:<class 'int'>, repr: 5490
            ____________ key:XDIM, value_type:<class 'int'>, repr: 20
            ____________ key:ULX, value_type:<class 'int'>, repr: 600000
            ____________ key:NROWS, value_type:<class 'int'>, repr: 5490
            ____________ key:YDIM, value_type:<class 'int'>, repr: -20
            ______ key:SENSING_TIME, value_type:<class 'datetime.datetime'>, repr: 2016-03-26 10:34:06.538000+00:00
        """
466
467
        # TODO add SRF object
        metadata = {}
468
        if not self._metadata:
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
            del self.logger # otherwise each input object would have multiple fileHandlers

            metadata['U']                         = self.inObjs[0].meta_odict['EarthSunDist']
            metadata['SENSING_TIME']              = self.inObjs[0].acq_datetime
            #metadata['SENSING_TIME']              = datetime.strptime('2015-08-12 10:40:21 +0000', '%Y-%m-%d %H:%M:%S %z')
            metadata['viewing_zenith']            = self._meta_get_viewing_zenith()
            metadata['viewing_azimuth']           = self._meta_get_viewing_azimuth()
            metadata['relative_viewing_azimuth']  = self._meta_get_relative_viewing_azimuth()
            metadata['sun_mean_azimuth']          = self.inObjs[0].meta_odict['SunAzimuth']
            metadata['sun_mean_zenith']           = 90-self.inObjs[0].meta_odict['SunElevation']
            metadata['solar_irradiance']          = self._meta_get_solar_irradiance()
            metadata['aux_data']                  = self._meta_get_aux_data()
            metadata['spatial_samplings']         = self._meta_get_spatial_samplings()

            self._metadata = metadata
484
485
486

        return self._metadata

487
488
489
490
491
492
493
494
    @property
    def options(self):
        """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
495
            self._options["AC"]['bands'] = [b for b in self.data.keys() if b in self._options["AC"]['bands']]
496
            self._options["report"]["reporting"] = self.reporting
497
498
            return self._options

499

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

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

        for inObj in self.inObjs:

528
529
530
531
532
            # 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.')
533

534
535
            # set spatial information
            spatial_samplings[inObj.arr.xgsd] = dict(
536
537
538
539
540
541
                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)
542

543
544
545
546
        return spatial_samplings


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

        :return:
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
        {'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 = {}

        for inObj in self.inObjs:
583
            for bandN, bandIdx in inObj.arr.bandnames.items():
584
                if bandN not in viewing_zenith:
585
586
587
                    arr2pass = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim == 3 else inObj.VZA_arr
                    viewing_zenith[bandN] = arr2pass.astype(np.float16)
                    #viewing_zenith[bandN] = inObj.VZA_arr[:, :, bandIdx] if inObj.VZA_arr.ndim==3 else inObj.VZA_arr
588
589
590
591
592
593
594
595
596
597
598
599
        return viewing_zenith


    def _meta_get_viewing_azimuth(self):
        """

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

        viewing_azimuth = {}

        for inObj in self.inObjs:
600
            for bandN, bandIdx in inObj.arr.bandnames.items():
601
                if bandN not in viewing_azimuth:
602
603
604
605
                    arr2pass = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr
                    viewing_azimuth[bandN] = arr2pass.astype(np.float16)
                    #viewing_azimuth[bandN] = inObj.VAA_arr[:, :, bandIdx] if inObj.VAA_arr.ndim==3 else inObj.VAA_arr

606
607
608
609
610
611
612
        return viewing_azimuth


    def _meta_get_relative_viewing_azimuth(self):
        """

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

615
616
        relative_viewing_azimuth = {}

617
        for inObj in self.inObjs:
618
            for bandN, bandIdx in inObj.arr.bandnames.items():
619
                if bandN not in relative_viewing_azimuth:
620
621
622
623
624
                    arr2pass = inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim == 3 else inObj.RAA_arr
                    relative_viewing_azimuth[bandN] = arr2pass.astype(np.float16)
                    #relative_viewing_azimuth[bandN] = \
                    #    inObj.RAA_arr[:, :, bandIdx] if inObj.RAA_arr.ndim==3 else inObj.RAA_arr

625
        return relative_viewing_azimuth
626
627


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

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

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

        return dem
671
672
673


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


685
686
687
688
689
690
    def _get_mask_clouds(self):
        """Returns an instance of S2Mask in case cloud mask is given by input GMS objects. Otherwise None is returned.

        :return:
        """

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

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

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

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

704
705
706
707
708
709
710
711
712
713
714
            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:
715
716
                    self.logger.error('\nAn error occurred during FMASK cloud masking. Error message was: ')
                    self.logger.error(traceback.format_exc())
717
                    return None
718

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

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

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

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

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

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


759
760
    def run_atmospheric_correction(self, dump_ac_input=False):
        # type: (bool) -> list
761
762
763
        """Collects all input data for atmospheric correction, runs the AC and returns the corrected L1C objects
        containing surface reflectance.

764
765
        :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
766
        """
767
768

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

771
772
773
774
775
776
        rs_data = dict(
            data                  = self.data,
            metadata              = self.metadata,
            nodata                = self.nodata,
            band_spatial_sampling = self.band_spatial_sampling,
            tile_name             = self.tile_name,
777
            dem                   = self._get_dem(),
778
779
            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
780
        ) # NOTE: all keys of this dict are later converted to attributes of RSImage
781

782
        script  = False
783

784
785
786
        # create an instance of RSImage
        rs_image = RSImage(**rs_data)

787
        self.ac_input = dict(
788
789
            rs_image = rs_image,
            options  = self.options,
790
            logger   = repr(self.logger), # only a string
791
792
            script   = script
        )
793

794
        # run AC
795
        self.logger.info('Atmospheric correction started.')
796
        try:
797
            rs_image.logger = self.logger
798
            self.results = ac_gms(rs_image, self.options, logger=self.logger, script=script)
799

800
        except Exception as e:
801
            # serialialize AC input
802
803
804
805
806
807
808
            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.'
                                  %path_dump)
809
810
811
812
813

            # delete AC input arrays
            for inObj in self.inObjs:
                inObj.delete_ac_input_arrays()

814
815
            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'
816
                              % (self.inObjs[0].scene_ID, self.inObjs[0].entity_ID, repr(e)))
817
            self.logger.error(traceback.format_exc())
818
            # TODO include that in the job summary
819

820
821
            return list(self.inObjs)

822
823

        # get processing infos
824
        self.proc_info = self.ac_input['options']['processing'] # FIXME this is not appended to GMS objects
825

826
827
        # join results
        self._join_results_to_inObjs() # sets self.outObjs
828

829
830
        # delete input arrays that are not needed anymore
        [inObj.delete_ac_input_arrays() for inObj in self.inObjs]
831

832
833
834
835
        return self.outObjs


    def _join_results_to_inObjs(self):
836
837
838
        """
        Join results of atmospheric correction to the input GMS objects.
        """
839

840
841
842
843
844
845
846
847
848
849
        self.logger.info('Joining results of atmospheric correction to input GMS objects.')
        del self.logger  # otherwise logging in inObjs would open a second FileHandler to the same file (which is permitted)

        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]
850

851
852
853
854
855
856
857
        self.outObjs = self.inObjs


    def _join_data_ac(self):
        """
        Join ATMOSPHERICALLY CORRECTED ARRAY as 3D int8 or int16 BOA reflectance array, scaled to scale factor from config.
        """
858

859
        if self.results.data_ac is not None:
860
            for inObj in self.inObjs:
861
                assert isinstance(inObj, L1B_object)
862
                nodata = self.results.nodata[inObj.arr.xgsd]  # 2D mask with True outside of image coverage
Daniel Scheffler's avatar
Daniel Scheffler committed
863
                ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
864
                out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
865

866
867
868
869
870
871
                # 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()
872
                inObj.meta_odict = inObj.MetaObj.to_odict()  # actually auto-updated by getter
873

874
                # join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config
875
876
                # FIXME AC output nodata values = 0 -> new nodata areas but mask not updated
                oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
877
                surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
878
879
880
881
882
                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
                surf_refl[np.array(inObj.mask_nodata) == False] = oF_refl  # apply the original nodata mask (indicating background values)

Daniel Scheffler's avatar
Daniel Scheffler committed
883
                if self.results.bad_data_value is np.nan:
884
                    surf_refl[np.isnan(surf_refl)] = oF_refl
Daniel Scheffler's avatar
Daniel Scheffler committed
885
                else:
886
                    surf_refl[surf_refl == self.results.bad_data_value] = oF_refl # FIXME meaningful to set AC nans to -9999?
887
888
889

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

892
893
894
        else:
            self.logger.warning('Atmospheric correction did not return a result for the input array. '
                                'Thus the output keeps NOT atmospherically corrected.')
895

896

897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
    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
927

928
929
            joined = False
            for inObj in self.inObjs:
930
931
                # delete all previous cloud masks
                del inObj.mask_clouds
932
933
934
935

                # 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
936
937
                    inObj.mask_clouds.legend = self.results.mask_clouds.mask_legend  # dict(value=string, string=value))
                    # FIXME legend is not used later
938
939

                    # set cloud mask nodata value
940
                    tgt_nodata = get_outFillZeroSaturated(mask_clouds_ac.dtype)[0]
941
942
                    ac_out_nodata = self.ac_input['options']['cld_mask']['nodata_value_mask']
                    if tgt_nodata not in self.results.mask_clouds.mask_legend.keys():
943
                        inObj.mask_clouds[inObj.mask_clouds[:] == ac_out_nodata] = tgt_nodata
944
945
946
947
                        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.'
948
                                      % (tgt_nodata, self.results.mask_clouds.mask_legend[tgt_nodata], ac_out_nodata))
949
950
951
952
                        mask_clouds_nodata = ac_out_nodata

                    inObj.mask_clouds.nodata = mask_clouds_nodata

953
                    joined = True
954

955
956
957
958
            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.')

959
        else:
960
961
            self.logger.warning("Atmospheric correction did not provide a 'mask_clouds.mask_array' array. "
                                "GMS_object.mask_clouds kept None.")
962
963


964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
    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
984

985
986
987
            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.')
988

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