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
88
        self._cloud_masking_algorithm = None
        self._meta_odict = None
89
        self._coreg_info = None
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
130
        # 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

131
        self.job_ID = CFG.ID
132
133
134
        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])
135
136
        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
137
        self.MGRS_info = None
Daniel Scheffler's avatar
Daniel Scheffler committed
138
139
        self.lonlat_arr = None  # set by self.write_tiles_to_ENVIfile
        self.trueDataCornerUTM = None  # set by self.from_tiles
140
141

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

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

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

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

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

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

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    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

181
    def set_pathes(self):
182
183
184
185
186
187
        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()
188

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

        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"
201
                             "has been found at %s." % (self.sensor, self.entity_ID, self.path_archive))
202
            self.logger.info('Trying to download the dataset...')
203
            self.path_archive_valid = self._data_downloader(self.sensor, self.entity_ID)
204
205
206
        else:
            self.path_archive_valid = True

207
        if not CFG.inmem_serialization and self.ExtractedFolder and not os.path.isdir(self.ExtractedFolder):
208
209
210
211
212
213
214
            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)
215
        if not CFG.inmem_serialization and self.ExtractedFolder:
216
217
            assert os.path.exists(self.path_archive), \
                'Invalid path for temporary files. Directory %s does not exist.' % self.ExtractedFolder
218
219
220
221
222
223
224
225

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

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

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

240
241
242
243
244
245
246
247
248
249
250
    @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

251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
    @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

267
268
    @property
    def GMS_identifier(self):
269
270
        return collections.OrderedDict(zip(
            ['image_type', 'Satellite', 'Sensor', 'Subsystem', 'proc_level', 'dataset_ID', 'logger'],
271
272
            [self.image_type, self.satellite, self.sensor, self.subsystem, self.proc_level, self.dataset_ID,
             self.logger]))
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288

    @property
    def MetaObj(self):
        if self._meta_odict:
            # if there is already a meta_odict -> create a new MetaObj from it (ensures synchronization!)
            self._MetaObj = METADATA(self.GMS_identifier).from_odict(self._meta_odict)
            del self.meta_odict
        elif not self._MetaObj:
            # if there is no meta_odict and no MetaObj -> create MetaObj by reading metadata from disk
            pass  # reading from disk should use L1A_P.L1A_object.import_metadata -> so just return None

        return self._MetaObj

    @MetaObj.setter
    def MetaObj(self, MetaObj):
        assert isinstance(MetaObj, METADATA), "'MetaObj' can only be set to an instance of METADATA class. " \
289
                                              "Got %s." % type(MetaObj)
290
291
292
        self._MetaObj = MetaObj

        # update meta_odict
293
        del self.meta_odict  # it is recreated if getter is used the next time
294
295
296

    @MetaObj.deleter
    def MetaObj(self):
297
298
299
300
301
        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

302
303
304
305
306
        self._MetaObj = None

    @property
    def meta_odict(self):
        if self._MetaObj:
307
            # if there is already a MetaObj -> create new meta_odict from it (ensures synchronization!)
308
309
310
311
            self._meta_odict = self._MetaObj.to_odict()
            del self.MetaObj
        elif not self._meta_odict:
            # if there is no MetaObj and no meta_odict -> use MetaObj getter to read metadata from disk
312
            pass  # reading from disk should use L1A_P.L1A_object.import_metadata -> so just return None
313
314
315
316
317
318
            self._meta_odict = None

        return self._meta_odict

    @meta_odict.setter
    def meta_odict(self, odict):
319
320
        assert isinstance(odict, (collections.OrderedDict, dict)), "'meta_odict' can only be set to an instance of " \
                                                                   "collections.OrderedDict. Got %s." % type(odict)
321
322
323
        self._meta_odict = odict

        # update MetaObj
324
        del self.MetaObj  # it is recreated if getter is used the next time
325
326
327
328
329

    @meta_odict.deleter
    def meta_odict(self):
        self._meta_odict = None

330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    @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

376
377
378
379
380
    @property
    def dict_LayerOptTherm(self):
        if self._dict_LayerOptTherm:
            return self._dict_LayerOptTherm
        elif self.LayerBandsAssignment:
381
            self._dict_LayerOptTherm = get_dict_LayerOptTherm(self.GMS_identifier, self.LayerBandsAssignment)
382
383
384
385
386
387
388
            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
389

390
391
392
393
        return True if self.image_type == 'RSD' and re.search('OLI', self.sensor, re.I) else False

    @property
    def coreg_needed(self):
394
        if self._coreg_needed is None:
395
            self._coreg_needed = not (self.dataset_ID == CFG.datasetid_spatial_ref)
396
        return self._coreg_needed
397
398
399
400
401

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

402
403
404
405
406
407
408
409
    @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},
                'original map info': self.meta_odict['map info'],
                'updated map info': None,
410
                'shift_reliability': None,
411
412
413
414
415
416
                '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
                'reference projection': self.meta_odict['coordinate system string'],
                'reference extent': {'rows': None, 'cols': None},
417
418
                'reference grid': [list(CFG.spatial_ref_gridx),
                                   list(CFG.spatial_ref_gridy)],
419
420
421
422
423
424
425
426
427
                'success': False
            }

        return self._coreg_info

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

428
429
430
431
    @property
    def resamp_needed(self):
        if self._resamp_needed is None:
            gt = mapinfo2geotransform(self.meta_odict['map info'])
432
433
            self._resamp_needed = not is_coord_grid_equal(gt, CFG.spatial_ref_gridx,
                                                          CFG.spatial_ref_gridy)
434
435
436
437
438
439
        return self._resamp_needed

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

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
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
583
584
585
    @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]
            if hasattr(self, 'meta_odict') and self.meta_odict:
                self._arr.gt = mapinfo2geotransform(self.meta_odict['map info'])
                self._arr.prj = self.meta_odict['coordinate system string']

        # 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

    @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

586
587
    @property
    def masks(self):
588
        # if self.mask_nodata is not None and self.mask_clouds is not None and \
589
590
591
592
        #     self._masks is not None and self._masks.bands==1:

        #     self.build_combined_masks_array()

593
594
595
        return self._masks

    @masks.setter
596
    def masks(self, *geoArr_initArgs):
597
598
599
600
        """
        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
601

602
        if geoArr_initArgs[0] is not None:
603
            self._masks = GeoArray(*geoArr_initArgs)
604
            self._masks.nodata = 0
605
606
            self._masks.gt = self.arr.gt
            self._masks.prj = self.arr.prj
607
608
        else:
            del self.masks
609

610
611
612
613
    @masks.deleter
    def masks(self):
        self._masks = None

614
615
616
617
618
619
620
621
622
623
624
625
626
    @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
627
            # noinspection PyProtectedMember
628
629
630
631
            if cnfArr._nodata is None:
                cnfArr.nodata = DEF_D.get_outFillZeroSaturated(cnfArr.dtype)[0]
            cnfArr.gt = self.arr.gt
            cnfArr.prj = self.arr.prj
632
            cnfArr.bandnames = ['confidence']
633
634
635

            self._mask_clouds_confidence = cnfArr
        else:
636
            del self.mask_clouds_confidence
637
638
639
640
641
642
643
644
645
646
647
648

    @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:
        """

649
        return self._ac_errors
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664

    @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
665
            # noinspection PyProtectedMember
666
667
668
669
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
670
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
671
672
673
674
675
676
677
678
679

            self._ac_errors = errArr
        else:
            del self.ac_errors

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

680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
    @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

698
699
700
701
702
703
704
    @property
    def spec_homo_errors(self):
        """Returns an instance of GeoArray containing error information calculated during spectral homogenization.

        :return:
        """

705
        return self._spec_homo_errors
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720

    @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
721
            # noinspection PyProtectedMember
722
723
724
725
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
726
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
727
728
729
730
731
732
733
734
735
736
737

            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):
738
739
740
741
        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
742
                return None
743

744
745
746
747
            self.logger.info('Generating combined accuracy layers array..')
            try:
                from ..algorithms.L2C_P import AccuracyCube
                self._accuracy_layers = AccuracyCube(self)
748
749
750
751

            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:
752
753
                        self.logger.warning('The given GMS_object contains no accuracy layers although computation '
                                            'of accurracy layers was enabled in job configuration.')
754
755
756
757
758
                    else:
                        pass  # self._accuracy_layers keeps None
                else:
                    raise

759
760
            except Exception as e:
                raise RuntimeError('Failed to generate AccuracyCube!', e)
761

762
        return self._accuracy_layers
763
764
765
766
767
768
769
770
771

    @accuracy_layers.setter
    def accuracy_layers(self, geoArr_initArgs):
        if geoArr_initArgs[0] is not None:
            acc_lay = GeoArray(geoArr_initArgs)
            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
772
            # noinspection PyProtectedMember
773
774
775
776
777
778
779
780
781
782
783
784
            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

785
786
787
788
    @accuracy_layers.deleter
    def accuracy_layers(self):
        self._accuracy_layers = None

789
790
791
792
793
794
795
796
797
798
799
    @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

800
801
802
    @property
    def cloud_masking_algorithm(self):
        if not self._cloud_masking_algorithm:
803
            self._cloud_masking_algorithm = CFG.cloud_masking_algorithm[self.satellite]
804
805
        return self._cloud_masking_algorithm

806
807
    @property
    def ac_options(self):
808
        # type: () -> dict
809
        """
810
811
        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.
812
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
813

814
        if not self._ac_options:
815
            path_ac_options = CFG.path_custom_sicor_options or PG.get_path_ac_options(self.GMS_identifier)
816

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

Daniel Scheffler's avatar
Daniel Scheffler committed
821
                # update some file paths depending on the current environment
822
823
                opt_dict['DEM']['fn'] = CFG.path_dem_proc_srtm_90m
                opt_dict['ECMWF']['path_db'] = CFG.path_ECMWF_db
824
825
826
                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
827
                opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]')
828
                if 'uncertainties' in opt_dict:
829
830
831
832
                    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
833

834
835
836
837
838
839
                # 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
840
                # opt_dict['AC']['n_cores'] = CFG.CPUs if CFG.allow_subMultiprocessing else 1
841

842
                self._ac_options = opt_dict
843
844
845
            else:
                self.logger.warning('There is no options file available for atmospheric correction. '
                                    'Atmospheric correction must be skipped.')
846

847
848
        return self._ac_options

849
    def get_copied_dict_and_props(self, remove_privates=False):
850
        # type: (bool) -> dict
851
        """Returns a copy of the current object dictionary including the current values of all object properties."""
852
853

        # loggers must be closed
854
        self.close_loggers()
855
856
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
857
858
859
860
861

        out_dict = self.__dict__.copy()

        # add properties
        property_names = [p for p in dir(self.__class__) if isinstance(getattr(self.__class__, p), property)]
862
        [out_dict.update({propK: copy.copy(getattr(self, propK))}) for propK in property_names]
863
864
865
866
867
868
869
870
871

        # 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

872
873
    def attributes2dict(self, remove_privates=False):
        # type: (bool) -> dict
874
        """Returns a copy of the current object dictionary including the current values of all object properties."""
875
876

        # loggers must be closed
877
        self.close_loggers()
878
879
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
880
881
882
883

        out_dict = self.__dict__.copy()

        # add some selected property values
884
885
        for i in ['GMS_identifier', 'LayerBandsAssignment', 'coreg_needed', 'coreg_info', 'resamp_needed',
                  'dict_LayerOptTherm', 'georef', 'meta_odict']:
886
            out_dict[i] = getattr(self, i)
887
888
889
890
891

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

892
        self._loggers_disabled = False  # enables automatic recreation of loggers
893
894
        return out_dict

895
    def _data_downloader(self, sensor, entity_ID):
896
897
898
899
        self.logger.info('Data downloader started.')
        success = False
        " > download source code for Landsat here < "
        if not success:
900
901
            self.logger.critical(
                "Download for %s dataset '%s' failed. No further processing possible." % (sensor, entity_ID))
902
            raise RuntimeError('Archive download failed.')
903
904
        return success

905
906
    @classmethod
    def from_disk(cls, tuple_GMS_subset):
907
908
909
910
        """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])
        """
911
        GMS_obj = cls()
912

913
        path_GMS_file = tuple_GMS_subset[0]
914
        GMSfileDict = INP_R.GMSfile2dict(path_GMS_file)
915
916

        # copy all attributes from GMS file (private attributes are not touched since they are not included in GMS file)
917
        GMS_obj.meta_odict = GMSfileDict['meta_odict']  # set that first in order to make some getters and setters work
918
919
        for key, value in GMSfileDict.items():
            if key in ['GMS_identifier', 'georef', 'dict_LayerOptTherm']:
920
                continue  # properties that should better be created on the fly
921
922
923
924
925
926
927
            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)
928

929
        GMS_obj.arr_shape, GMS_obj.arr_pos = tuple_GMS_subset[1]
930

931
        GMS_obj.arr = GMS_obj.pathGen.get_path_imagedata()
932
        # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters):
933
        GMS_obj.masks = GMS_obj.pathGen.get_path_maskdata()
934
        GMS_obj.accuracy_layers = GMS_obj.pathGen.get_path_accuracylayers()
935

936
        return GMS_obj
937

938
939
    @classmethod
    def from_sensor_subsystems(cls, list_GMS_objs):
940
        # type: (List[GMS_object]) -> GMS_object
941
942
943
944
945
        """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)
        """
946
        GMS_obj_merged = cls()
947

948
        # assertions
949
950
        assert len(list_GMS_objs) > 1, "'GMS_object.from_sensor_subsystems()' expects multiple input GMS objects. " \
                                       "Got %d." % len(list_GMS_objs)
951
        assert all([is_coord_grid_equal(list_GMS_objs[0].arr.gt, *obj.arr.xygrid_specs) for obj in list_GMS_objs[1:]]),\
952
953
954
            "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, \
955
956
957
            "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))), \
958
            "The input 'list_GMS_objs' contains duplicates: %s" % subsystems
959

960
961
962
963
964
        ##################
        # merge logfiles #
        ##################

        # read all logs into DataFrame, sort it by the first column
965
        [GMS_obj.close_loggers() for GMS_obj in list_GMS_objs]  # close the loggers of the input objects
966
967
968
969
970
        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
971
972
            # FIXME this will log e.g. atm. corr 3 times for S2A -> use captured streams instead?
            allLogs_df = allLogs_df.append(df)
973
974
975
976

        allLogs_df = allLogs_df.sort_values(0)

        # set common metadata, needed for logfile
977
978
979
        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
980
981

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

985
        # log
986
987
        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]))
988

989
990
991
992
        ##################
        # MERGE METADATA #
        ##################

993
994
995
        # 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']:
996
                continue  # properties that should better be created on the fly
997
998
            elif key in ['baseN', 'path_logfile', 'scene_ID', 'subsystem']:
                continue  # either previously set with common values or not needed for merged GMS_object
999
            try:
1000
                setattr(GMS_obj_merged, key, value)
1001
1002
            except Exception:
                raise AttributeError("Can't set attribute %s." % key)
1003

1004
        # update LayerBandsAssignment and get full list of output bandnames
1005
        from .metadata import get_LayerBandsAssignment
1006
1007
        # use identifier of first input GMS object for getting LBA (respects current proc_level):
        gms_idf = list_GMS_objs[0].GMS_identifier
1008
1009
        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]
1010
1011

        # update layer-dependent metadata with respect to remaining input GMS objects
1012
1013
1014
        GMS_obj_merged.meta_odict.update({
            'band names': [('Band %s' % i) for i in GMS_obj_merged.LayerBandsAssignment],
            'LayerBandsAssignment': GMS_obj_merged.LayerBandsAssignment,
1015
            'Subsystem': '',
1016
            'PhysUnit': GMS_obj_merged.meta_odict['PhysUnit'],  # TODO can contain val_optical / val_thermal
1017
        })
1018
1019
        GMS_obj_merged.subsystem = ''
        del GMS_obj_merged.pathGen  # must be refreshed because subsystem is now ''
1020
        GMS_obj_merged.close_loggers()  # must also be refreshed because it depends on pathGen
1021

1022
        for attrN in layerdependent_metadata:
1023
1024
1025
1026
1027
1028
            # 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)))
1029
                elif isinstance(attr_val, (dict, collections.OrderedDict)):
1030
1031
1032
1033
1034
1035
                    attrDic_fullLBA.update(attr_val)
                else:
                    raise ValueError(attrN)

            # update the attribute in self.MetaObj
            if attrDic_fullLBA:
1036
                val2set = [attrDic_fullLBA[bN] for bN in GMS_obj_merged.LayerBandsAssignment] \
1037
                    if isinstance(getattr(list_GMS_objs[0].MetaObj, attrN), list) else attrDic_fullLBA
1038
                setattr(GMS_obj_merged.MetaObj, attrN, val2set)
1039

1040
1041
1042
        ####################
        # MERGE ARRAY DATA #
        ####################
1043

1044
1045
1046
1047
        # 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]))

1048
        # overwrite array data with merged arrays, clipped to common_extent and reordered according to FullLayerBandsAss
1049
1050
1051
        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
1052
            all_arrays = [getattr(GMS_obj, '_%s' % attrname) for GMS_obj in list_GMS_objs]
1053
1054
1055
1056
1057
1058
1059
1060
1061

            # 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:
1062
                        # FIXME mask_clouds_confidence is no GeoArray until here
1063
                        # FIXME -> has no nodata value -> calculation throughs warning
1064
1065
                        geoArr_same_extent = \
                            GeoArray(*geoArr.get_mapPos(
1066
1067
1068
1069
                                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)
1070
1071
                        geoArrs_same_extent.append(geoArr_same_extent)
                    else:
1072
1073
                        # e.g. in case of cloud mask that is only appended to the GMS object with the same
                        # spatial resolution)
1074
1075
1076
1077
1078
1079
                        geoArrs_same_extent.append(None)

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

1080
1081
            # validate output GeoArrays #
            #############################
1082

1083
1084
            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
1085
                                    for gA in geoArrs_same_extent[1:]])
1086
1087
                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:]])
1088
1089
1090
1091
1092
                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.')

1093
1094
            # set output arrays #
            #####################
1095

1096
1097
            # 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]:
1098
1099
                # check that each desired band name for the current attribute is provided by one of the input
                # GMS objects
1100
1101
                available_bandNs = list(chain.from_iterable([list(gA.bandnames) for gA in geoArrs_same_extent]))
                for bN in bandnames:
1102
                    if bN not in available_bandNs:
1103
                        raise ValueError("The given input GMS objects (subsystems) do not provide a bandname '%s' for "
1104
1105
                                         "the attribute '%s'. Available band names amongst all input GMS objects are: "
                                         "%s" % (bN, attrname, str(available_bandNs)))
1106
1107

                # merge arrays
1108
1109
                def get_band(bandN):
                    return [gA[bandN] for gA in geoArrs_same_extent if gA and bandN in gA.bandnames][0]
1110
1111
                full_geoArr = GeoArray(np.dstack((get_band(bandN) for bandN in bandnames)),
                                       geoArrs_same_extent[0].gt, geoArrs_same_extent[0].prj,
1112
1113
                                       bandnames=bandnames,
                                       nodata=geoArrs_same_extent[0].nodata)
1114
                setattr(GMS_obj_merged, attrname, full_geoArr)
1115

1116
            # handle the remaining arrays
1117
1118
1119
            else:
                # masks, dem, mask_nodata, mask_clouds, mask_clouds_confidence
                if attrname == 'dem':
1120
1121
                    # use the DEM of the first input object
                    # (if the grid is the same, the DEMs should be the same anyway)
1122
                    GMS_obj_merged.dem = geoArrs_same_extent[0]
1123

1124
1125
                elif attrname == 'mask_nodata':
                    # must not be merged -> self.arr is already merged, so just recalculate it (np.all)
1126
                    GMS_obj_merged.mask_nodata = GMS_obj_merged.calc_mask_nodata(overwrite=True)
1127

1128
1129
1130
                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]
1131
1132
                    if len(mask_clouds) > 1:
                        raise ValueError('Expected mask clouds in only one subsystem. Got %s.' % len(mask_clouds))
1133
                    GMS_obj_merged.mask_clouds = mask_clouds[0] if mask_clouds else None
1134

1135
1136
1137
                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]
1138
1139
1140
                    if len(mask_clouds_conf) > 1:
                        raise ValueError(
                            'Expected mask_clouds_conf in only one subsystem. Got %s.' % len(mask_clouds_conf))
1141
                    GMS_obj_merged.mask_clouds_confidence = mask_clouds_conf[0] if mask_clouds_conf else None
1142

1143
                elif attrname == 'masks':
1144
                    # self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
1145
                    GMS_obj_merged.masks = None
1146

1147
1148
1149
1150
1151
1152
1153
1154
        # 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]
1155
1156
1157
1158

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

1159
        if re.search('BOA_Reflectance', GMS_obj_merged.MetaObj.PhysUnit, re.I) and not CFG.ac_fillnonclear_areas: