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

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

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

28
29
30
31
try:
    from osgeo import gdalnumeric
except ImportError:
    import gdalnumeric
32

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

42
43
from ..misc.logging import GMS_logger as DatasetLogger
from ..misc import database_tools as DB_T
44
45
46
47
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
48
from ..options.config import GMS_config as CFG
49
50
51
from ..algorithms import geoprocessing as GEOP
from ..io import input_reader as INP_R
from ..io import output_writer as OUT_W
52
53
from ..misc import helper_functions as HLP_F
from ..misc import definition_dicts as DEF_D
54
from ..misc.locks import IOLock
55

56
57
58
if TYPE_CHECKING:
    from ..algorithms.L1C_P import L1C_object  # noqa F401  # flake8 issue

59
__author__ = 'Daniel Scheffler'
60
61


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

67
    def __init__(self, pathImage=''):
68
        # add private attributes
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
        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
88
        self._dict_LayerOptTherm = None
89
        self._cloud_masking_algorithm = None
90
        self._coreg_info = None
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
131
        # 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

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

142
143
        self.mem_usage = {}

144
        # set pathes
145
146
147
148
        self.path_cloud_class_obj = ''

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

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

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

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

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
163
164
        return self.__dict__

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

184
    def set_pathes(self):
185
186
187
188
189
190
        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()
191

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

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

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

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

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

        # save prior logs
239
        # if logger is None and self._logger is not None:
240
        #    self.log += self.logger.captured_stream
241
242
        self._logger = logger

243
244
245
246
247
248
249
250
251
252
253
    @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

254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
    @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

270
271
    @property
    def GMS_identifier(self):
272
273
        return GMS_identifier(self.image_type, self.satellite, self.sensor, self.subsystem, self.proc_level,
                              self.dataset_ID, self.logger)
274
275
276

    @property
    def MetaObj(self):
277
278
        # 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
279
280
281
282
283
284

        return self._MetaObj

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

    @MetaObj.deleter
    def MetaObj(self):
290
291
292
293
294
        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

295
296
        self._MetaObj = None

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
339
340
341
342
    @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

343
344
345
346
347
    @property
    def dict_LayerOptTherm(self):
        if self._dict_LayerOptTherm:
            return self._dict_LayerOptTherm
        elif self.LayerBandsAssignment:
348
            self._dict_LayerOptTherm = get_dict_LayerOptTherm(self.GMS_identifier, self.LayerBandsAssignment)
349
350
351
352
353
354
355
            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
356

357
358
359
360
        return True if self.image_type == 'RSD' and re.search('OLI', self.sensor, re.I) else False

    @property
    def coreg_needed(self):
361
        if self._coreg_needed is None:
362
            self._coreg_needed = not (self.dataset_ID == CFG.datasetid_spatial_ref)
363
        return self._coreg_needed
364
365
366
367
368

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

369
370
371
372
373
374
    @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},
375
                'original map info': self.MetaObj.map_info,
376
                'updated map info': None,
377
                'shift_reliability': None,
378
379
380
381
                '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
382
                'reference projection': self.MetaObj.projection,
383
                'reference extent': {'rows': None, 'cols': None},
384
385
                'reference grid': [list(CFG.spatial_ref_gridx),
                                   list(CFG.spatial_ref_gridy)],
386
387
388
389
390
391
392
393
394
                'success': False
            }

        return self._coreg_info

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

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

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

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
436
437
438
439
    @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

440
441
442
443
    @property
    def arr_meta(self):
        return self.MetaObj.to_odict()

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

554
555
    @property
    def masks(self):
556
        # if self.mask_nodata is not None and self.mask_clouds is not None and \
557
558
559
560
        #     self._masks is not None and self._masks.bands==1:

        #     self.build_combined_masks_array()

561
562
563
        return self._masks

    @masks.setter
564
    def masks(self, *geoArr_initArgs):
565
566
567
568
        """
        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
569

570
        if geoArr_initArgs[0] is not None:
571
            self._masks = GeoArray(*geoArr_initArgs)
572
            self._masks.nodata = 0
573
574
            self._masks.gt = self.arr.gt
            self._masks.prj = self.arr.prj
575
576
        else:
            del self.masks
577

578
579
580
581
    @masks.deleter
    def masks(self):
        self._masks = None

582
583
584
585
586
587
588
589
590
591
592
593
594
    @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
595
            # noinspection PyProtectedMember
596
597
598
599
            if cnfArr._nodata is None:
                cnfArr.nodata = DEF_D.get_outFillZeroSaturated(cnfArr.dtype)[0]
            cnfArr.gt = self.arr.gt
            cnfArr.prj = self.arr.prj
600
            cnfArr.bandnames = ['confidence']
601
602
603

            self._mask_clouds_confidence = cnfArr
        else:
604
            del self.mask_clouds_confidence
605
606
607
608
609
610
611
612
613
614
615
616

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

617
        return self._ac_errors
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632

    @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
633
            # noinspection PyProtectedMember
634
635
636
637
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
638
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
639
640
641
642
643
644
645
646
647

            self._ac_errors = errArr
        else:
            del self.ac_errors

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

648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    @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

666
667
668
669
670
671
672
    @property
    def spec_homo_errors(self):
        """Returns an instance of GeoArray containing error information calculated during spectral homogenization.

        :return:
        """

673
        return self._spec_homo_errors
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688

    @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
689
            # noinspection PyProtectedMember
690
691
692
693
            if errArr._nodata is None:
                errArr.nodata = DEF_D.get_outFillZeroSaturated(errArr.dtype)[0]
            errArr.gt = self.arr.gt
            errArr.prj = self.arr.prj
694
            errArr.bandnames = self.LBA2bandnames(self.LayerBandsAssignment) if errArr.ndim == 3 else ['median']
695
696
697
698
699
700
701
702
703
704
705

            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):
706
707
708
709
        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
710
                return None
711

712
713
714
            try:
                from ..algorithms.L2C_P import AccuracyCube
                self._accuracy_layers = AccuracyCube(self)
715
716
717

            except ValueError as e:
                if str(e) == 'The given GMS_object contains no accuracy layers for combination.':
Daniel Scheffler's avatar
Bugfix.    
Daniel Scheffler committed
718
719
720
721
722
723
724
725
                    if any([CFG.ac_estimate_accuracy, CFG.spathomo_estimate_accuracy, CFG.spechomo_estimate_accuracy]):
                        if self.proc_level == 'L2C':
                            self.logger.warning('The given GMS_object contains no accuracy layers although computation '
                                                'of accurracy layers was enabled in job configuration.')
                        else:
                            # Suppress the warning if accuracy_layers is just called by GMS_object.to_tiles() within
                            # L2A or L2B processing.
                            pass
726
727
728
729
730
                    else:
                        pass  # self._accuracy_layers keeps None
                else:
                    raise

731
732
            except Exception as e:
                raise RuntimeError('Failed to generate AccuracyCube!', e)
733

734
        return self._accuracy_layers
735
736

    @accuracy_layers.setter
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
737
    def accuracy_layers(self, *geoArr_initArgs):
738
        if geoArr_initArgs[0] is not None:
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
739
            acc_lay = GeoArray(*geoArr_initArgs)
740
741
742
743
            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
744
            # noinspection PyProtectedMember
745
746
747
748
749
750
751
752
753
754
755
756
            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

757
758
759
760
    @accuracy_layers.deleter
    def accuracy_layers(self):
        self._accuracy_layers = None

761
762
763
764
765
766
767
768
769
770
771
    @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

772
773
774
    @property
    def cloud_masking_algorithm(self):
        if not self._cloud_masking_algorithm:
775
            self._cloud_masking_algorithm = CFG.cloud_masking_algorithm[self.satellite]
776
777
        return self._cloud_masking_algorithm

778
779
    @property
    def ac_options(self):
780
        # type: () -> dict
781
        """
782
783
        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.
784
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
785

786
        if not self._ac_options:
787
            path_ac_options = CFG.path_custom_sicor_options or PG.get_path_ac_options(self.GMS_identifier)
788

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

Daniel Scheffler's avatar
Daniel Scheffler committed
793
                # update some file paths depending on the current environment
794
795
                opt_dict['DEM']['fn'] = CFG.path_dem_proc_srtm_90m
                opt_dict['ECMWF']['path_db'] = CFG.path_ECMWF_db
796
797
798
                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
799
                opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]')
800
                if 'uncertainties' in opt_dict:
801
802
803
804
                    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
805

806
807
808
809
810
811
                # 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
812
                # opt_dict['AC']['n_cores'] = CFG.CPUs if CFG.allow_subMultiprocessing else 1
813

814
                self._ac_options = opt_dict
815
816
817
            else:
                self.logger.warning('There is no options file available for atmospheric correction. '
                                    'Atmospheric correction must be skipped.')
818

819
820
        return self._ac_options

821
    def get_copied_dict_and_props(self, remove_privates=False):
822
        # type: (bool) -> dict
823
        """Returns a copy of the current object dictionary including the current values of all object properties."""
824
825

        # loggers must be closed
826
        self.close_loggers()
827
828
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
829
830
831
832
833

        out_dict = self.__dict__.copy()

        # add properties
        property_names = [p for p in dir(self.__class__) if isinstance(getattr(self.__class__, p), property)]
834
        [out_dict.update({propK: copy.copy(getattr(self, propK))}) for propK in property_names]
835
836
837
838
839
840
841
842
843

        # 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

844
845
    def attributes2dict(self, remove_privates=False):
        # type: (bool) -> dict
846
        """Returns a copy of the current object dictionary including the current values of all object properties."""
847
848

        # loggers must be closed
849
        self.close_loggers()
850
851
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
852
853
854
855

        out_dict = self.__dict__.copy()

        # add some selected property values
856
        for i in ['GMS_identifier', 'LayerBandsAssignment', 'coreg_needed', 'coreg_info', 'resamp_needed',
857
858
859
860
861
                  'dict_LayerOptTherm', 'georef', 'MetaObj']:
            if i == 'MetaObj':
                out_dict['meta_odict'] = self.MetaObj.to_odict()
            else:
                out_dict[i] = getattr(self, i)
862
863
864
865
866

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

867
        self._loggers_disabled = False  # enables automatic recreation of loggers
868
869
        return out_dict

870
    def _data_downloader(self, sensor, entity_ID):
871
872
873
874
        self.logger.info('Data downloader started.')
        success = False
        " > download source code for Landsat here < "
        if not success:
875
876
            self.logger.critical(
                "Download for %s dataset '%s' failed. No further processing possible." % (sensor, entity_ID))
877
            raise RuntimeError('Archive download failed.')
878
879
        return success

880
881
    @classmethod
    def from_disk(cls, tuple_GMS_subset):
882
883
884
885
        """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])
        """
886
        GMS_obj = cls()
887

888
        path_GMS_file = tuple_GMS_subset[0]
889
        GMSfileDict = INP_R.GMSfile2dict(path_GMS_file)
890
891

        # copy all attributes from GMS file (private attributes are not touched since they are not included in GMS file)
892
893
894
895
896
897
898

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

899
900
        for key, value in GMSfileDict.items():
            if key in ['GMS_identifier', 'georef', 'dict_LayerOptTherm']:
901
                continue  # properties that should better be created on the fly
902
903
904
905
906
907
908
            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)
909

910
        GMS_obj.arr_shape, GMS_obj.arr_pos = tuple_GMS_subset[1]
911

912
        GMS_obj.arr = GMS_obj.pathGen.get_path_imagedata()
913
        # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters):
914
        GMS_obj.masks = GMS_obj.pathGen.get_path_maskdata()
915
        GMS_obj.accuracy_layers = GMS_obj.pathGen.get_path_accuracylayers()
916

917
        return GMS_obj
918

919
920
    @classmethod
    def from_sensor_subsystems(cls, list_GMS_objs):
921
        # type: (List[GMS_object]) -> GMS_object
922
923
924
925
926
        """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)
        """
927
        GMS_obj_merged = cls()
928

929
        # assertions
930
931
        assert len(list_GMS_objs) > 1, "'GMS_object.from_sensor_subsystems()' expects multiple input GMS objects. " \
                                       "Got %d." % len(list_GMS_objs)
932
        assert all([is_coord_grid_equal(list_GMS_objs[0].arr.gt, *obj.arr.xygrid_specs) for obj in list_GMS_objs[1:]]),\
933
934
935
            "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, \
936
937
938
            "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))), \
939
            "The input 'list_GMS_objs' contains duplicates: %s" % subsystems
940

941
942
943
944
945
        ##################
        # merge logfiles #
        ##################

        # read all logs into DataFrame, sort it by the first column
946
        [GMS_obj.close_loggers() for GMS_obj in list_GMS_objs]  # close the loggers of the input objects
947
948
949
        paths_inLogs = [GMS_obj.pathGen.get_path_logfile() for GMS_obj in list_GMS_objs]
        allLogs_df = DataFrame()
        for log in paths_inLogs:
950
            df = read_csv(log, sep='\n', delimiter=None, header=None,  # no delimiter needed
951
                          engine='python')  # engine suppresses a pandas warning
952
            allLogs_df = allLogs_df.append(df)
953

954
955
        allLogs_df = allLogs_df.sort_values(0)  # sorting uses timestamps that appear on first position in logs
        allLogs_df = allLogs_df.drop_duplicates()  # otherwise, e.g., logs from AC would appear 3 times for S2A
956
957

        # set common metadata, needed for logfile
958
959
960
        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
961
962

        # write the merged logfile and flush previous logger
963
        np.savetxt(GMS_obj_merged.path_logfile, np.array(allLogs_df), delimiter=None, fmt="%s")
964
        GMS_obj_merged.close_loggers()
965

966
        # log
967
968
        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]))
969

970
971
972
973
        ##################
        # MERGE METADATA #
        ##################

974
975
976
        # 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']:
977
                continue  # properties that should better be created on the fly
978
979
            elif key in ['baseN', 'path_logfile', 'scene_ID', 'subsystem']:
                continue  # either previously set with common values or not needed for merged GMS_object
980
            try:
981
                setattr(GMS_obj_merged, key, value)
982
983
            except Exception:
                raise AttributeError("Can't set attribute %s." % key)
984

985
        # update LayerBandsAssignment and get full list of output bandnames
986
        from .metadata import get_LayerBandsAssignment
987
988
        # use identifier of first input GMS object for getting LBA (respects current proc_level):
        gms_idf = list_GMS_objs[0].GMS_identifier
989
990
        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]
991
992

        # update layer-dependent metadata with respect to remaining input GMS objects
993
994
995
        GMS_obj_merged.MetaObj.LayerBandsAssignment = GMS_obj_merged.LayerBandsAssignment
        GMS_obj_merged.MetaObj.Subsystem = ''

996
997
        GMS_obj_merged.subsystem = ''
        del GMS_obj_merged.pathGen  # must be refreshed because subsystem is now ''
998
        GMS_obj_merged.close_loggers()  # must also be refreshed because it depends on pathGen
999

1000
        for attrN in layerdependent_metadata:
1001
1002
1003
1004
1005
1006
            # 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)))
1007
                elif isinstance(attr_val, (dict, OrderedDict)):
1008
1009
1010
1011
1012
1013
                    attrDic_fullLBA.update(attr_val)
                else:
                    raise ValueError(attrN)

            # update the attribute in self.MetaObj
            if attrDic_fullLBA:
1014
                val2set = [attrDic_fullLBA[bN] for bN in GMS_obj_merged.LayerBandsAssignment] \
1015
                    if isinstance(getattr(list_GMS_objs[0].MetaObj, attrN), list) else attrDic_fullLBA
1016
                setattr(GMS_obj_merged.MetaObj, attrN, val2set)
1017

1018
1019
1020
        ####################
        # MERGE ARRAY DATA #
        ####################
1021

1022
1023
1024
1025
        # 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]))

1026
        # overwrite array data with merged arrays, clipped to common_extent and reordered according to FullLayerBandsAss
1027
1028
1029
        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
1030
            all_arrays = [getattr(GMS_obj, '_%s' % attrname) for GMS_obj in list_GMS_objs]
1031
1032
1033
1034
1035
1036
1037
1038
1039

            # 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:
1040
                        # FIXME mask_clouds_confidence is no GeoArray until here
1041
                        # FIXME -> has no nodata value -> calculation throughs warning
1042
1043
                        geoArr_same_extent = \
                            GeoArray(*geoArr.get_mapPos(
1044
1045
1046
1047
                                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)
1048
1049
                        geoArrs_same_extent.append(geoArr_same_extent)
                    else:
1050
1051
                        # e.g. in case of cloud mask that is only appended to the GMS object with the same
                        # spatial resolution)
1052
1053
1054
1055
1056
1057
                        geoArrs_same_extent.append(None)

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

1058
1059
            # validate output GeoArrays #
            #############################
1060

1061
1062
            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
1063
                                    for gA in geoArrs_same_extent[1:]])
1064
1065
                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:]])
1066
1067
1068
1069
1070
                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.')

1071
1072
            # set output arrays #
            #####################
1073

1074
1075
            # 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]:
1076
1077
                # check that each desired band name for the current attribute is provided by one of the input
                # GMS objects
1078
1079
                available_bandNs = list(chain.from_iterable([list(gA.bandnames) for gA in geoArrs_same_extent]))
                for bN in bandnames:
1080
                    if bN not in available_bandNs:
1081
                        raise ValueError("The given input GMS objects (subsystems) do not provide a bandname '%s' for "
1082
1083
                                         "the attribute '%s'. Available band names amongst all input GMS objects are: "
                                         "%s" % (bN, attrname, str(available_bandNs)))
1084
1085

                # merge arrays
1086
1087
                def get_band(bandN):
                    return [gA[bandN] for gA in geoArrs_same_extent if gA and bandN in gA.bandnames][0]
1088
1089
                full_geoArr = GeoArray(np.dstack((get_band(bandN) for bandN in bandnames)),
                                       geoArrs_same_extent[0].gt, geoArrs_same_extent[0].prj,
1090
1091
                                       bandnames=bandnames,
                                       nodata=geoArrs_same_extent[0].nodata)
1092
                setattr(GMS_obj_merged, attrname, full_geoArr)
1093

1094
            # handle the remaining arrays
1095
1096
1097
            else:
                # masks, dem, mask_nodata, mask_clouds, mask_clouds_confidence
                if attrname == 'dem':
1098
1099
                    # use the DEM of the first input object
                    # (if the grid is the same, the DEMs should be the same anyway)
1100
                    GMS_obj_merged.dem = geoArrs_same_extent[0]
1101

1102
1103
                elif attrname == 'mask_nodata':
                    # must not be merged -> self.arr is already merged, so just recalculate it (np.all)
1104
                    GMS_obj_merged.mask_nodata = GMS_obj_merged.calc_mask_nodata(overwrite=True)
1105

1106
1107
1108
                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]
1109
1110
                    if len(mask_clouds) > 1:
                        raise ValueError('Expected mask clouds in only one subsystem. Got %s.' % len(mask_clouds))
1111
                    GMS_obj_merged.mask_clouds = mask_clouds[0] if mask_clouds else None
1112

1113
1114
1115
                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]
1116
1117
1118
                    if len(mask_clouds_conf) > 1:
                        raise ValueError(
                            'Expected mask_clouds_conf in only one subsystem. Got %s.' % len(mask_clouds_conf))
1119
                    GMS_obj_merged.mask_clouds_confidence = mask_clouds_conf[0] if mask_clouds_conf else None
1120

1121
                elif attrname == 'masks':
1122
                    # self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
1123
                    GMS_obj_merged.masks = None
1124

1125
1126
1127
1128
1129
1130
1131
1132
        # 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]
1133
1134
1135
1136

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