gms_object.py 115 KB
Newer Older
1
2
3
4
5
# -*- coding: utf-8 -*-

import collections
import copy
import datetime
6
import functools
7
8
9
10
11
12
13
import glob
import json
import os
import re
import shutil
import sys
import warnings
14
import logging
15
from collections import OrderedDict
16
from itertools import chain
17
from typing import Iterable, List, Union, TYPE_CHECKING  # noqa F401  # flake8 issue
18
19
20

import numpy as np
import spectral
21
from spectral.io import envi
22
from numba import jit
23
from pandas import DataFrame, read_csv
24
from nested_dict import nested_dict
25

26
27
28
29
try:
    from osgeo import gdalnumeric
except ImportError:
    import gdalnumeric
30

31
from geoarray import GeoArray, NoDataMask, CloudMask
32
from py_tools_ds.geo.coord_grid import is_coord_grid_equal
33
from py_tools_ds.geo.projection import EPSG2WKT, WKT2EPSG, get_UTMzone, isProjectedOrGeographic
34
35
from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
36
37
from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX, imXY2mapXY, imYX2mapYX
from py_tools_ds.numeric.array import get_array_tilebounds
38
from sicor.options import get_options as get_ac_options
39

40
41
from ..misc.logging import GMS_logger as DatasetLogger
from ..misc import database_tools as DB_T
42
43
44
45
from ..misc import path_generator as PG
from ..model.mgrs_tile import MGRS_tile
from ..model.metadata import \
    METADATA, get_dict_LayerOptTherm, metaDict_to_metaODict, get_LayerBandsAssignment, layerdependent_metadata
46
from ..options.config import GMS_config as CFG
47
48
49
from ..algorithms import geoprocessing as GEOP
from ..io import input_reader as INP_R
from ..io import output_writer as OUT_W
50
51
from ..misc import helper_functions as HLP_F
from ..misc import definition_dicts as DEF_D
52
from ..misc.locks import IOLock
53

54
55
56
if TYPE_CHECKING:
    from ..algorithms.L1C_P import L1C_object  # noqa F401  # flake8 issue

57
__author__ = 'Daniel Scheffler'
58
59


60
class GMS_object(object):
61
62
63
64
    # class attributes
    # NOTE: these attributes can be modified and seen by ALL GMS_object instances
    proc_status_all_GMSobjs = nested_dict()

65
    def __init__(self, pathImage=''):
66
        # add private attributes
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
        self._logger = None
        self._loggers_disabled = False
        self._log = ''
        self._MetaObj = None
        self._LayerBandsAssignment = None
        self._coreg_needed = None
        self._resamp_needed = None
        self._arr = None
        self._mask_nodata = None
        self._mask_clouds = None
        self._mask_clouds_confidence = None
        self._masks = None
        self._dem = None
        self._pathGen = None
        self._ac_options = {}
        self._ac_errors = None
        self._spat_homo_errors = None
        self._spec_homo_errors = None
        self._accuracy_layers = None
86
        self._dict_LayerOptTherm = None
87
        self._cloud_masking_algorithm = None
88
        self._coreg_info = None
89

90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
        # defaults
        self.proc_level = 'init'
        self.image_type = ''
        self.satellite = ''
        self.sensor = ''
        self.subsystem = ''
        self.sensormode = ''
        self.acq_datetime = None  # also includes time, set by from_disk() and within L1A_P
        self.entity_ID = ''
        self.scene_ID = -9999
        self.filename = ''
        self.dataset_ID = -9999
        self.outInterleave = 'bsq'
        self.VAA_mean = None  # set by self.calc_mean_VAA()
        self.corner_lonlat = None
        self.trueDataCornerPos = []  # set by self.calc_corner_positions()
        self.trueDataCornerLonLat = []  # set by self.calc_corner_positions()
        self.fullSceneCornerPos = []  # set by self.calc_corner_positions()
        self.fullSceneCornerLonLat = []  # set by self.calc_corner_positions()
        # rows,cols,bands of the full scene (not of the subset as possibly represented by self.arr.shape)
        self.shape_fullArr = [None, None, None]
        self.arr_shape = 'cube'
        self.arr_desc = ''  # description of data units for self.arr
        self.arr_pos = None  # <tuple> in the form ((row_start,row_end),(col_start,col_end))
        self.tile_pos = None  # <list>, filled by self.get_tilepos()
        self.GeoTransProj_ok = True  # set by self.validate_GeoTransProj_GeoAlign()
        self.GeoAlign_ok = True  # set by self.validate_GeoTransProj_GeoAlign()
        self.masks_meta = {}  # set by self.build_L1A_masks()
        # self.CLD_obj               = CLD_P.GmsCloudClassifier(classifier=self.path_cloud_class_obj)

        # set pathes
        self.path_archive = ''
        self.path_procdata = ''
        self.ExtractedFolder = ''
        self.baseN = ''
        self.path_logfile = ''
        self.path_archive_valid = False
        self.path_InFilePreprocessor = None
        self.path_MetaPreprocessor = None

130
        self.job_ID = CFG.ID
131
132
133
        self.dataset_ID = -9999 if self.scene_ID == -9999 else \
            int(DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['datasetid'],
                                                {'id': self.scene_ID})[0][0])
134
135
        self.scenes_proc_ID = None  # set by Output writer after creation/update of db record in table scenes_proc
        self.mgrs_tiles_proc_ID = None  # set by Output writer after creation/update of db rec in table mgrs_tiles_proc
136
        self.MGRS_info = None
Daniel Scheffler's avatar
Daniel Scheffler committed
137
138
        self.lonlat_arr = None  # set by self.write_tiles_to_ENVIfile
        self.trueDataCornerUTM = None  # set by self.from_tiles
139
140

        # set pathes
141
142
143
144
        self.path_cloud_class_obj = ''

        # handle initialization arguments
        if pathImage:
145
            self.arr = pathImage  # run the setter for 'arr' which creates an Instance of GeoArray
146
147
148
149

    def __getstate__(self):
        """Defines how the attributes of GMS object are pickled."""

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
150
        self.close_loggers()
151
        del self.pathGen  # path generator can only be used for the current processing level
152
153

        # delete arrays if their in-mem size is to big to be pickled
154
        # => (avoids MaybeEncodingError: Error sending result: '[<gms_preprocessing.algorithms.L2C_P.L2C_object
155
        #    object at 0x7fc44f6399e8>]'. Reason: 'error("'i' format requires -2147483648 <= number <= 2147483647",)')
156
        if self.proc_level == 'L2C' and CFG.inmem_serialization:
157
            self.flush_array_data()
158

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
159
160
        return self.__dict__

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
    def __setstate__(self, ObjDict):
        """Defines how the attributes of GMS object are unpickled."""

        self.__dict__ = ObjDict
        # TODO unpickle meta to MetaObj

    def __deepcopy__(self, memodict={}):
        """Returns a deepcopy of the object excluding loggers because loggers are not serializable."""

        cls = self.__class__
        result = cls.__new__(cls)
        self.close_loggers()
        del self.pathGen  # has a logger
        memodict[id(self)] = result

        for k, v in self.__dict__.items():
            setattr(result, k, copy.deepcopy(v, memodict))
        return result

180
    def set_pathes(self):
181
182
183
184
185
186
        self.baseN = self.pathGen.get_baseN()
        self.path_procdata = self.pathGen.get_path_procdata()
        self.ExtractedFolder = self.pathGen.get_path_tempdir()
        self.path_logfile = self.pathGen.get_path_logfile()
        self.pathGen = PG.path_generator(self.__dict__)  # passes a logger in addition to previous attributes
        self.path_archive = self.pathGen.get_local_archive_path_baseN()
187

188
        if not CFG.inmem_serialization:
Daniel Scheffler's avatar
Daniel Scheffler committed
189
            self.path_InFilePreprocessor = os.path.join(self.ExtractedFolder, '%s%s_DN.bsq'
190
191
                                                        % (self.entity_ID,
                                                           ('_%s' % self.subsystem if self.subsystem else '')))
192
        else:  # keep data in memory
193
            self.path_InFilePreprocessor = None  # None: keeps all produced data in memory (numpy array attributes)
194
195
196
197
198
199

        self.path_MetaPreprocessor = self.path_archive

    def validate_pathes(self):
        if not os.path.isfile(self.path_archive) and not os.path.isdir(self.path_archive):
            self.logger.info("The %s dataset '%s' has not been processed earlier and no corresponding raw data archive"
200
                             "has been found at %s." % (self.sensor, self.entity_ID, self.path_archive))
201
            self.logger.info('Trying to download the dataset...')
202
            self.path_archive_valid = self._data_downloader(self.sensor, self.entity_ID)
203
204
205
        else:
            self.path_archive_valid = True

206
        if not CFG.inmem_serialization and self.ExtractedFolder and not os.path.isdir(self.ExtractedFolder):
207
208
209
210
211
212
213
            os.makedirs(self.ExtractedFolder)

        assert os.path.exists(self.path_archive), 'Invalid path to RAW data. File %s does not exist at %s.' \
                                                  % (os.path.basename(self.path_archive),
                                                     os.path.dirname(self.path_archive))
        assert isinstance(self.path_archive, str), 'Invalid path to RAW data. Got %s instead of string or unicode.' \
                                                   % type(self.path_archive)
214
        if not CFG.inmem_serialization and self.ExtractedFolder:
215
216
            assert os.path.exists(self.path_archive), \
                'Invalid path for temporary files. Directory %s does not exist.' % self.ExtractedFolder
217
218
219
220
221
222
223
224

    @property
    def logger(self):
        if self._loggers_disabled:
            return None
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
225
            self._logger = DatasetLogger('log__' + self.baseN, fmt_suffix=self.scene_ID, path_logfile=self.path_logfile,
226
                                         log_level=CFG.log_level, append=True)
227
228
229
230
            return self._logger

    @logger.setter
    def logger(self, logger):
231
        assert isinstance(logger, logging.Logger) or logger in ['not set', None], \
232
            "GMS_obj.logger can not be set to %s." % logger
233
234

        # save prior logs
235
        # if logger is None and self._logger is not None:
236
        #    self.log += self.logger.captured_stream
237
238
        self._logger = logger

239
240
241
242
243
244
245
246
247
248
249
    @property  # FIXME does not work yet
    def log(self):
        """Returns a string of all logged messages until now."""

        return self._log

    @log.setter
    def log(self, string):
        assert isinstance(string, str), "'log' can only be set to a string. Got %s." % type(string)
        self._log = string

250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
    @property
    def proc_status(self):
        # type: () -> str
        """
        Get the processing status of the current GMS_object (subclass) instance for the current processing level.

        Possible values: 'initialized', 'running', 'finished', 'failed'
        """
        # NOTE: self.proc_status_all_GMSobjs is a class attribute (visible and modifyable from any other subsystem)
        return self.proc_status_all_GMSobjs[self.scene_ID][self.subsystem][self.proc_level]

    @proc_status.setter
    def proc_status(self, new_status):
        # type: (str) -> None
        self.proc_status_all_GMSobjs[self.scene_ID][self.subsystem][self.proc_level] = new_status

266
267
    @property
    def GMS_identifier(self):
268
269
        return GMS_identifier(self.image_type, self.satellite, self.sensor, self.subsystem, self.proc_level,
                              self.dataset_ID, self.logger)
270
271
272

    @property
    def MetaObj(self):
273
274
        # TODO if there is no MetaObj -> create MetaObj by reading metadata from disk
        # reading from disk should use L1A_P.L1A_object.import_metadata -> so just return None
275
276
277
278
279
280

        return self._MetaObj

    @MetaObj.setter
    def MetaObj(self, MetaObj):
        assert isinstance(MetaObj, METADATA), "'MetaObj' can only be set to an instance of METADATA class. " \
281
                                              "Got %s." % type(MetaObj)
282
283
284
285
        self._MetaObj = MetaObj

    @MetaObj.deleter
    def MetaObj(self):
286
287
288
289
290
        if hasattr(self, '_MetaObj') and self._MetaObj and hasattr(self._MetaObj, 'logger') and \
                self._MetaObj.logger not in [None, 'not set']:
            self._MetaObj.logger.close()
            self._MetaObj.logger = None

291
292
        self._MetaObj = None

293
294
295
296
297
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
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
    @property
    def pathGen(self):
        # type: () -> PG.path_generator
        """
        Returns the path generator object for generating file pathes belonging to the GMS object.
        """

        if self._pathGen and self._pathGen.proc_level == self.proc_level:
            return self._pathGen
        else:
            self._pathGen = PG.path_generator(self.__dict__.copy(), MGRS_info=self.MGRS_info)

        return self._pathGen

    @pathGen.setter
    def pathGen(self, pathGen):
        assert isinstance(pathGen, PG.path_generator), 'GMS_object.pathGen can only be set to an instance of ' \
                                                       'path_generator. Got %s.' % type(pathGen)
        self._pathGen = pathGen

    @pathGen.deleter
    def pathGen(self):
        self._pathGen = None

    @property
    def subset(self):
        return [self.arr_shape, self.arr_pos]

    @property
    def LayerBandsAssignment(self):
        # FIXME merge that with self.MetaObj.LayerBandsAssignment
        # FIXME -> otherwise a change of LBA in MetaObj is not recognized here
        if self._LayerBandsAssignment:
            return self._LayerBandsAssignment
        elif self.image_type == 'RSD':
            self._LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier) \
                if self.sensormode != 'P' else get_LayerBandsAssignment(self.GMS_identifier, nBands=1)
            return self._LayerBandsAssignment
        else:
            return ''

    @LayerBandsAssignment.setter
    def LayerBandsAssignment(self, LBA_list):
        self._LayerBandsAssignment = LBA_list
        self.MetaObj.LayerBandsAssignment = LBA_list

339
340
341
342
343
    @property
    def dict_LayerOptTherm(self):
        if self._dict_LayerOptTherm:
            return self._dict_LayerOptTherm
        elif self.LayerBandsAssignment:
344
            self._dict_LayerOptTherm = get_dict_LayerOptTherm(self.GMS_identifier, self.LayerBandsAssignment)
345
346
347
348
349
350
351
            return self._dict_LayerOptTherm
        else:
            return None

    @property
    def georef(self):
        """Returns True if the current dataset can serve as spatial reference."""
Daniel Scheffler's avatar
Daniel Scheffler committed
352

353
354
355
356
        return True if self.image_type == 'RSD' and re.search('OLI', self.sensor, re.I) else False

    @property
    def coreg_needed(self):
357
        if self._coreg_needed is None:
358
            self._coreg_needed = not (self.dataset_ID == CFG.datasetid_spatial_ref)
359
        return self._coreg_needed
360
361
362
363
364

    @coreg_needed.setter
    def coreg_needed(self, value):
        self._coreg_needed = value

365
366
367
368
369
370
    @property
    def coreg_info(self):
        if not self._coreg_info:
            self._coreg_info = {
                'corrected_shifts_px': {'x': 0, 'y': 0},
                'corrected_shifts_map': {'x': 0, 'y': 0},
371
                'original map info': self.MetaObj.map_info,
372
                'updated map info': None,
373
                'shift_reliability': None,
374
375
376
377
                'reference scene ID': None,
                'reference entity ID': None,
                'reference geotransform': None,
                # reference projection must be the own projection in order to avoid overwriting with a wrong EPSG
378
                'reference projection': self.MetaObj.projection,
379
                'reference extent': {'rows': None, 'cols': None},
380
381
                'reference grid': [list(CFG.spatial_ref_gridx),
                                   list(CFG.spatial_ref_gridy)],
382
383
384
385
386
387
388
389
390
                'success': False
            }

        return self._coreg_info

    @coreg_info.setter
    def coreg_info(self, val):
        self._coreg_info = val

391
392
393
    @property
    def resamp_needed(self):
        if self._resamp_needed is None:
394
            gt = mapinfo2geotransform(self.MetaObj.map_info)
395
396
            self._resamp_needed = not is_coord_grid_equal(gt, CFG.spatial_ref_gridx,
                                                          CFG.spatial_ref_gridy)
397
398
399
400
401
402
        return self._resamp_needed

    @resamp_needed.setter
    def resamp_needed(self, value):
        self._resamp_needed = value

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
    @property
    def arr(self):
        # type: () -> GeoArray
        # TODO this must return a subset if self.subset is not None
        return self._arr

    @arr.setter
    def arr(self, *geoArr_initArgs):
        # TODO this must be able to handle subset inputs in tiled processing
        self._arr = GeoArray(*geoArr_initArgs)

        # set nodata value and geoinfos
        # NOTE: MetaObj is NOT gettable before import_metadata has been executed!
        if hasattr(self, 'MetaObj') and self.MetaObj:
            self._arr.nodata = self.MetaObj.spec_vals['fill']
            self._arr.gt = mapinfo2geotransform(self.MetaObj.map_info) if self.MetaObj.map_info else [0, 1, 0, 0, 0, -1]
            self._arr.prj = self.MetaObj.projection if self.MetaObj.projection else self._arr.projection
        else:
            self._arr.nodata = DEF_D.get_outFillZeroSaturated(self._arr.dtype)[0]

        # set bandnames like this: [B01, .., B8A,]
        if self.LayerBandsAssignment:
            if len(self.LayerBandsAssignment) == self._arr.bands:
                self._arr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment)
            else:
                warnings.warn("Cannot set 'bandnames' attribute of GMS_object.arr because LayerBandsAssignment has %s "
                              "bands and GMS_object.arr has %s bands."
                              % (len(self.LayerBandsAssignment), self._arr.bands))

    @arr.deleter
    def arr(self):
        self._arr = None

436
437
438
439
    @property
    def arr_meta(self):
        return self.MetaObj.to_odict()

440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
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
    @property
    def mask_nodata(self):
        if self._mask_nodata is not None:
            if not self._mask_nodata.is_inmem and self._mask_nodata.bands > 1:
                # NoDataMask object self._mask_nodata points to multi-band image file (bands mask_nodata/mask_clouds)
                # -> read processes of not needed bands need to be avoided
                self._mask_nodata = NoDataMask(self._mask_nodata[:, :, 0],
                                               geotransform=self._mask_nodata.gt,
                                               projection=self._mask_nodata.prj)
            return self._mask_nodata

        elif self._masks:
            # return nodata mask from self.masks
            self._mask_nodata = NoDataMask(self._masks[:, :, 0],  # TODO use band names
                                           geotransform=self._masks.gt,
                                           projection=self._masks.prj)
            return self._mask_nodata

        elif isinstance(self.arr, GeoArray):
            self.logger.info('Calculating nodata mask...')
            self._mask_nodata = self.arr.mask_nodata  # calculates mask nodata if not already present
            return self._mask_nodata
        else:
            return None

    @mask_nodata.setter
    def mask_nodata(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            self._mask_nodata = NoDataMask(*geoArr_initArgs)
            self._mask_nodata.nodata = False
            self._mask_nodata.gt = self.arr.gt
            self._mask_nodata.prj = self.arr.prj
        else:
            del self.mask_nodata

    @mask_nodata.deleter
    def mask_nodata(self):
        self._mask_nodata = None

    @property
    def mask_clouds(self):
        if self._mask_clouds is not None:
            if not self._mask_clouds.is_inmem and self._mask_clouds.bands > 1:
                # CloudMask object self._mask_clouds points to multi-band image file on disk
                # (bands mask_nodata/mask_clouds) -> read processes of not needed bands need to be avoided
                self._mask_clouds = CloudMask(self._mask_clouds[:, :, 1],
                                              geotransform=self._mask_clouds.gt,
                                              projection=self._mask_clouds.prj)  # TODO add legend
            return self._mask_clouds

        elif self._masks and self._masks.bands > 1:  # FIXME this will not be secure if there are more than 2 bands
            # return cloud mask from self.masks
            self._mask_clouds = CloudMask(self._masks[:, :, 1],  # TODO use band names
                                          geotransform=self._masks.gt,
                                          projection=self._masks.prj)
            return self._mask_clouds
        else:
            return None  # TODO don't calculate cloud mask?

    @mask_clouds.setter
    def mask_clouds(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:  # FIXME shape validation?
            self._mask_clouds = CloudMask(*geoArr_initArgs)
            self._mask_clouds.nodata = 0
            self._mask_clouds.gt = self.arr.gt
            self._mask_clouds.prj = self.arr.prj
        else:
            del self.mask_clouds

    @mask_clouds.deleter
    def mask_clouds(self):
        self._mask_clouds = None

    @property
    def dem(self):
        """
        Returns an SRTM DEM in the exact dimension an pixel grid of self.arr as an instance of GeoArray.
        """

        if self._dem is None:
            self.logger.info('Generating DEM...')
            DC_args = (self.arr.box.boxMapXY, self.arr.prj, self.arr.xgsd, self.arr.ygsd)
            try:
                DC = INP_R.DEM_Creator(dem_sensor='SRTM', logger=self.logger)
                self._dem = DC.from_extent(*DC_args)
            except RuntimeError:
                self.logger.warning('SRTM DEM generation failed. Falling back to ASTER...')
                DC = INP_R.DEM_Creator(dem_sensor='ASTER', logger=self.logger)
                self._dem = DC.from_extent(*DC_args)

            self._dem.nodata = DEF_D.get_outFillZeroSaturated(self._dem.dtype)[0]
        return self._dem

    @dem.setter
    def dem(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            geoArr = GeoArray(*geoArr_initArgs)
            assert self._dem.shape[:2] == self.arr.shape[:2]

            self._dem = geoArr
            self._dem.nodata = DEF_D.get_outFillZeroSaturated(self._dem.dtype)[0]
            self._dem.gt = self.arr.gt
            self._dem.prj = self.arr.prj
        else:
            del self.dem

    @dem.deleter
    def dem(self):
        self._dem = None

550
551
    @property
    def masks(self):
552
        # if self.mask_nodata is not None and self.mask_clouds is not None and \
553
554
555
556
        #     self._masks is not None and self._masks.bands==1:

        #     self.build_combined_masks_array()

557
558
559
        return self._masks

    @masks.setter
560
    def masks(self, *geoArr_initArgs):
561
562
563
564
        """
        NOTE: This does not automatically update mask_nodata and mask_clouds BUT if mask_nodata and mask_clouds are
        None their getters will automatically synchronize!
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
565

566
        if geoArr_initArgs[0] is not None:
567
            self._masks = GeoArray(*geoArr_initArgs)
568
            self._masks.nodata = 0
569
570
            self._masks.gt = self.arr.gt
            self._masks.prj = self.arr.prj
571
572
        else:
            del self.masks
573

574
575
576
577
    @masks.deleter
    def masks(self):
        self._masks = None

578
579
580
581
582
583
584
585
586
587
588
589
590
    @property
    def mask_clouds_confidence(self):
        return self._mask_clouds_confidence

    @mask_clouds_confidence.setter
    def mask_clouds_confidence(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            cnfArr = GeoArray(*geoArr_initArgs)

            assert cnfArr.shape == self.arr.shape[:2], \
                "The 'mask_clouds_confidence' GeoArray can only be instanced with an array of the same dimensions " \
                "like GMS_obj.arr. Got %s." % str(cnfArr.shape)

Daniel Scheffler's avatar
Daniel Scheffler committed
591
            # noinspection PyProtectedMember
592
593
594
595
            if cnfArr._nodata is None:
                cnfArr.nodata = DEF_D.get_outFillZeroSaturated(cnfArr.dtype)[0]
            cnfArr.gt = self.arr.gt
            cnfArr.prj = self.arr.prj
596
            cnfArr.bandnames = ['confidence']
597
598
599

            self._mask_clouds_confidence = cnfArr
        else:
600
            del self.mask_clouds_confidence
601
602
603
604
605
606
607
608
609
610
611
612

    @mask_clouds_confidence.deleter
    def mask_clouds_confidence(self):
        self._mask_clouds_confidence = None

    @property
    def ac_errors(self):
        """Returns an instance of GeoArray containing error information calculated by the atmospheric correction.

        :return:
        """

613
        return self._ac_errors
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628

    @ac_errors.setter
    def ac_errors(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            errArr = GeoArray(*geoArr_initArgs)

            if CFG.ac_bandwise_accuracy:
                assert errArr.shape == self.arr.shape, \
                    "The 'ac_errors' GeoArray can only be instanced with an array of the same dimensions like " \
                    "GMS_obj.arr. Got %s." % str(errArr.shape)
            else:
                assert errArr.shape[:2] == self.arr.shape[:2], \
                    "The 'ac_errors' GeoArray can only be instanced with an array of the same X/Y dimensions like " \
                    "GMS_obj.arr. Got %s." % str(errArr.shape)

Daniel Scheffler's avatar
Daniel Scheffler committed
629
            # noinspection PyProtectedMember
630
631
632
633
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
634
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
635
636
637
638
639
640
641
642
643

            self._ac_errors = errArr
        else:
            del self.ac_errors

    @ac_errors.deleter
    def ac_errors(self):
        self._ac_errors = None

644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
    @property
    def spat_homo_errors(self):
        if not self._spat_homo_errors and self.coreg_info['shift_reliability'] is not None:
            errArr = GeoArray(np.full(self.arr.shape[:2], self.coreg_info['shift_reliability'], dtype=np.uint8),
                              geotransform=self.arr.geotransform,
                              projection=self.arr.projection,
                              bandnames=['shift_reliability'],
                              nodata=DEF_D.get_outFillZeroSaturated(np.uint8)[0])
            errArr[self.arr.mask_nodata.astype(np.int8) == 0] = errArr.nodata

            self._spat_homo_errors = errArr

        return self._spat_homo_errors

    @spat_homo_errors.deleter
    def spat_homo_errors(self):
        self._spat_homo_errors = None

662
663
664
665
666
667
668
    @property
    def spec_homo_errors(self):
        """Returns an instance of GeoArray containing error information calculated during spectral homogenization.

        :return:
        """

669
        return self._spec_homo_errors
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684

    @spec_homo_errors.setter
    def spec_homo_errors(self, *geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            errArr = GeoArray(*geoArr_initArgs)

            if CFG.spechomo_bandwise_accuracy:
                assert errArr.shape == self.arr.shape, \
                    "The 'spec_homo_errors' GeoArray can only be instanced with an array of the same dimensions like " \
                    "GMS_obj.arr. Got %s." % str(errArr.shape)
            else:
                assert errArr.shape[:2] == self.arr.shape[:2], \
                    "The 'spec_homo_errors' GeoArray can only be instanced with an array of the same X/Y dimensions " \
                    "like GMS_obj.arr. Got %s." % str(errArr.shape)

Daniel Scheffler's avatar
Daniel Scheffler committed
685
            # noinspection PyProtectedMember
686
687
688
689
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
690
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
691
692
693
694
695
696
697
698
699
700
701

            self._spec_homo_errors = errArr
        else:
            del self.spec_homo_errors

    @spec_homo_errors.deleter
    def spec_homo_errors(self):
        self._spec_homo_errors = None

    @property
    def accuracy_layers(self):
702
703
704
705
        if not self._accuracy_layers:
            if not self.proc_level.startswith('L2'):
                self.logger.warning('Attempt to get %s accuracy layers failed - they are a Level 2 feature only.'
                                    % self.proc_level)
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
706
                return None
707

708
709
710
            try:
                from ..algorithms.L2C_P import AccuracyCube
                self._accuracy_layers = AccuracyCube(self)
711
712
713
714

            except ValueError as e:
                if str(e) == 'The given GMS_object contains no accuracy layers for combination.':
                    if CFG.ac_estimate_accuracy or CFG.spechomo_estimate_accuracy:
715
716
                        self.logger.warning('The given GMS_object contains no accuracy layers although computation '
                                            'of accurracy layers was enabled in job configuration.')
717
718
719
720
721
                    else:
                        pass  # self._accuracy_layers keeps None
                else:
                    raise

722
723
            except Exception as e:
                raise RuntimeError('Failed to generate AccuracyCube!', e)
724

725
        return self._accuracy_layers
726
727

    @accuracy_layers.setter
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
728
    def accuracy_layers(self, *geoArr_initArgs):
729
        if geoArr_initArgs[0] is not None:
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
730
            acc_lay = GeoArray(*geoArr_initArgs)
731
732
733
734
            assert acc_lay.shape[:2] == self.arr.shape[:2],\
                "The 'accuracy_layers' GeoArray can only be instanced with an array of the same dimensions like " \
                "GMS_obj.arr. Got %s." % str(acc_lay.shape)

Daniel Scheffler's avatar
Daniel Scheffler committed
735
            # noinspection PyProtectedMember
736
737
738
739
740
741
742
743
744
745
746
747
            if acc_lay._nodata is None:
                acc_lay.nodata = DEF_D.get_outFillZeroSaturated(acc_lay.dtype)[0]
            acc_lay.gt = self.arr.gt
            acc_lay.prj = self.arr.prj

            if not acc_lay.bandnames:
                raise ValueError

            self._accuracy_layers = acc_lay
        else:
            del self._accuracy_layers

748
749
750
751
    @accuracy_layers.deleter
    def accuracy_layers(self):
        self._accuracy_layers = None

752
753
754
755
756
757
758
759
760
761
762
    @property
    def accuracy_layers_meta(self):
        if self._accuracy_layers is not None:
            return {'map info': geotransform2mapinfo(self._accuracy_layers.gt, self._accuracy_layers.projection),
                    'coordinate system string': self._accuracy_layers.projection,
                    'bands': self._accuracy_layers.bands,
                    'band names': list(self._accuracy_layers.bandnames),
                    'data ignore value': self._accuracy_layers.nodata}
        else:
            return None

763
764
765
    @property
    def cloud_masking_algorithm(self):
        if not self._cloud_masking_algorithm:
766
            self._cloud_masking_algorithm = CFG.cloud_masking_algorithm[self.satellite]
767
768
        return self._cloud_masking_algorithm

769
770
    @property
    def ac_options(self):
771
        # type: () -> dict
772
        """
773
774
        Returns the options dictionary needed as input for atmospheric correction. If an empty dictionary is returned,
        atmospheric correction is not yet available for the current sensor and will later be skipped.
775
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
776

777
        if not self._ac_options:
778
            path_ac_options = CFG.path_custom_sicor_options or PG.get_path_ac_options(self.GMS_identifier)
779

780
            if path_ac_options and os.path.exists(path_ac_options):
781
782
                # don't validate because options contain pathes that do not exist on another server:
                opt_dict = get_ac_options(path_ac_options, validation=False)
783

Daniel Scheffler's avatar
Daniel Scheffler committed
784
                # update some file paths depending on the current environment
785
786
                opt_dict['DEM']['fn'] = CFG.path_dem_proc_srtm_90m
                opt_dict['ECMWF']['path_db'] = CFG.path_ECMWF_db
787
788
789
                opt_dict['S2Image'][
                    'S2_MSI_granule_path'] = None  # only a placeholder -> will always be None for GMS usage
                opt_dict['output'] = []  # outputs are not needed for GMS -> so
790
                opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]')
791
                if 'uncertainties' in opt_dict:
792
793
794
795
                    if CFG.ac_estimate_accuracy:
                        opt_dict['uncertainties']['snr_model'] = PG.get_path_snr_model(self.GMS_identifier)
                    else:
                        del opt_dict['uncertainties']  # SICOR will not compute uncertainties if that key is missing
796

797
798
799
800
801
802
                # apply custom configuration
                opt_dict["logger"]['level'] = CFG.log_level
                opt_dict["ram"]['upper_limit'] = CFG.ac_max_ram_gb
                opt_dict["ram"]['unit'] = 'GB'
                opt_dict["AC"]['fill_nonclear_areas'] = CFG.ac_fillnonclear_areas
                opt_dict["AC"]['clear_area_labels'] = CFG.ac_clear_area_labels
803
                # opt_dict['AC']['n_cores'] = CFG.CPUs if CFG.allow_subMultiprocessing else 1
804

805
                self._ac_options = opt_dict
806
807
808
            else:
                self.logger.warning('There is no options file available for atmospheric correction. '
                                    'Atmospheric correction must be skipped.')
809

810
811
        return self._ac_options

812
    def get_copied_dict_and_props(self, remove_privates=False):
813
        # type: (bool) -> dict
814
        """Returns a copy of the current object dictionary including the current values of all object properties."""
815
816

        # loggers must be closed
817
        self.close_loggers()
818
819
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
820
821
822
823
824

        out_dict = self.__dict__.copy()

        # add properties
        property_names = [p for p in dir(self.__class__) if isinstance(getattr(self.__class__, p), property)]
825
        [out_dict.update({propK: copy.copy(getattr(self, propK))}) for propK in property_names]
826
827
828
829
830
831
832
833
834

        # remove private attributes
        if remove_privates:
            out_dict = {k: v for k, v in out_dict.items() if not k.startswith('_')}

        self._loggers_disabled = False  # enables automatic recreation of loggers

        return out_dict

835
836
    def attributes2dict(self, remove_privates=False):
        # type: (bool) -> dict
837
        """Returns a copy of the current object dictionary including the current values of all object properties."""
838
839

        # loggers must be closed
840
        self.close_loggers()
841
842
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
843
844
845
846

        out_dict = self.__dict__.copy()

        # add some selected property values
847
        for i in ['GMS_identifier', 'LayerBandsAssignment', 'coreg_needed', 'coreg_info', 'resamp_needed',
848
849
850
851
852
                  'dict_LayerOptTherm', 'georef', 'MetaObj']:
            if i == 'MetaObj':
                out_dict['meta_odict'] = self.MetaObj.to_odict()
            else:
                out_dict[i] = getattr(self, i)
853
854
855
856
857

        # remove private attributes
        if remove_privates:
            out_dict = {k: v for k, v in out_dict.items() if not k.startswith('_')}

858
        self._loggers_disabled = False  # enables automatic recreation of loggers
859
860
        return out_dict

861
    def _data_downloader(self, sensor, entity_ID):
862
863
864
865
        self.logger.info('Data downloader started.')
        success = False
        " > download source code for Landsat here < "
        if not success:
866
867
            self.logger.critical(
                "Download for %s dataset '%s' failed. No further processing possible." % (sensor, entity_ID))
868
            raise RuntimeError('Archive download failed.')
869
870
        return success

871
872
    @classmethod
    def from_disk(cls, tuple_GMS_subset):
873
874
875
876
        """Fills an already instanced GMS object with data from disk. Excludes array attributes in Python mode.

        :param tuple_GMS_subset:    <tuple> e.g. ('/path/gms_file.gms', ['cube', None])
        """
877
        GMS_obj = cls()
878

879
        path_GMS_file = tuple_GMS_subset[0]
880
        GMSfileDict = INP_R.GMSfile2dict(path_GMS_file)
881
882

        # copy all attributes from GMS file (private attributes are not touched since they are not included in GMS file)
883
884
885
886
887
888
889

        # set MetaObj first in order to make some getters and setters work
        GMS_id = GMS_identifier(image_type=GMSfileDict['image_type'], satellite=GMSfileDict['satellite'],
                                sensor=GMSfileDict['sensor'], subsystem=GMSfileDict['subsystem'],
                                proc_level=GMSfileDict['proc_level'], dataset_ID=GMSfileDict['dataset_ID'], logger=None)
        GMS_obj.MetaObj = METADATA(GMS_id).from_odict(GMSfileDict['meta_odict'])

890
891
        for key, value in GMSfileDict.items():
            if key in ['GMS_identifier', 'georef', 'dict_LayerOptTherm']:
892
                continue  # properties that should better be created on the fly
893
894
895
896
897
898
899
            elif key == 'acq_datetime':
                GMS_obj.acq_datetime = datetime.datetime.strptime(value, '%Y-%m-%d %H:%M:%S.%f%z')
            else:
                try:
                    setattr(GMS_obj, key, value)
                except Exception:
                    raise AttributeError("Can't set attribute %s." % key)
900

901
        GMS_obj.arr_shape, GMS_obj.arr_pos = tuple_GMS_subset[1]
902

903
        GMS_obj.arr = GMS_obj.pathGen.get_path_imagedata()
904
        # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters):
905
        GMS_obj.masks = GMS_obj.pathGen.get_path_maskdata()
906
        GMS_obj.accuracy_layers = GMS_obj.pathGen.get_path_accuracylayers()
907

908
        return GMS_obj
909

910
911
    @classmethod
    def from_sensor_subsystems(cls, list_GMS_objs):
912
        # type: (List[GMS_object]) -> GMS_object
913
914
915
916
917
        """Merge separate GMS objects belonging to the same scene-ID into ONE GMS object.

        :param list_GMS_objs:   <list> of GMS objects covering the same geographic area but representing different
                                sensor subsystems (e.g. 3 GMS_objects for Sentinel-2 10m/20m/60m bands)
        """
918
        GMS_obj_merged = cls()
919

920
        # assertions
921
922
        assert len(list_GMS_objs) > 1, "'GMS_object.from_sensor_subsystems()' expects multiple input GMS objects. " \
                                       "Got %d." % len(list_GMS_objs)
923
        assert all([is_coord_grid_equal(list_GMS_objs[0].arr.gt, *obj.arr.xygrid_specs) for obj in list_GMS_objs[1:]]),\
924
925
926
            "The input GMS objects must have the same pixel grid. Received: %s" \
            % np.array([obj.arr.xygrid_specs for obj in list_GMS_objs])
        assert len(list(set([GMS_obj.proc_level for GMS_obj in list_GMS_objs]))) == 1, \
927
928
929
            "The input GMS objects for GMS_object.from_sensor_subsystems() must have the same processing level."
        subsystems = [GMS_obj.subsystem for GMS_obj in list_GMS_objs]
        assert len(subsystems) == len(list(set(subsystems))), \
930
            "The input 'list_GMS_objs' contains duplicates: %s" % subsystems
931

932
933
934
935
936
        ##################
        # merge logfiles #
        ##################

        # read all logs into DataFrame, sort it by the first column
937
        [GMS_obj.close_loggers() for GMS_obj in list_GMS_objs]  # close the loggers of the input objects
938
939
940
941
942
        paths_inLogs = [GMS_obj.pathGen.get_path_logfile() for GMS_obj in list_GMS_objs]
        allLogs_df = DataFrame()
        for log in paths_inLogs:
            df = read_csv(log, sep='\n', delimiter=':   ', header=None,
                          engine='python')  # engine suppresses a pandas warning
943
944
            # FIXME this will log e.g. atm. corr 3 times for S2A -> use captured streams instead?
            allLogs_df = allLogs_df.append(df)
945
946
947
948

        allLogs_df = allLogs_df.sort_values(0)

        # set common metadata, needed for logfile
949
950
951
        GMS_obj_merged.baseN = list_GMS_objs[0].pathGen.get_baseN(merged_subsystems=True)
        GMS_obj_merged.path_logfile = list_GMS_objs[0].pathGen.get_path_logfile(merged_subsystems=True)
        GMS_obj_merged.scene_ID = list_GMS_objs[0].scene_ID
952
953

        # write the merged logfile and flush previous logger
954
        np.savetxt(GMS_obj_merged.path_logfile, np.array(allLogs_df), delimiter=':   ', fmt="%s")
955
        GMS_obj_merged.close_loggers()
956

957
        # log
958
959
        GMS_obj_merged.logger.info('Merging the subsystems %s to a single GMS object...'
                                   % ', '.join([GMS_obj.subsystem for GMS_obj in list_GMS_objs]))
960

961
962
963
964
        ##################
        # MERGE METADATA #
        ##################

965
966
967
        # copy all attributes from the first input GMS file (private attributes are not touched)
        for key, value in list_GMS_objs[0].__dict__.copy().items():
            if key in ['GMS_identifier', 'georef', 'dict_LayerOptTherm']:
968
                continue  # properties that should better be created on the fly
969
970
            elif key in ['baseN', 'path_logfile', 'scene_ID', 'subsystem']:
                continue  # either previously set with common values or not needed for merged GMS_object
971
            try:
972
                setattr(GMS_obj_merged, key, value)
973
974
            except Exception:
                raise AttributeError("Can't set attribute %s." % key)
975

976
        # update LayerBandsAssignment and get full list of output bandnames
977
        from .metadata import get_LayerBandsAssignment
978
979
        # use identifier of first input GMS object for getting LBA (respects current proc_level):
        gms_idf = list_GMS_objs[0].GMS_identifier
980
981
        GMS_obj_merged.LayerBandsAssignment = get_LayerBandsAssignment(gms_idf, return_fullLBA=True)
        bandnames = ['B%s' % i if len(i) == 2 else 'B0%s' % i for i in GMS_obj_merged.LayerBandsAssignment]
982
983

        # update layer-dependent metadata with respect to remaining input GMS objects
984
985
986
        GMS_obj_merged.MetaObj.LayerBandsAssignment = GMS_obj_merged.LayerBandsAssignment
        GMS_obj_merged.MetaObj.Subsystem = ''

987
988
        GMS_obj_merged.subsystem = ''
        del GMS_obj_merged.pathGen  # must be refreshed because subsystem is now ''
989
        GMS_obj_merged.close_loggers()  # must also be refreshed because it depends on pathGen
990

991
        for attrN in layerdependent_metadata:
992
993
994
995
996
997
            # combine values from separate subsystems to a single value
            attrDic_fullLBA = {}
            for GMS_obj in list_GMS_objs:
                attr_val = getattr(GMS_obj.MetaObj, attrN)
                if isinstance(attr_val, list):
                    attrDic_fullLBA.update(dict(zip(GMS_obj.LayerBandsAssignment, attr_val)))
998
                elif isinstance(attr_val, (dict, collections.OrderedDict)):
999
1000
1001
1002
1003
1004
                    attrDic_fullLBA.update(attr_val)
                else:
                    raise ValueError(attrN)

            # update the attribute in self.MetaObj
            if attrDic_fullLBA:
1005
                val2set = [attrDic_fullLBA[bN] for bN in GMS_obj_merged.LayerBandsAssignment] \
1006
                    if isinstance(getattr(list_GMS_objs[0].MetaObj, attrN), list) else attrDic_fullLBA
1007
                setattr(GMS_obj_merged.MetaObj, attrN, val2set)
1008

1009
1010
1011
        ####################
        # MERGE ARRAY DATA #
        ####################
1012

1013
1014
1015
1016
        # find the common extent. NOTE: boundsMap is expected in the order [xmin,xmax,ymin,ymax]
        geoExtents = np.array([GMS_obj.arr.box.boundsMap for GMS_obj in list_GMS_objs])
        common_extent = (min(geoExtents[:, 0]), max(geoExtents[:, 1]), min(geoExtents[:, 2]), max(geoExtents[:, 3]))

1017
        # overwrite array data with merged arrays, clipped to common_extent and reordered according to FullLayerBandsAss
1018
1019
1020
        for attrname in ['arr', 'ac_errors', 'dem', 'mask_nodata', 'mask_clouds', 'mask_clouds_confidence', 'masks']:

            # get current attribute of each subsystem without running property getters
1021
            all_arrays = [getattr(GMS_obj, '_%s' % attrname) for GMS_obj in list_GMS_objs]
1022
1023
1024
1025
1026
1027
1028
1029
1030

            # get the same geographical extent for each input GMS object
            if len(set(tuple(ext) for ext in geoExtents.tolist())) > 1:
                # in case of different extents
                geoArrs_same_extent = []

                for geoArr in all_arrays:

                    if geoArr is not None:
1031
                        # FIXME mask_clouds_confidence is no GeoArray until here
1032
                        # FIXME -> has no nodata value -> calculation throughs warning
1033
1034
                        geoArr_same_extent = \
                            GeoArray(*geoArr.get_mapPos(
1035
1036
1037
1038
                                mapBounds=np.array(common_extent)[[0, 2, 1, 3]],  # pass (xmin, ymin, xmax, ymax)
                                mapBounds_prj=geoArr.prj),
                                     bandnames=list(geoArr.bandnames.keys()),
                                     nodata=geoArr.nodata)
1039
1040
                        geoArrs_same_extent.append(geoArr_same_extent)
                    else:
1041
1042
                        # e.g. in case of cloud mask that is only appended to the GMS object with the same
                        # spatial resolution)
1043
1044
1045
1046
1047
1048
                        geoArrs_same_extent.append(None)

            else:
                # skip get_mapPos() if all input GMS objects have the same extent
                geoArrs_same_extent = all_arrays

1049
1050
            # validate output GeoArrays #
            #############################
1051

1052
1053
            if len([gA for gA in geoArrs_same_extent if gA is not None]) > 1:
                equal_bounds = all([geoArrs_same_extent[0].box.boundsMap == gA.box.boundsMap
1054
                                    for gA in geoArrs_same_extent[1:]])
1055
1056
                equal_epsg = all([geoArrs_same_extent[0].epsg == gA.epsg for gA in geoArrs_same_extent[1:]])
                equal_xydims = all([geoArrs_same_extent[0].shape[:2] == gA.shape[:2] for gA in geoArrs_same_extent[1:]])
1057
1058
1059
1060
1061
                if not all([equal_bounds, equal_epsg, equal_xydims]):
                    raise RuntimeError('Something went wrong during getting the same geographical extent for all the '
                                       'input GMS objects. The extents, projections or pixel dimensions of the '
                                       'calculated input GMS objects are not equal.')

1062
1063
            # set output arrays #
            #####################
1064

1065
1066
            # handle those arrays where bands have to be reordered according to FullLayerBandsAssignment
            if attrname in ['arr', 'ac_errors'] and list(set(geoArrs_same_extent)) != [None]:
1067
1068
                # check that each desired band name for the current attribute is provided by one of the input
                # GMS objects
1069
1070
                available_bandNs = list(chain.from_iterable([list(gA.bandnames) for gA in geoArrs_same_extent]))
                for bN in bandnames:
1071
                    if bN not in available_bandNs:
1072
                        raise ValueError("The given input GMS objects (subsystems) do not provide a bandname '%s' for "
1073
1074
                                         "the attribute '%s'. Available band names amongst all input GMS objects are: "
                                         "%s" % (bN, attrname, str(available_bandNs)))
1075
1076

                # merge arrays
1077
1078
                def get_band(bandN):
                    return [gA[bandN] for gA in geoArrs_same_extent if gA and bandN in gA.bandnames][0]
1079
1080
                full_geoArr = GeoArray(np.dstack((get_band(bandN) for bandN in bandnames)),
                                       geoArrs_same_extent[0].gt, geoArrs_same_extent[0].prj,
1081
1082
                                       bandnames=bandnames,
                                       nodata=geoArrs_same_extent[0].nodata)
1083
                setattr(GMS_obj_merged, attrname, full_geoArr)
1084

1085
            # handle the remaining arrays
1086
1087
1088
            else:
                # masks, dem, mask_nodata, mask_clouds, mask_clouds_confidence
                if attrname == 'dem':
1089
1090
                    # use the DEM of the first input object
                    # (if the grid is the same, the DEMs should be the same anyway)
1091
                    GMS_obj_merged.dem = geoArrs_same_extent[0]
1092

1093
1094
                elif attrname == 'mask_nodata':
                    # must not be merged -> self.arr is already merged, so just recalculate it (np.all)
1095
                    GMS_obj_merged.mask_nodata = GMS_obj_merged.calc_mask_nodata(overwrite=True)
1096

1097
1098
1099
                elif attrname == 'mask_clouds':
                    # possibly only present in ONE subsystem (set by atm. Corr.)
                    mask_clouds = [msk for msk in geoArrs_same_extent if msk is not None]
1100
1101
                    if len(mask_clouds) > 1:
                        raise ValueError('Expected mask clouds in only one subsystem. Got %s.' % len(mask_clouds))
1102
                    GMS_obj_merged.mask_clouds = mask_clouds[0] if mask_clouds else None
1103

1104
1105
1106
                elif attrname == 'mask_clouds_confidence':
                    # possibly only present in ONE subsystem (set by atm. Corr.)
                    mask_clouds_conf = [msk for msk in geoArrs_same_extent if msk is not None]
1107
1108
1109
                    if len(mask_clouds_conf) > 1:
                        raise ValueError(
                            'Expected mask_clouds_conf in only one subsystem. Got %s.' % len(mask_clouds_conf))
1110
                    GMS_obj_merged.mask_clouds_confidence = mask_clouds_conf[0] if mask_clouds_conf else None
1111

1112
                elif attrname == 'masks':
1113
                    # self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
1114
                    GMS_obj_merged.masks = None
1115

1116
1117
1118
1119
1120
1121
1122
1123
        # avoid unequal nodata edges between individual layers (resampling artifacts etc.) #
        ####################################################################################

        # get pixel values of areas that have not been atmospherically corrected (non-clear)
        nonclear_labels = [lbl for lbl in ["Clear", "Snow", "Water", "Shadow", "Cirrus", "Cloud"]
                           if lbl not in CFG.ac_clear_area_labels]
        cloud_mask_legend = DEF_D.get_mask_classdefinition('mask_clouds', GMS_obj_merged.satellite)
        nonclear_pixVals = [cloud_mask_legend[lbl] for lbl in nonclear_labels]
1124
1125
1126
1127

        # apply cloud mask to image data and all products derived from image data
        # (only if image data represents BOA-Ref and cloud areas are not to be filled with TOA-Ref)

1128
        if re.search('BOA_Reflectance', GMS_obj_merged.MetaObj.PhysUnit, re.I) and not CFG.ac_fillnonclear_areas:
1129
1130
            # fill non-clear areas with no data values (for all bands)
            for pixVal in nonclear_pixVals:
1131
1132
                mask_nonclear = GMS_obj_merged.mask_clouds[:] == pixVal
                GMS_obj_merged.arr[mask_nonclear] = DEF_D.get_outFillZeroSaturated(GMS_obj_merged.arr.dtype)[0]
1133

1134
1135
1136
                if GMS_obj_merged.ac_errors:
                    GMS_obj_merged.ac_errors[mask_nonclear] = \
                        DEF_D.get_outFillZeroSaturated(GMS_obj_merged.ac_errors.dtype)[0]
1137

1138
        # update no data mask (adds all pixels to no data where at least one band has no data)