gms_object.py 71.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# -*- coding: utf-8 -*-

import collections
import copy
import datetime
import glob
import json
import os
import re
import shutil
import sys
import warnings
13
import logging
14
from itertools import chain
15
16
17

import numpy as np
import spectral
18
from spectral.io import envi
19
from numba import jit
20
from pandas import DataFrame, read_csv
21

22
23
24
25
try:
    from osgeo import gdalnumeric
except ImportError:
    import gdalnumeric
26

27
from geoarray import GeoArray
28
from py_tools_ds.geo.coord_grid import is_coord_grid_equal
29
from py_tools_ds.geo.projection import EPSG2WKT
30
31
32
from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX
33
from sicor.options import get_options as get_ac_options
34

35
36
from ..misc.logging import GMS_logger as DatasetLogger
from ..misc.mgrs_tile import MGRS_tile
37
from ..model.metadata import METADATA, get_dict_LayerOptTherm, metaDict_to_metaODict
38
39
40
41
from ..model.dataset import Dataset
from ..misc import path_generator as PG
from ..misc import database_tools as DB_T
from ..config import GMS_config as CFG
42
43
44
from ..algorithms import geoprocessing as GEOP
from ..io import input_reader as INP_R
from ..io import output_writer as OUT_W
45
46
from ..misc import helper_functions as HLP_F
from ..misc import definition_dicts as DEF_D
47

48
__author__ = 'Daniel Scheffler'
49
50


51
class GMS_object(Dataset):
52
53
54
55
    def __init__(self, pathImage=''):
        # get all attributes of base class "Dataset"
        super(GMS_object, self).__init__()

56
        # add private attributes
57
        self._dict_LayerOptTherm = None
58
59
        self._cloud_masking_algorithm = None
        self._meta_odict = None
60
61

        self.job_ID = CFG.job.ID
62
63
64
65
66
        # FIXME not needed anymore?:
        # self.dataset_ID = int(DB_T.get_info_from_postgreSQLdb(CFG.job.conn_database, 'scenes', ['datasetid'],
        #                                {'id': self.scene_ID})[0][0]) if self.scene_ID !=-9999 else -9999
        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
67
        self.MGRS_info = None
68
69

        # set pathes
70
71
72
73
        self.path_cloud_class_obj = ''

        # handle initialization arguments
        if pathImage:
74
75
            # run the setter for 'arr' of the base class 'Dataset' which creates an Instance of GeoArray
            self.arr = pathImage
76
77
78
79

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

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
80
        self.close_loggers()
81
        del self.pathGen  # path generator can only be used for the current processing level
82
83

        # delete arrays if their in-mem size is to big to be pickled
84
        # => (avoids MaybeEncodingError: Error sending result: '[<gms_preprocessing.algorithms.L2C_P.L2C_object
85
        #    object at 0x7fc44f6399e8>]'. Reason: 'error("'i' format requires -2147483648 <= number <= 2147483647",)')
86
87
88
        if self.proc_level == 'L2C' and CFG.job.exec_mode == 'Flink':
            # FIXME check by bandname
            if self.mask_nodata is not None and self.masks.bands > 1 and self.mask_clouds is not None:
89
                del self.masks
90

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
91
92
        return self.__dict__

93
    def set_pathes(self):
94
95
96
97
98
99
        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()
100

101
        if CFG.job.exec_mode == 'Python':
Daniel Scheffler's avatar
Daniel Scheffler committed
102
            self.path_InFilePreprocessor = os.path.join(self.ExtractedFolder, '%s%s_DN.bsq'
103
104
105
106
                                                        % (self.entity_ID,
                                                           ('_%s' % self.subsystem if self.subsystem else '')))
        else:  # Flink
            self.path_InFilePreprocessor = None  # None: keeps all produced data in memory (numpy array attributes)
107
108
109
110
111
112

        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"
113
                             "has been found at %s." % (self.sensor, self.entity_ID, self.path_archive))
114
            self.logger.info('Trying to download the dataset...')
115
            self.path_archive_valid = self._data_downloader(self.sensor, self.entity_ID)
116
117
118
        else:
            self.path_archive_valid = True

119
        if CFG.job.exec_mode == 'Python' and self.ExtractedFolder and not os.path.isdir(self.ExtractedFolder):
120
121
122
123
124
125
126
            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)
127
128
129
        if CFG.job.exec_mode == 'Python' and self.ExtractedFolder:
            assert os.path.exists(self.path_archive), \
                'Invalid path for temporary files. Directory %s does not exist.' % self.ExtractedFolder
130
131
132
133
134
135
136
137

    @property
    def logger(self):
        if self._loggers_disabled:
            return None
        if self._logger and self._logger.handlers[:]:
            return self._logger
        else:
138
139
            self._logger = DatasetLogger('log__' + self.baseN, fmt_suffix=self.scene_ID, path_logfile=self.path_logfile,
                                         log_level=CFG.job.log_level, append=True)
140
141
142
143
            return self._logger

    @logger.setter
    def logger(self, logger):
144
        assert isinstance(logger, logging.Logger) or logger in ['not set', None], \
145
            "GMS_obj.logger can not be set to %s." % logger
146
147

        # save prior logs
148
        # if logger is None and self._logger is not None:
149
        #    self.log += self.logger.captured_stream
150
151
152
153
        self._logger = logger

    @property
    def GMS_identifier(self):
154
155
        return collections.OrderedDict(zip(
            ['image_type', 'Satellite', 'Sensor', 'Subsystem', 'proc_level', 'dataset_ID', 'logger'],
156
157
            [self.image_type, self.satellite, self.sensor, self.subsystem, self.proc_level, self.dataset_ID,
             self.logger]))
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

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

        return self._MetaObj

    @MetaObj.setter
    def MetaObj(self, MetaObj):
        assert isinstance(MetaObj, METADATA), "'MetaObj' can only be set to an instance of METADATA class. " \
174
                                              "Got %s." % type(MetaObj)
175
176
177
        self._MetaObj = MetaObj

        # update meta_odict
178
        del self.meta_odict  # it is recreated if getter is used the next time
179
180
181
182
183
184
185
186

    @MetaObj.deleter
    def MetaObj(self):
        self._MetaObj = None

    @property
    def meta_odict(self):
        if self._MetaObj:
187
            # if there is already a MetaObj -> create new meta_odict from it (ensures synchronization!)
188
189
190
191
            self._meta_odict = self._MetaObj.to_odict()
            del self.MetaObj
        elif not self._meta_odict:
            # if there is no MetaObj and no meta_odict -> use MetaObj getter to read metadata from disk
192
            pass  # reading from disk should use L1A_P.L1A_object.import_metadata -> so just return None
193
194
195
196
197
198
            self._meta_odict = None

        return self._meta_odict

    @meta_odict.setter
    def meta_odict(self, odict):
199
200
        assert isinstance(odict, (collections.OrderedDict, dict)), "'meta_odict' can only be set to an instance of " \
                                                                   "collections.OrderedDict. Got %s." % type(odict)
201
202
203
        self._meta_odict = odict

        # update MetaObj
204
        del self.MetaObj  # it is recreated if getter is used the next time
205
206
207
208
209

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

210
211
212
213
214
    @property
    def dict_LayerOptTherm(self):
        if self._dict_LayerOptTherm:
            return self._dict_LayerOptTherm
        elif self.LayerBandsAssignment:
215
            self._dict_LayerOptTherm = get_dict_LayerOptTherm(self.identifier, self.LayerBandsAssignment)
216
217
218
219
220
221
222
            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
223

224
225
226
227
        return True if self.image_type == 'RSD' and re.search('OLI', self.sensor, re.I) else False

    @property
    def coreg_needed(self):
228
229
230
        if self._coreg_needed is None:
            self._coreg_needed = not (self.dataset_ID == CFG.usecase.datasetid_spatial_ref)
        return self._coreg_needed
231
232
233
234
235

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

236
237
238
239
    @property
    def resamp_needed(self):
        if self._resamp_needed is None:
            gt = mapinfo2geotransform(self.meta_odict['map info'])
240
241
            self._resamp_needed = not is_coord_grid_equal(gt, CFG.usecase.spatial_ref_gridx,
                                                          CFG.usecase.spatial_ref_gridy)
242
243
244
245
246
247
248
249
        return self._resamp_needed

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

    @property
    def masks(self):
250
        # if self.mask_nodata is not None and self.mask_clouds is not None and \
251
252
253
254
        #     self._masks is not None and self._masks.bands==1:

        #     self.build_combined_masks_array()

255
256
257
        return self._masks

    @masks.setter
258
    def masks(self, *geoArr_initArgs):
259
260
261
262
        """
        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
263

264
        if geoArr_initArgs[0] is not None:
265
            self._masks = GeoArray(*geoArr_initArgs)
266
            self._masks.nodata = 0
267
268
            self._masks.gt = self.arr.gt
            self._masks.prj = self.arr.prj
269
270
        else:
            del self.masks
271

272
273
274
275
    @masks.deleter
    def masks(self):
        self._masks = None

276
277
278
279
280
281
    @property
    def cloud_masking_algorithm(self):
        if not self._cloud_masking_algorithm:
            self._cloud_masking_algorithm = CFG.job.cloud_masking_algorithm[self.satellite]
        return self._cloud_masking_algorithm

282
283
284
    @property
    def ac_options(self):
        """
285
286
        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.
287
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
288

289
        if not self._ac_options:
290
291
            path_ac_options = PG.get_path_ac_options(self.GMS_identifier)

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

Daniel Scheffler's avatar
Daniel Scheffler committed
296
                # update some file paths depending on the current environment
297
298
299
300
301
                opt_dict['DEM']['fn'] = CFG.job.path_dem_proc_srtm_90m
                opt_dict['ECMWF']['path_db'] = CFG.job.path_ECMWF_db
                for key in opt_dict['RTFO']:
                    if 'atm_tables_fn' in opt_dict['RTFO'][key]:
                        opt_dict['RTFO'][key]['atm_tables_fn'] = PG.get_path_ac_table(key)
302
303
                opt_dict['S2Image'][
                    'S2_MSI_granule_path'] = None  # only a placeholder -> will always be None for GMS usage
304
305
306
                if CFG.job.cloud_masking_algorithm[self.satellite] == 'SICOR':
                    opt_dict['cld_mask']['persistence_file'] = PG.get_path_cloud_class_obj(self.GMS_identifier)
                    opt_dict['cld_mask']['novelty_detector'] = None  # FIXME update this after switching to SICOR
307
                opt_dict['output'] = []  # outputs are not needed for GMS -> so
308
                opt_dict['report']['report_path'] = os.path.join(self.pathGen.get_path_procdata(), '[TYPE]')
309
310
                if 'uncertainties' in opt_dict:
                    opt_dict['uncertainties']['snr_model'] = PG.get_path_snr_model(self.GMS_identifier)
311

Daniel Scheffler's avatar
Daniel Scheffler committed
312
                # opt_dict['AC']['n_cores'] = CFG.job.CPUs if CFG.job.allow_subMultiprocessing else 1
313

314
                self._ac_options = opt_dict
315
316
317
            else:
                self.logger.warning('There is no options file available for atmospheric correction. '
                                    'Atmospheric correction must be skipped.')
318

319
320
        return self._ac_options

321
    def get_copied_dict_and_props(self, remove_privates=False):
322
        # type: (bool) -> dict
323
        """Returns a copy of the current object dictionary including the current values of all object properties."""
324
325
326

        # loggers must be closed
        self.close_GMS_loggers()
327
328
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
329
330
331
332
333

        out_dict = self.__dict__.copy()

        # add properties
        property_names = [p for p in dir(self.__class__) if isinstance(getattr(self.__class__, p), property)]
334
        [out_dict.update({propK: copy.copy(getattr(self, propK))}) for propK in property_names]
335
336
337
338
339
340
341
342
343

        # 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

344
345
    def attributes2dict(self, remove_privates=False):
        # type: (bool) -> dict
346
        """Returns a copy of the current object dictionary including the current values of all object properties."""
347
348
349

        # loggers must be closed
        self.close_GMS_loggers()
350
351
        # this disables automatic recreation of loggers (otherwise loggers are created by using getattr()):
        self._loggers_disabled = True
352
353
354
355

        out_dict = self.__dict__.copy()

        # add some selected property values
356
        for i in ['GMS_identifier', 'LayerBandsAssignment', 'coreg_needed', 'resamp_needed', 'dict_LayerOptTherm',
357
                  'georef', 'meta_odict']:
358
            out_dict[i] = getattr(self, i)
359
360
361
362
363

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

364
        self._loggers_disabled = False  # enables automatic recreation of loggers
365
366
        return out_dict

367
    def _data_downloader(self, sensor, entity_ID):
368
369
370
371
        self.logger.info('Data downloader started.')
        success = False
        " > download source code for Landsat here < "
        if not success:
372
373
            self.logger.critical(
                "Download for %s dataset '%s' failed. No further processing possible." % (sensor, entity_ID))
374
            raise RuntimeError('Archive download failed.')
375
376
        return success

377
378
379
380
381
    def from_disk(self, tuple_GMS_subset):
        """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])
        """
382

383
        path_GMS_file = tuple_GMS_subset[0]
384
        GMSfileDict = INP_R.GMSfile2dict(path_GMS_file)
385
386

        # copy all attributes from GMS file (private attributes are not touched since they are not included in GMS file)
387
        self.meta_odict = GMSfileDict['meta_odict']  # set that first in order to make some getters and setters work
388
389
        for key, value in GMSfileDict.items():
            if key in ['GMS_identifier', 'georef', 'dict_LayerOptTherm']:
390
                continue  # properties that should better be created on the fly
391
392
            try:
                setattr(self, key, value)
393
394
            except Exception:
                raise AttributeError("Can't set attribute %s." % key)
395

396
        self.acq_datetime = datetime.datetime.strptime(self.acq_datetime, '%Y-%m-%d %H:%M:%S.%f%z')
397
398
        self.arr_shape, self.arr_pos = tuple_GMS_subset[1]

399
400
401
        self.arr = self.pathGen.get_path_imagedata()
        # self.mask_nodata and self.mask_clouds are auto-synchronized via self.masks (see their getters):
        self.masks = self.pathGen.get_path_maskdata()
402

403
404
        return copy.copy(self)

405
406
407
408
409
410
411
412
    def from_sensor_subsystems(self, list_GMS_objs):
        # type: (list) -> GMS_object
        """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)
        """

413
        # assertions
414
415
        assert len(list_GMS_objs) > 1, "'GMS_object.from_sensor_subsystems()' expects multiple input GMS objects. " \
                                       "Got %d." % len(list_GMS_objs)
416
        assert all([is_coord_grid_equal(list_GMS_objs[0].arr.gt, *obj.arr.xygrid_specs) for obj in list_GMS_objs[1:]]),\
417
418
419
            "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, \
420
421
422
            "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))), \
423
            "The input 'list_GMS_objs' contains duplicates: %s" % subsystems
424

425
        # log
426
        list_GMS_objs[0].logger.info('Merging the subsystems %s to a single GMS object...'
427
                                     % ', '.join([GMS_obj.subsystem for GMS_obj in list_GMS_objs]))
428
429

        # find the common extent. NOTE: boundsMap is expected in the order [xmin,xmax,ymin,ymax]
430
431
        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]))
432
433
434
435
436

        # MERGE METADATA
        # 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']:
437
                continue  # properties that should better be created on the fly
438
439
            try:
                setattr(self, key, value)
440
441
            except Exception:
                raise AttributeError("Can't set attribute %s." % key)
442

443
        # update LayerBandsAssignment and get full list of output bandnames
444
        from .metadata import get_LayerBandsAssignment
445
446
        # use identifier of first input GMS object for getting LBA (respects current proc_level):
        gms_idf = list_GMS_objs[0].GMS_identifier
447
        self.LayerBandsAssignment = get_LayerBandsAssignment(gms_idf, return_fullLBA=True)
448
        bandnames = ['B%s' % i if len(i) == 2 else 'B0%s' % i for i in self.LayerBandsAssignment]
449
450
451

        # update layer-dependent metadata with respect to remaining input GMS objects
        self.meta_odict.update({
452
453
454
455
            'band names': [('Band %s' % i) for i in self.LayerBandsAssignment],
            'LayerBandsAssignment': self.LayerBandsAssignment,
            'Subsystem': '',
            'PhysUnit': self.meta_odict['PhysUnit'],  # TODO can contain val_optical / val_thermal
456
457
        })
        self.subsystem = ''
458
459
        del self.pathGen  # must be refreshed because subsystem is now ''
        self.close_GMS_loggers()  # must also be refreshed because it depends on pathGen
460

461
462
        for attrN in ['SolIrradiance', 'CWL', 'FWHM', 'Offsets', 'OffsetsRef', 'Gains', 'GainsRef',
                      'ThermalConstK1', 'ThermalConstK2', 'ViewingAngle_arrProv', 'IncidenceAngle_arrProv']:
463
464
465
466
467
468
469

            # 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)))
470
                elif isinstance(attr_val, (dict, collections.OrderedDict)):
471
472
473
474
475
476
477
                    attrDic_fullLBA.update(attr_val)
                else:
                    raise ValueError(attrN)

            # update the attribute in self.MetaObj
            if attrDic_fullLBA:
                val2set = [attrDic_fullLBA[bN] for bN in self.LayerBandsAssignment] \
478
                    if isinstance(getattr(list_GMS_objs[0].MetaObj, attrN), list) else attrDic_fullLBA
479
480
481
                setattr(self.MetaObj, attrN, val2set)

        # merge logfiles (read all logs into DataFrame, sort it by the first column and write to new logfile
482
        [GMS_obj.close_GMS_loggers() for GMS_obj in list_GMS_objs]  # close the loggers of the input objects
483
        paths_inLogs = [GMS_obj.pathGen.get_path_logfile() for GMS_obj in list_GMS_objs]
484
        allLogs_df = DataFrame()
485
        for log in paths_inLogs:
486
487
488
489
            df = read_csv(log, sep='\n', delimiter=':   ', header=None,
                          engine='python')  # engine suppresses a pandas warning
            allLogs_df = allLogs_df.append(
                df)  # FIXME this will log e.g. atm. corr 3 times for S2A -> use captured streams instead?
490
491
492
493
494

        allLogs_df = allLogs_df.sort_values(0)
        np.savetxt(self.pathGen.get_path_logfile(), np.array(allLogs_df), delimiter=':   ', fmt="%s")

        # MERGE ARRAY DATA
495
        # overwrite array data with merged arrays, clipped to common_extent and reordered according to FullLayerBandsAss
496
497
498
        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
499
            all_arrays = [getattr(GMS_obj, '_%s' % attrname) for GMS_obj in list_GMS_objs]
500
501
502
503
504
505
506
507
508

            # 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:
509
510
                        # FIXME mask_clouds_confidence is until here no GeoArray
                        # FIXME -> has no nodata value -> calculation throughs warning
511
512
                        geoArr_same_extent = \
                            GeoArray(*geoArr.get_mapPos(
513
514
515
516
                                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)
517
518
                        geoArrs_same_extent.append(geoArr_same_extent)
                    else:
519
520
                        # e.g. in case of cloud mask that is only appended to the GMS object with the same
                        # spatial resolution)
521
522
523
524
525
526
527
                        geoArrs_same_extent.append(None)

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

            # validate output GeoArrays
528
529
            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
530
                                    for gA in geoArrs_same_extent[1:]])
531
532
                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:]])
533
534
535
536
537
538
539
540
541
                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.')

            # set output arrays
            if attrname in ['arr', 'ac_errors'] and list(set(geoArrs_same_extent)) != [None]:
                # the bands of these arrays have to be reordered according to FullLayerBandsAssignment

542
543
                # check that each desired band name for the current attribute is provided by one of the input
                # GMS objects
544
545
                available_bandNs = list(chain.from_iterable([list(gA.bandnames) for gA in geoArrs_same_extent]))
                for bN in bandnames:
546
                    if bN not in available_bandNs:
547
                        raise ValueError("The given input GMS objects (subsystems) do not provide a bandname '%s' for "
548
549
                                         "the attribute '%s'. Available band names amongst all input GMS objects are: "
                                         "%s" % (bN, attrname, str(available_bandNs)))
550
551

                # merge arrays
552
                def get_band(bandN): return [gA[bandN] for gA in geoArrs_same_extent if gA and bandN in gA.bandnames][0]
553
554
                full_geoArr = GeoArray(np.dstack((get_band(bandN) for bandN in bandnames)),
                                       geoArrs_same_extent[0].gt, geoArrs_same_extent[0].prj,
555
556
                                       bandnames=bandnames,
                                       nodata=geoArrs_same_extent[0].nodata)
557
558
559
560
561
                setattr(self, attrname, full_geoArr)

            else:
                # masks, dem, mask_nodata, mask_clouds, mask_clouds_confidence
                if attrname == 'dem':
562
563
                    # use the DEM of the first input object
                    # (if the grid is the same, the DEMs should be the same anyway)
564
565
566
567
568
569
570
                    self.dem = geoArrs_same_extent[0]
                elif attrname == 'mask_nodata':
                    # must not be merged -> self.arr is already merged, so just recalculate it (np.all)
                    self.mask_nodata = self.calc_mask_nodata(overwrite=True)
                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]
571
572
                    if len(mask_clouds) > 1:
                        raise ValueError('Expected mask clouds in only one subsystem. Got %s.' % len(mask_clouds))
573
574
575
576
                    self.mask_clouds = mask_clouds[0] if mask_clouds else None
                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]
577
578
579
                    if len(mask_clouds_conf) > 1:
                        raise ValueError(
                            'Expected mask_clouds_conf in only one subsystem. Got %s.' % len(mask_clouds_conf))
580
                    self.mask_clouds_confidence = mask_clouds_conf[0] if mask_clouds_conf else None
581
                elif attrname == 'masks':
582
583
584
585
586
587
588
589
                    # self.mask_nodata and self.mask_clouds will already be set here -> so just recreate it from there
                    self.masks = None

        # recreate self.masks
        self.build_combined_masks_array()

        # update array-dependent metadata
        self.meta_odict.update({
590
591
            'samples': self.arr.cols, 'lines': self.arr.rows, 'bands': self.arr.bands,
            'map info': geotransform2mapinfo(self.arr.gt, self.arr.prj), 'coordinate system string': self.arr.prj, })
592
593

        # set shape of full array
594
        self.shape_fullArr = self.arr.shape
595
596
597

        return copy.copy(self)

598
599
600
601
602
603
    def from_tiles(self, list_GMS_tiles):
        # type: (list) -> self
        """Merge separate GMS objects with different spatial coverage but belonging to the same scene-ID to ONE GMS object.

        :param list_GMS_tiles: <list> of GMS objects that have been created by cut_GMS_obj_into_blocks()
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
604

605
606
607
608
609
        if 'IMapUnorderedIterator' in str(type(list_GMS_tiles)):
            list_GMS_tiles = list(list_GMS_tiles)

        # copy all attributes except of array attributes
        tile1 = list_GMS_tiles[0]
610
611
        [setattr(self, i, getattr(tile1, i)) for i in tile1.__dict__
         if not callable(getattr(tile1, i)) and not isinstance(getattr(tile1, i), (np.ndarray, GeoArray))]
612
613

        # MERGE ARRAY-ATTRIBUTES
614
        list_arraynames = [i for i in tile1.__dict__ if not callable(getattr(tile1, i)) and
615
                           isinstance(getattr(tile1, i), (np.ndarray, GeoArray))]
616
617
        list_arraynames = ['_arr'] + [i for i in list_arraynames if
                                      i != '_arr']  # list must start with _arr, otherwise setters will not work
618
619
620
621

        for arrname in list_arraynames:
            samplearray = getattr(tile1, arrname)
            assert isinstance(samplearray, (np.ndarray, GeoArray)), \
622
                'Received a %s object for attribute %s. Expected a numpy array or an instance of GeoArray.' \
623
                % (type(samplearray), arrname)
624
625
            is_3d = samplearray.ndim == 3
            bands = (samplearray.shape[2],) if is_3d else ()  # dynamic -> works for arr, cld_arr,...
626
627
628
629
            target_shape = tuple(self.shape_fullArr[:2]) + bands
            target_dtype = samplearray.dtype
            merged_array = self._numba_array_merger(list_GMS_tiles, arrname, target_shape, target_dtype)

630
631
            setattr(self, arrname if not arrname.startswith('_') else arrname[1:],
                    merged_array)  # use setters if possible
632
633
634
635
            # NOTE: this asserts that each attribute starting with '_' has also a property with a setter!

        # UPDATE ARRAY-DEPENDENT ATTRIBUTES
        self.arr_shape = 'cube'
636
        self.arr_pos = None
637
638
639

        # update MetaObj attributes
        self.meta_odict.update({
640
641
            'samples': self.arr.cols, 'lines': self.arr.rows, 'bands': self.arr.bands,
            'map info': geotransform2mapinfo(self.arr.gt, self.arr.prj), 'coordinate system string': self.arr.prj, })
642
643
644
645

        # calculate data_corners_imXY (mask_nodata is always an array here because get_mapPos always returns an array)
        corners_imYX = calc_FullDataset_corner_positions(
            self.mask_nodata, assert_four_corners=False, algorithm='shapely')
646
        self.trueDataCornerPos = [(YX[1], YX[0]) for YX in corners_imYX]  # [UL, UR, LL, LR]
647
648
649
650
651
652

        # calculate trueDataCornerLonLat
        data_corners_LatLon = pixelToLatLon(self.trueDataCornerPos, geotransform=self.arr.gt, projection=self.arr.prj)
        self.trueDataCornerLonLat = [(YX[1], YX[0]) for YX in data_corners_LatLon]

        # calculate trueDataCornerUTM
653
654
        data_corners_utmYX = pixelToMapYX(self.trueDataCornerPos, geotransform=self.arr.gt,
                                          projection=self.arr.prj)  # FIXME asserts gt in UTM coordinates
655
656
657
658
        self.trueDataCornerUTM = [(YX[1], YX[0]) for YX in data_corners_utmYX]

        return copy.copy(self)

659
660
661
    @staticmethod
    @jit
    def _numba_array_merger(list_GMS_tiles, arrname2merge, target_shape, target_dtype):
Daniel Scheffler's avatar
Daniel Scheffler committed
662
        # type: (list, str, tuple, np.dtype) -> np.ndarray
663
664
665
666
667
668
669
670
671
        """
        private function, e.g. called by merge_GMS_tiles_to_GMS_obj() in order to fasten array merging

        :param list_GMS_tiles:
        :param arrname2merge:
        :param target_shape:
        :param target_dtype:
        :return:
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
672

673
674
675
676
677
678
679
        out_arr = np.empty(target_shape, dtype=target_dtype)
        for idx, tile in enumerate(list_GMS_tiles):
            rowStart, rowEnd = tile.arr_pos[0]
            colStart, colEnd = tile.arr_pos[1]
            out_arr[rowStart:rowEnd + 1, colStart:colEnd + 1] = getattr(tile, arrname2merge)
        return out_arr

Daniel Scheffler's avatar
Daniel Scheffler committed
680
    def log_for_fullArr_or_firstTile(self, log_msg, subset=None):
681
682
683
684
685
686
687
        """Send a message to the logger only if full array or the first tile is currently processed.
        This function can be called when processing any tile but log message will only be sent from first tile.

        :param log_msg:  the log message to be logged
        :param subset:   subset argument as sent to e.g. DN2TOARadRefTemp that indicates which tile is to be processed.
                         Not needed if self.arr_pos is not None.
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
688

689
690
691
692
        if subset is None and \
            (self.arr_shape == 'cube' or self.arr_pos is None or [self.arr_pos[0][0], self.arr_pos[1][0]] == [0, 0]) or\
                subset == ['cube', None] or (subset and [subset[1][0][0], subset[1][1][0]] == [0, 0]) or \
                hasattr(self, 'logAtThisTile') and getattr(self, 'logAtThisTile'):  # cube or 1st tile
Daniel Scheffler's avatar
Daniel Scheffler committed
693
            self.logger.info(log_msg)
694
695
696
697
        else:
            pass

    def apply_nodata_mask_to_ObjAttr(self, attrname, out_nodata_val=None):
698
        # type: (str,int) -> None
699
        """Applies self.mask_nodata to the specified array attribute by setting all values where mask_nodata is 0 to the
700
701
702
703
704
705
706
        given nodata value.

        :param attrname:         The attribute to apply the nodata mask to. Must be an array attribute or
                                 a string path to a previously saved ENVI-file.
        :param out_nodata_val:   set the values of the given attribute to this value.
        """

707
        assert hasattr(self, attrname)
708

709
        if getattr(self, attrname) is not None:
710

711
712
713
            if isinstance(getattr(self, attrname), str):
                update_spec_vals = True if attrname == 'arr' else False
                self.apply_nodata_mask_to_saved_ENVIfile(getattr(self, attrname), out_nodata_val, update_spec_vals)
714
            else:
715
                assert isinstance(getattr(self, attrname), (np.ndarray, GeoArray)), \
716
                    'L1A_obj.%s must be a numpy array or an instance of GeoArray. Got type %s.' \
717
718
                    % (attrname, type(getattr(self, attrname)))
                assert hasattr(self, 'mask_nodata') and self.mask_nodata is not None
719

720
                self.log_for_fullArr_or_firstTile('Applying nodata mask to L1A_object.%s...' % attrname)
721

722
                nodata_val = out_nodata_val if out_nodata_val else \
723
                    DEF_D.get_outFillZeroSaturated(getattr(self, attrname).dtype)[0]
724
                getattr(self, attrname)[self.mask_nodata.astype(np.int8) == 0] = nodata_val
725

726
727
                if attrname == 'arr':
                    self.MetaObj.spec_vals['fill'] = nodata_val
728
729
730
731
732
733

    def build_combined_masks_array(self):
        # type: () -> dict
        """Generates self.masks attribute (unsigned integer 8bit) from by concatenating all masks included in GMS obj.
        The corresponding metadata is assigned to L1A_obj.masks_meta. Empty mask attributes are skipped."""

734
        arrays2combine = [aN for aN in ['mask_nodata', 'mask_clouds']
735
                          if hasattr(self, aN) and isinstance(getattr(self, aN), (GeoArray, np.ndarray))]
736
737
        if arrays2combine:
            self.log_for_fullArr_or_firstTile('Combining masks...')
738
739

            def get_data(arrName): return getattr(self, arrName).astype(np.uint8)[:, :, None]
740
741

            for aN in arrays2combine:
742
                if False in np.equal(getattr(self, aN), getattr(self, aN).astype(np.uint8)):
743
744
745
                    warnings.warn('Some pixel values of attribute %s changed during data type '
                                  'conversion within build_combined_masks_array().')

746
            # set self.masks
747
748
749
            self.masks = get_data(arrays2combine[0]) if len(arrays2combine) == 1 else \
                np.concatenate([get_data(aN) for aN in arrays2combine], axis=2)
            self.masks.bandnames = arrays2combine  # set band names of GeoArray (allows later indexing by band name)
750

751
            # set self.masks_meta
752
            nodataVal = DEF_D.get_outFillZeroSaturated(self.masks.dtype)[0]
753
            self.masks_meta = {'map info': self.MetaObj.map_info, 'coordinate system string': self.MetaObj.projection,
754
755
                               'bands': len(arrays2combine), 'band names': arrays2combine,
                               'data ignore value': nodataVal}
756
757

            return {'desc': 'masks', 'row_start': 0, 'row_end': self.shape_fullArr[0],
758
                    'col_start': 0, 'col_end': self.shape_fullArr[1], 'data': self.masks}  # usually not needed
759
760

    def apply_nodata_mask_to_saved_ENVIfile(self, path_saved_ENVIhdr, custom_nodata_val=None, update_spec_vals=False):
761
        # type: (str,int,bool) -> None
762
763
        """Applies self.mask_nodata to a saved ENVI file with the same X/Y dimensions like self.mask_nodata by setting all
         values where mask_nodata is 0 to the given nodata value.
764
765
766
767
768
769
770
771

        :param path_saved_ENVIhdr:  <str> The path of the ENVI file to apply the nodata mask to.
        :param custom_nodata_val:   <int> set the values of the given attribute to this value.
        :param update_spec_vals:    <bool> whether to update self.MetaObj.spec_vals['fill']
        """

        self.log_for_fullArr_or_firstTile('Applying nodata mask to saved ENVI file...')
        assert os.path.isfile(path_saved_ENVIhdr)
772
773
774
        assert hasattr(self, 'mask_nodata') and self.mask_nodata is not None
        if not path_saved_ENVIhdr.endswith('.hdr') and os.path.isfile(os.path.splitext(path_saved_ENVIhdr)[0] + '.hdr'):
            path_saved_ENVIhdr = os.path.splitext(path_saved_ENVIhdr)[0] + '.hdr'
775
        if custom_nodata_val is None:
776
            dtype_IDL = int(INP_R.read_ENVIhdr_to_dict(path_saved_ENVIhdr)['data type'])
777
            nodata_val = DEF_D.get_outFillZeroSaturated(DEF_D.dtype_lib_IDL_Python[dtype_IDL])[0]
778
779
        else:
            nodata_val = custom_nodata_val
780
        FileObj = spectral.open_image(path_saved_ENVIhdr)
781
        File_memmap = FileObj.open_memmap(writable=True)
782
        File_memmap[self.mask_nodata == 0] = nodata_val
783
784
        if update_spec_vals:
            self.MetaObj.spec_vals['fill'] = nodata_val
785
786

    def combine_tiles_to_ObjAttr(self, tiles, target_attr):
787
        # type: (list,str) -> None
788
789
790
791
792
793
794
        """Combines tiles, e.g. produced by L1A_P.L1A_object.DN2TOARadRefTemp() to a single attribute.
        If CFG.usecase.CFG.job.exec_mode == 'Python' the produced attribute is additionally written to disk.

        :param tiles:           <list> a list of dictionaries with the keys 'desc', 'data', 'row_start','row_end',
                                'col_start' and 'col_end'
        :param target_attr:     <str> the name of the attribute to be produced
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
795

796
        warnings.warn("'combine_tiles_to_ObjAttr' is deprecated.", DeprecationWarning)
797
        assert tiles[0] and isinstance(tiles, list) and isinstance(tiles[0], dict), \
798
799
800
            "The 'tiles' argument has to be list of dictionaries with the keys 'desc', 'data', 'row_start'," \
            "'row_end', 'col_start' and 'col_end'."
        self.logger.info("Building L1A attribute '%s' by combining given tiles..." % target_attr)
801
802
        tiles = [tiles] if not isinstance(tiles, list) else tiles
        sampleTile = dict(tiles[0])
803
804
805
        target_shape = self.shape_fullArr if len(sampleTile['data'].shape) == 3 else self.shape_fullArr[:2]
        setattr(self, target_attr, np.empty(target_shape, dtype=sampleTile['data'].dtype))
        for tile in tiles:
806
807
808
809
810
            rS, rE, cS, cE = tile['row_start'], tile['row_end'], tile['col_start'], tile['col_end']
            if len(target_shape) == 3:
                getattr(self, target_attr)[rS:rE + 1, cS:cE + 1, :] = tile['data']
            else:
                getattr(self, target_attr)[rS:rE + 1, cS:cE + 1] = tile['data']
811
        if target_attr == 'arr':
812
813
            self.arr_desc = sampleTile['desc']
            self.arr_shape = 'cube' if len(self.arr.shape) == 3 else 'band' if len(self.arr.shape) == 2 else 'unknown'
814

815
816
            if CFG.job.exec_mode == 'Python':  # and not 'Flink'
                path_radref_file = os.path.join(self.ExtractedFolder, self.baseN + '__' + self.arr_desc)
817
818
819
820
821
822
823
824
825
826
827
                # path_radref_file = os.path.abspath('./testing/out/%s_TOA_Ref' % self.baseN)
                while not os.path.isdir(os.path.dirname(path_radref_file)):
                    try:
                        os.makedirs(os.path.dirname(path_radref_file))
                    except OSError as e:  # maybe not neccessary anymore in python 3
                        if not e.errno == 17:
                            raise
                GEOP.ndarray2gdal(self.arr, path_radref_file, importFile=self.MetaObj.Dataname, direction=3)
                self.MetaObj.Dataname = path_radref_file

    def write_tiles_to_ENVIfile(self, tiles, overwrite=True):
828
        # type: (list,bool) -> None
829
830
831
832
833
834
835
836
        """Writes tiles, e.g. produced by L1A_P.L1A_object.DN2TOARadRefTemp() to a single output ENVI file.

        :param tiles:           <list> a list of dictionaries with the keys 'desc', 'data', 'row_start','row_end',
                                'col_start' and 'col_end'
        :param overwrite:       whether to overwrite files that have been produced earlier
        """

        self.logger.info("Writing tiles '%s' temporarily to disk..." % tiles[0]['desc'])
837
838
839
840
        outpath = os.path.join(self.ExtractedFolder, '%s__%s.%s' % (self.baseN, tiles[0]['desc'], self.outInterleave))
        if CFG.usecase.conversion_type_optical in tiles[0]['desc'] or \
           CFG.usecase.conversion_type_thermal in tiles[0]['desc']:
            self.meta_odict = self.MetaObj.to_odict()  # important in order to keep geotransform/projection
841
842
843
844
            self.arr_desc = tiles[0]['desc']
            self.arr = outpath
            # self.arr = os.path.abspath('./testing/out/%s_TOA_Ref.bsq' % self.baseN)
            self.MetaObj.Dataname = self.arr
845
            self.arr_shape = \
846
847
                'cube' if len(tiles[0]['data'].shape) == 3 else 'band' if len(
                    tiles[0]['data'].shape) == 2 else 'unknown'
848
849
850
        elif tiles[0]['desc'] == 'masks':
            self.masks = outpath
        elif tiles[0]['desc'] == 'lonlat_arr':
851
852
853
854
855
            # outpath = os.path.join(os.path.abspath('./testing/out/'),'%s__%s.%s'
            #     %(self.baseN, tiles[0]['desc'], self.outInterleave))
            self.lonlat_arr = outpath  # FIXME
        outpath = os.path.splitext(outpath)[0] + '.hdr' if not outpath.endswith('.hdr') else outpath
        out_shape = self.shape_fullArr[:2] + ([tiles[0]['data'].shape[2]] if len(tiles[0]['data'].shape) == 3 else [1])
856
857
        OUT_W.Tiles_Writer(tiles, outpath, out_shape, tiles[0]['data'].dtype, self.outInterleave, self.meta_odict,
                           overwrite=overwrite)
858
859

    def to_MGRS_tiles(self, pixbuffer=10, v=False):
860
        # type: (int) -> self
861
        """Returns a generator object where items represent the MGRS tiles for the GMS object.
862
863
864
865
866

        :param pixbuffer:   <int> a buffer in pixel values used to generate an overlap between the returned MGRS tiles
        :param v:           <bool> verbose mode
        :return:            <list> of MGRS_tile objects
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
867

868
        assert self.arr_shape == 'cube', "Only 'cube' objects can be cut into MGRS tiles. Got %s." % self.arr_shape
869
        self.logger.info('Cutting scene %s (entity ID %s) into MGRS tiles...' % (self.scene_ID, self.entity_ID))
870

871
872
        # get GeoDataFrame containing all overlapping MGRS tiles
        # (MGRS geometries completely within nodata area are excluded)
Daniel Scheffler's avatar
Daniel Scheffler committed
873
        GDF_MGRS_tiles = DB_T.get_overlapping_MGRS_tiles(CFG.job.conn_database,
874
                                                         tgt_corners_lonlat=self.trueDataCornerLonLat)
875
876

        # calculate image coordinate bounds of the full GMS object for each MGRS tile within the GeoDataFrame
877
878
879
880
        gt = mapinfo2geotransform(self.meta_odict['map info'])

        def get_arrBounds(MGRStileObj):
            return list(np.array(MGRStileObj.poly_utm.buffer(pixbuffer * gt[1]).bounds)[[0, 2, 1, 3]])
881

882
883
        GDF_MGRS_tiles['MGRStileObj'] = list(GDF_MGRS_tiles['granuleid'].map(lambda mgrsTileID: MGRS_tile(mgrsTileID)))
        GDF_MGRS_tiles['map_bounds_MGRS'] = list(GDF_MGRS_tiles['MGRStileObj'].map(get_arrBounds))  # xmin,xmax,ymin,yma
884
885
886
887
888
889

        # find first tile to log and assign 'logAtThisTile' later
        dictIDxminymin = {(b[0] + b[2]): ID for ID, b in
                          zip(GDF_MGRS_tiles['granuleid'], GDF_MGRS_tiles['map_bounds_MGRS'])}
        firstTile_ID = dictIDxminymin[min(dictIDxminymin.keys())]

890
        # ensure self.masks exists (does not exist in Flink mode because in that case self.fill_from_disk() is skipped)
891
892
893
        if not hasattr(self, 'masks') or self.masks is None:
            self.build_combined_masks_array()  # creates self.masks and self.masks_meta

894
895
896
        # read whole dataset into RAM in order to fasten subsetting
        self.arr.to_mem()
        self.masks.to_mem()  # to_mem ensures that the whole dataset is present in memory
897
898
899

        # produce data for each MGRS tile in loop
        for GDF_idx, GDF_row in GDF_MGRS_tiles.iterrows():
900
901
902
903
904
905
906
            tileObj = self.get_subset_obj(mapBounds=GDF_row.map_bounds_MGRS,
                                          mapBounds_prj=EPSG2WKT(GDF_row['MGRStileObj'].EPSG),
                                          out_prj=EPSG2WKT(GDF_row['MGRStileObj'].EPSG),
                                          logmsg='Producing MGRS tile %s from scene %s (entity ID %s).'
                                                 % (GDF_row.granuleid, self.scene_ID, self.entity_ID),
                                          v=v)
            MGRS_tileID = GDF_row['granuleid']
907
908
909
910
911
912
913
914

            # set MGRS info
            tileObj.arr_shape = 'MGRS_tile'
            tileObj.MGRS_info = {'tile_ID': MGRS_tileID, 'grid1mil': MGRS_tileID[:3], 'grid100k': MGRS_tileID[3:]}

            # set logAtThisTile
            tileObj.logAtThisTile = MGRS_tileID == firstTile_ID

915
916
917
918
            # close logger of tileObj and of self in order to avoid logging permission errors
            tileObj.close_GMS_loggers()
            self.close_GMS_loggers()

919
920
921
            yield tileObj

        # set array attributes back to file path if they had been a filePath before
922
923
924
925
        if self.arr.filePath:
            self.arr.to_disk()
        if self.masks.filePath:
            self.masks.to_disk()
926
927
928
929

    def to_GMS_file(self, path_gms_file=None):
        self.close_GMS_loggers()

930
931
932
        # make sure meta_odict is present
        self.meta_odict = self.meta_odict

933
        dict2write = self.attributes2dict(remove_privates=True)
934
        dict2write['arr_shape'], dict2write['arr_pos'] = ['cube', None]
935
        path_gms_file = path_gms_file if path_gms_file else self.pathGen.get_path_gmsfile()
936
937
938

        for k, v in list(dict2write.items()):
            # if isinstance(v,np.ndarray) or isinstance(v,dict) or hasattr(v,'__dict__'):
939
940
            # so, wenn meta-dicts nicht ins gms-file sollen. ist aber vllt ni schlecht
            # -> erlaubt lesen der metadaten direkt aus gms
941

942
943
            if k == 'MetaObj':
                continue  # make sure MetaObj getter is not called -> would delete meta_odict
944
            elif isinstance(v, datetime.datetime):
945
                dict2write[k] = v.strftime('%Y-%m-%d %H:%M:%S.%f%z')  # FIXME
946
            elif isinstance(v, DatasetLogger):
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
                if hasattr(v, 'handlers') and v.handlers[:]:
                    warnings.warn('Not properly closed logger at GMS_obj.logger pointing to %s.' % v.path_logfile)
                dict2write[k] = 'not set'
            elif isinstance(v, collections.OrderedDict) or isinstance(v, dict):
                dict2write[k] = dict2write[k].copy()
                if 'logger' in v:
                    if hasattr(dict2write[k]['logger'], 'handlers') and dict2write[k]['logger'].handlers[:]:
                        warnings.warn("Not properly closed logger at %s['logger'] pointing to %s."
                                      % (k, dict2write[k]['logger'].path_logfile))
                    dict2write[k]['logger'] = 'not set'
            elif isinstance(v, np.ndarray):
                # delete every 3D Array larger than 100x100
                if len(v.shape) == 2 and sum(v.shape) <= 200:
                    dict2write[k] = v.tolist()  # numpy arrays are not jsonable
                else:
                    del dict2write[k]
            elif hasattr(v, '__dict__'):
                # löscht Instanzen von Objekten und Arrays, aber keine OrderedDicts
                if hasattr(v, 'logger'):
                    if hasattr(dict2write[k].logger, 'handlers') and dict2write[k].logger.handlers[:]:
                        warnings.warn("Not properly closed logger at %s.logger pointing to %s."
                                      % (k, dict2write[k].logger.path_logfile))
                    dict2write[k].logger = 'not set'
                del dict2write[k]

                # class customJSONEncoder(json.JSONEncoder):
                #     def default(self, obj):
                #         if isinstance(obj, np.ndarray):
                #             return '> numpy array <'
                #         if isinstance(obj, dict): # funktioniert nicht
                #             return '> python dictionary <'
                #         if hasattr(obj,'__dict__'):
                #             return '> python object <'
                #         # Let the base class default method raise the TypeError
                #         return json.JSONEncoder.default(self, obj)
982
983
                # json.dump(In, open(path_out_baseN,'w'), skipkeys=True,
                #           sort_keys=True, cls=customJSONEncoder, separators=(',', ': '), indent =4)
984
985
986
        with open(path_gms_file, 'w') as outF:
            json.dump(dict2write