gms_object.py 118 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
        return self._accuracy_layers
707
708

    @accuracy_layers.setter
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
709
    def accuracy_layers(self, *geoArr_initArgs):
710
        if geoArr_initArgs[0] is not None:
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
711
            acc_lay = GeoArray(*geoArr_initArgs)
712
713
714
715
            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
716
            # noinspection PyProtectedMember
717
718
719
720
721
722
723
724
725
726
727
728
            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

729
730
731
732
    @accuracy_layers.deleter
    def accuracy_layers(self):
        self._accuracy_layers = None

733
734
735
736
737
738
739
740
741
742
743
    @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

744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
    def generate_accuracy_layers(self):
        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)
            return None

        try:
            from ..algorithms.L2C_P import AccuracyCube
            self._accuracy_layers = AccuracyCube(self)  # don't use setter because it sets it to a GeoArray instance

        except ValueError as e:
            if str(e) == 'The given GMS_object contains no accuracy layers for combination.':
                if any([CFG.ac_estimate_accuracy, CFG.spathomo_estimate_accuracy, CFG.spechomo_estimate_accuracy]):
                    self.logger.warning('The given GMS_object contains no accuracy layers although computation '
                                        'of accurracy layers was enabled in job configuration.')
                else:
                    pass  # self._accuracy_layers keeps None
            else:
                raise

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

767
768
769
    @property
    def cloud_masking_algorithm(self):
        if not self._cloud_masking_algorithm:
770
            self._cloud_masking_algorithm = CFG.cloud_masking_algorithm[self.satellite]
771
772
        return self._cloud_masking_algorithm

773
774
    @property
    def ac_options(self):
775
        # type: () -> dict
776
        """
777
778
        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.
779
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
780

781
        if not self._ac_options:
782
            path_ac_options = CFG.path_custom_sicor_options or PG.get_path_ac_options(self.GMS_identifier)
783

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

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

801
802
803
804
805
806
                # 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
807
                # opt_dict['AC']['n_cores'] = CFG.CPUs if CFG.allow_subMultiprocessing else 1
808

809
                self._ac_options = opt_dict
810
811
812
            else:
                self.logger.warning('There is no options file available for atmospheric correction. '
                                    'Atmospheric correction must be skipped.')
813

814
815
        return self._ac_options

816
    def get_copied_dict_and_props(self, remove_privates=False):
817
        # type: (bool) -> dict
818
        """Returns a copy of the current object dictionary including the current values of all object properties."""
819
820

        # loggers must be closed
821
        self.close_loggers()
822
823
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
824
825
826
827
828

        out_dict = self.__dict__.copy()

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

        # 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

839
840
    def attributes2dict(self, remove_privates=False):
        # type: (bool) -> dict
841
        """Returns a copy of the current object dictionary including the current values of all object properties."""
842
843

        # loggers must be closed
844
        self.close_loggers()
845
846
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
847
848
849
850

        out_dict = self.__dict__.copy()

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

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

862
        self._loggers_disabled = False  # enables automatic recreation of loggers
863
864
        return out_dict

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

875
876
    @classmethod
    def from_disk(cls, tuple_GMS_subset):
877
878
879
880
        """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])
        """
881
        GMS_obj = cls()
882

883
        path_GMS_file = tuple_GMS_subset[0]
884
        GMSfileDict = INP_R.GMSfile2dict(path_GMS_file)
885
886

        # copy all attributes from GMS file (private attributes are not touched since they are not included in GMS file)
887
888
889
890
891
892
893

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

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

905
        GMS_obj.arr_shape, GMS_obj.arr_pos = tuple_GMS_subset[1]
906

907
        GMS_obj.arr = GMS_obj.pathGen.get_path_imagedata()
908
        # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters):
909
        GMS_obj.masks = GMS_obj.pathGen.get_path_maskdata()
910
        GMS_obj.accuracy_layers = GMS_obj.pathGen.get_path_accuracylayers()
911

912
        return GMS_obj
913

914
915
    @classmethod
    def from_sensor_subsystems(cls, list_GMS_objs):
916
        # type: (List[GMS_object]) -> GMS_object
917
918
919
920
921
        """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)
        """
922
        GMS_obj_merged = cls()
923

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

936
937
938
939
940
        ##################
        # merge logfiles #
        ##################

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

949
950
        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
951
952

        # set common metadata, needed for logfile
953
954
955
        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
956
957

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

961
        # log
962
963
        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]))
964

965
966
967
968
        ##################
        # MERGE METADATA #
        ##################

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

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

        # update layer-dependent metadata with respect to remaining input GMS objects
988
989
990
        GMS_obj_merged.MetaObj.LayerBandsAssignment = GMS_obj_merged.LayerBandsAssignment
        GMS_obj_merged.MetaObj.Subsystem = ''

991
992
        GMS_obj_merged.subsystem = ''
        del GMS_obj_merged.pathGen  # must be refreshed because subsystem is now ''
993
        GMS_obj_merged.close_loggers()  # must also be refreshed because it depends on pathGen
994

995
        for attrN in layerdependent_metadata:
996
997
998
999
1000
1001
            # 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)))
1002
                elif isinstance(attr_val, (dict, OrderedDict)):
1003
1004
1005
1006
1007
1008
                    attrDic_fullLBA.update(attr_val)
                else:
                    raise ValueError(attrN)

            # update the attribute in self.MetaObj
            if attrDic_fullLBA:
1009
                val2set = [attrDic_fullLBA[bN] for bN in GMS_obj_merged.LayerBandsAssignment] \
1010
                    if isinstance(getattr(list_GMS_objs[0].MetaObj, attrN), list) else attrDic_fullLBA
1011
                setattr(GMS_obj_merged.MetaObj, attrN, val2set)
1012

1013
1014
1015
        ####################
        # MERGE ARRAY DATA #
        ####################
1016

1017
1018
1019
1020
        # 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]))

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

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

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

1053
1054
            # validate output GeoArrays #
            #############################
1055

1056
1057
            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
1058
                                    for gA in geoArrs_same_extent[1:]])
1059
1060
                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:]])
1061
1062
1063
1064
1065
                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.')

1066
1067
            # set output arrays #
            #####################
1068

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

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

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

1097
1098
                elif attrname == 'mask_nodata':
                    # must not be merged -> self.arr is already merged, so just recalculate it (np.all)
1099
                    GMS_obj_merged.mask_nodata = GMS_obj_merged.calc_mask_nodata(overwrite=True)
1100

1101
1102
1103
                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]
1104
1105
                    if len(mask_clouds) > 1:
                        raise ValueError('Expected mask clouds in only one subsystem. Got %s.' % len(mask_clouds))
1106
                    GMS_obj_merged.mask_clouds = mask_clouds[0] if mask_clouds else None
1107

1108
1109
1110
                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]
1111
1112
1113
                    if len(mask_clouds_conf) > 1:
                        raise ValueError(
                            'Expected mask_clouds_conf in only one subsystem. Got %s.' % len(mask_clouds_conf))
1114
                    GMS_obj_merged.mask_clouds_confidence = mask_clouds_conf[0] if mask_clouds_conf else None
1115

1116
                elif attrname == 'masks':
1117
                    # self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
1118
                    GMS_obj_merged.masks = None
1119

1120
1121
1122
1123
1124
1125
1126
1127
        # 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]
1128
1129
1130
1131

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

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