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

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

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

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

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

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

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

58
__author__ = 'Daniel Scheffler'
59
60


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

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

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

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

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

141
142
        self.mem_usage = {}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return self._MetaObj

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

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

294
295
        self._MetaObj = None

296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
    @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

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

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

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

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

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

        return self._coreg_info

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

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

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

406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
    @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

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

443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
    @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

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

        #     self.build_combined_masks_array()

560
561
562
        return self._masks

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

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

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

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

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

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

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

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

            self._ac_errors = errArr
        else:
            del self.ac_errors

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

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

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

        :return:
        """

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

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

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

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

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

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

728
        return self._accuracy_layers
729
730

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

751
752
753
754
    @accuracy_layers.deleter
    def accuracy_layers(self):
        self._accuracy_layers = None

755
756
757
758
759
760
761
762
763
764
765
    @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

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

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

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

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

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

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

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

813
814
        return self._ac_options

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

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

        out_dict = self.__dict__.copy()

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

        # 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

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

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

        out_dict = self.__dict__.copy()

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

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

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

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

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

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

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

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

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

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

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

911
        return GMS_obj
912

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

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

935
936
937
938
939
        ##################
        # merge logfiles #
        ##################

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

        allLogs_df = allLogs_df.sort_values(0)

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

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

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

964
965
966
967
        ##################
        # MERGE METADATA #
        ##################

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

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

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

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

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

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

1012
1013
1014
        ####################
        # MERGE ARRAY DATA #
        ####################
1015

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

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

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

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

1052
1053
            # validate output GeoArrays #
            #############################
1054

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

1065
1066
            # set output arrays #
            #####################
1067

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

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

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

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

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

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

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

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

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

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

1137
1138
1139
                if GMS_obj_merged.ac_errors:
                    GMS_obj_merged.ac_errors[mask_nonclear] = \
                        DEF_D.get_outFillZeroSaturated(GMS_obj_merged.ac_errors.dtype)[0]