L1A_P.py 38.2 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2
3
4

# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data
#
5
# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6
7
8
9
10
11
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
12
# the terms of the GNU General Public License as published by the Free Software
13
14
# Foundation, either version 3 of the License, or (at your option) any later version.
# Please note the following exception: `gms_preprocessing` depends on tqdm, which
15
16
17
# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
18
19
20
21
22
23
24
25
26
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

27
import datetime
28
29
import os
import re
30
import warnings
31

32
33
34
import gdal
import gdalnumeric
import numpy as np
35
from natsort import natsorted
36

37
38
39
40
try:
    from pyhdf import SD
except ImportError:
    SD = None
41

42
from geoarray import GeoArray
43
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
44
from py_tools_ds.geo.coord_trafo import pixelToLatLon, imYX2mapYX
45
from py_tools_ds.geo.map_info import mapinfo2geotransform
46
from py_tools_ds.geo.projection import EPSG2WKT, isProjectedOrGeographic
47

48
from ..options.config import GMS_config as CFG
49
from . import geoprocessing as GEOP
50
from ..misc.definition_dicts import get_outFillZeroSaturated, is_dataset_provided_as_fullScene
51
from ..misc.locks import IOLock
52
from ..model.gms_object import GMS_object
53
from ..model import metadata as META
54

55
56
__author__ = 'Daniel Scheffler'

Daniel Scheffler's avatar
Daniel Scheffler committed
57

58
class L1A_object(GMS_object):
59
    """Features input reader and raster-/metadata homogenization."""
60

61
    def __init__(self, image_type='', satellite='', sensor='', subsystem='', sensormode='', acq_datetime=None,
62
                 entity_ID='', scene_ID=-9999, filename='', dataset_ID=-9999, proc_status='', **kwargs):
63
        """:param :  instance of gms_object.GMS_object or None
64
        """
65
66
        # TODO docstring

67
68
        # NOTE: kwargs is in here to allow instanciating with additional arg 'proc_level'

69
        # load default attribute values and methods
70
        super(L1A_object, self).__init__()
71
72

        # unpack kwargs
73
74
75
76
77
78
        self.proc_level = 'L1A'
        self.image_type = image_type  # FIXME not needed anymore?
        self.satellite = satellite
        self.sensor = sensor
        self.subsystem = subsystem
        self.sensormode = sensormode
79
        self.acq_datetime = acq_datetime
80
81
82
83
        self.entity_ID = entity_ID
        self.scene_ID = scene_ID
        self.filename = filename
        self.dataset_ID = dataset_ID
84
85
86

        # set pathes (only if there are valid init args)
        if image_type and satellite and sensor:
87
            self.set_pathes()  # also overwrites logfile
88
89
90
91
92
93
            self.validate_pathes()

            if self.path_archive_valid:
                self.logger.info('L1A object for %s %s%s (entity-ID %s) successfully initialized.'
                                 % (self.satellite, self.sensor,
                                    (' ' + self.subsystem) if self.subsystem not in [None, ''] else '', self.entity_ID))
94

95
96
97
98
99
100
        # (re)set the processing status
        if self.scene_ID in self.proc_status_all_GMSobjs:
            del self.proc_status_all_GMSobjs[self.scene_ID]

        self.proc_status = proc_status or 'initialized'  # if proc_status = 'running' is given by L1A_map

101
    def import_rasterdata(self):
102
        if re.search(r"ALOS", self.satellite, re.I):
103
104
105
            '''First 22 lines are nodata: = maybe due to an issue of the GDAL CEOS driver.
            But: UL of metadata refers to [row =0, col=21]! So the imported GeoTransform is correct when
            the first 21 columns are deleted.'''
106
107
            self.archive_to_rasObj(self.path_archive, self.path_InFilePreprocessor,
                                   subset=['block', [[None, None], [21, None]]])
108
        elif re.search(r"Terra", self.satellite, re.I):
109
            self.ASTER_HDF_to_rasObj(self.path_archive, path_output=self.path_InFilePreprocessor)
110
        else:
111
            self.archive_to_rasObj(self.path_archive, path_output=self.path_InFilePreprocessor)
112
113
114

        # set shape_fullArr
        self.shape_fullArr = self.arr.shape
115

116
    def archive_to_rasObj(self, path_archive, path_output=None, subset=None):
117
        from ..misc.helper_functions import convert_absPathArchive_to_GDALvsiPath
118

119
        assert subset is None or isinstance(subset, list) and len(subset) == 2, \
Daniel Scheffler's avatar
Daniel Scheffler committed
120
            "Subset argument has be a list with 2 elements."
121
122
        if subset:
            assert subset[0] == 'block',\
123
124
                "The function 'archive_to_rasObj' only supports block subsetting. Attempted to perform '%s' " \
                "subsetting." % subset[0]
125

126
        self.logger.info('Reading %s %s %s image data...' % (self.satellite, self.sensor, self.subsystem))
127
        gdal_path_archive = convert_absPathArchive_to_GDALvsiPath(path_archive)
128
        project_dir = os.path.abspath(os.curdir)
129
        os.chdir(os.path.dirname(path_archive))
130
        files_in_archive = gdal.ReadDirRecursive(gdal_path_archive)  # needs ~12sek for Landsat-8
131
132
        assert files_in_archive, 'No files in archive %s for scene %s (entity ID: %s). ' \
                                 % (os.path.basename(path_archive), self.scene_ID, self.entity_ID)
133
        full_LayerBandsAssignment = META.get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False)
134

135
136
137
138
        ####################################################
        # get list of raster files to be read from archive #
        ####################################################

139
140
        image_files = []
        is_ALOS_Landsat_S2 = \
141
142
            re.search(r'ALOS', self.satellite) or re.search(r'Landsat', self.satellite) or \
            re.search(r'Sentinel-2', self.satellite)
143
        n_files2search = len(full_LayerBandsAssignment) if is_ALOS_Landsat_S2 else 1
144

145
        for File in natsorted(files_in_archive):
146
            search_res = \
147
148
149
150
151
152
                re.search(r"IMG-0[0-9]-[\s\S]*", File) if re.search(r'ALOS', self.satellite) else \
                re.search(r"[\S]*_B[1-9][0-9]?[\S]*.TIF", File) if re.search(r'Landsat', self.satellite) else \
                re.search(r"[0-9]*.tif", File) if re.search(r'RapidEye', self.satellite) else \
                re.search(r"imagery.tif", File) if re.search(r'SPOT', self.satellite) else \
                re.search(r"[\S]*.SAFE/GRANULE/%s/IMG_DATA/[\S]*_B[0-9][\S]*.jp2"
                          % self.entity_ID, File) if re.search(r'Sentinel-2', self.satellite) else None
153

Daniel Scheffler's avatar
Daniel Scheffler committed
154
            if search_res:
155
                if re.search(r'Sentinel-2', self.satellite):
156
                    # add only those files that are corresponding to subsystem (e.g. S2A10: fullLBA = ['2','3','4','8'])
157
                    if 1 in [1 if re.search(r"[\S]*_B[0]?%s.jp2" % LBAn, os.path.basename(File)) else 0
158
159
160
161
                             for LBAn in full_LayerBandsAssignment]:
                        image_files.append(File)
                else:
                    image_files.append(File)
Daniel Scheffler's avatar
Daniel Scheffler committed
162
163

        # create and fill raster object
164
        if n_files2search > 1:
165
166
167
            #####################################
            # validate number of expected files #
            #####################################
168

169
            if re.search(r'ETM+', self.sensor) and self.acq_datetime > datetime.datetime(year=2003, month=5, day=31):
170
                expected_files_count = 2 * len(full_LayerBandsAssignment)
171
172
            else:
                expected_files_count = len(full_LayerBandsAssignment)
173

174
            assert len(image_files) == expected_files_count, "Expected %s image files in '%s'. Found %s." \
175
176
                                                             % (len(full_LayerBandsAssignment), path_archive,
                                                                len(image_files))
177

178
179
180
            ###############################
            # get paths of files to stack #
            ###############################
181

182
183
184
            # NOTE: image_files is a SORTED list of image filenames; self.LayerBandsAssignment may be sorted by CWL
            filtered_files = []
            for bN in self.LayerBandsAssignment:  # unsorted, e.g., ['1', '2', '3', '4', '5', '9', '6', '7']
185
                for fN, b in zip(image_files, natsorted(full_LayerBandsAssignment)):  # both sorted nicely
186
187
188
                    if b == bN:
                        filtered_files.append(fN)

189
            paths_files2stack = [os.path.join(gdal_path_archive, i) for i in filtered_files]
190
191
192
193

            #########################
            # read the raster data  #
            #########################
194

195
            rasObj = GEOP.GEOPROCESSING(paths_files2stack[0], self.logger)
196

197
            # in case a subset is to be read: prepare rasObj instance to read a subset
Daniel Scheffler's avatar
Daniel Scheffler committed
198
            if subset:
199
200
201
202
203
                full_dim = [0, rasObj.rowEnd, 0, rasObj.colEnd]
                sub_dim = [subset[1][0][0], subset[1][0][1], subset[1][1][0], subset[1][0][1]]
                sub_dim = [sub_dim[i] if sub_dim[i] else full_dim[i] for i in range(len(sub_dim))]
                subset = ['block', [[sub_dim[0], sub_dim[1] + 1], [sub_dim[2], sub_dim[3] + 1]]]
                rasObj = GEOP.GEOPROCESSING(paths_files2stack[0], self.logger, subset=subset)
204

205
            # perform layer stack
206
            with IOLock(allowed_slots=CFG.max_parallel_reads_writes, logger=self.logger):
207
208
209
210
211
212
                if CFG.inmem_serialization and path_output is None:  # numpy array output
                    self.arr = rasObj.Layerstacking(paths_files2stack)
                    self.path_InFilePreprocessor = paths_files2stack[0]
                else:  # 'MEMORY' or physical output
                    rasObj.Layerstacking(paths_files2stack, path_output=path_output)  # writes an output (gdal_merge)
                    self.arr = path_output
213

214
        else:
215
            assert image_files != [], 'No image files found within the archive matching the expected naming scheme.'
216
            path_file2load = os.path.join(gdal_path_archive, image_files[0])
Daniel Scheffler's avatar
Daniel Scheffler committed
217
            rasObj = GEOP.GEOPROCESSING(path_file2load, self.logger)
218

Daniel Scheffler's avatar
Daniel Scheffler committed
219
            if subset:
220
221
                full_dim = [0, rasObj.rowEnd, 0, rasObj.colEnd]
                sub_dim = [subset[1][0][0], subset[1][0][1], subset[1][1][0], subset[1][0][1]]
Daniel Scheffler's avatar
Daniel Scheffler committed
222
                sub_dim = [sub_dim[i] if sub_dim[i] else full_dim[i] for i in range(len(sub_dim))]
223
224
                subset = ['block', [[sub_dim[0], sub_dim[1] + 1], [sub_dim[2], sub_dim[3] + 1]]]
                rasObj = GEOP.GEOPROCESSING(path_file2load, self.logger, subset=subset)
225

226
            # read a single file
227
            with IOLock(allowed_slots=CFG.max_parallel_reads_writes, logger=self.logger):
228
229
230
231
232
233
234
235
                if CFG.inmem_serialization and path_output is None:  # numpy array output
                    self.arr = gdalnumeric.LoadFile(path_file2load) if subset is None else \
                        gdalnumeric.LoadFile(path_file2load, rasObj.colStart, rasObj.rowStart, rasObj.cols, rasObj.rows)
                    self.path_InFilePreprocessor = path_file2load
                else:  # 'MEMORY' or physical output
                    GEOP.ndarray2gdal(rasObj.tondarray(), path_output,
                                      geotransform=rasObj.geotransform, projection=rasObj.projection)
                    self.arr = path_output
236

237
238
        os.chdir(project_dir)

239
240
241
    def ASTER_HDF_to_rasObj(self, path_archive, path_output=None):
        subsystem_identifier = 'VNIR' if self.subsystem in ['VNIR1', 'VNIR2'] else 'SWIR' \
            if self.subsystem == 'SWIR' else 'TIR'
242

243
        ds = gdal.Open(path_archive)
244

Daniel Scheffler's avatar
Daniel Scheffler committed
245
        if ds:
246
            sds_md = ds.GetMetadata('SUBDATASETS')
247
            data_arr = None
248
            for bidx, b in enumerate(self.LayerBandsAssignment):
249
250
251
                sds_name = [i for i in sds_md.values() if '%s_Band%s:ImageData' % (subsystem_identifier, b) in str(i) or
                            '%s_Swath:ImageData%s' % (subsystem_identifier, b) in str(i)][0]
                data = gdalnumeric.LoadFile(sds_name)
252
253
                if bidx == 0:
                    data_arr = np.empty(data.shape + (len(self.LayerBandsAssignment),), data.dtype)
254
                data_arr[:, :, bidx] = data
255

256
            if CFG.inmem_serialization and path_output is None:  # numpy array output
257
258
                self.arr = data_arr
            else:
259
260
                GEOP.ndarray2gdal(data_arr, path_output, geotransform=ds.GetGeoTransform(),
                                  projection=ds.GetProjection(), direction=3)
261
                self.arr = path_output
262
263

        elif SD is not None:
264
            self.logger.info('Missing HDF4 support within GDAL. Reading HDF file using alternative reader.')
265
            hdfFile = SD.SD(path_archive, SD.SDC.READ)
266
            i, list_matching_dsIdx = 0, []
267

268
            while True:  # Get dataset indices within HDF file
269
                # noinspection PyBroadException
270
271
272
273
                try:
                    ds = hdfFile.select(i)
                    if subsystem_identifier in str(ds.dimensions()) and 'ImagePixel' in str(ds.dimensions()):
                        list_matching_dsIdx.append(i)
274
                    i += 1
275
                except Exception:
276
277
                    break

278
279
            list_matching_dsIdx = list_matching_dsIdx[:3] if self.subsystem == 'VNIR1' else \
                [list_matching_dsIdx[-1]] if self.subsystem == 'VNIR2' else list_matching_dsIdx
280
            data_arr = None
281
282
            for i, dsIdx in enumerate(list_matching_dsIdx):
                data = hdfFile.select(dsIdx)[:]
283
284
                if i == 0:
                    data_arr = np.empty(data.shape + (len(self.LayerBandsAssignment),), data.dtype)
285
                data_arr[:, :, i] = data
286

287
            if CFG.inmem_serialization and path_output is None:  # numpy array output
Daniel Scheffler's avatar
Daniel Scheffler committed
288
                self.arr = data_arr
289
            else:
290
                GEOP.ndarray2gdal(data_arr, path_output, direction=3)
291
                self.arr = path_output
292

293
294
295
296
        else:
            self.logger.error('Missing HDF4 support. Reading HDF file failed.')
            raise ImportError('No suitable library for reading HDF4 data available.')

297
        del ds
298

299
    def import_metadata(self):
300
301
302
        """Reads metainformation of the given file from the given ASCII metafile.
        Works for: RapidEye (metadata.xml),SPOT(metadata.dim),LANDSAT(mtl.txt),ASTER(downloaded coremetadata),
        ALOS(summary.txt & Leader file)
303
304

        :return:
305
        """
306

307
        self.logger.info('Reading %s %s %s metadata...' % (self.satellite, self.sensor, self.subsystem))
308
309
310
        self.MetaObj = META.METADATA(self.GMS_identifier)
        self.MetaObj.read_meta(self.scene_ID, self.path_InFilePreprocessor,
                               self.path_MetaPreprocessor, self.LayerBandsAssignment)
311

312
313
        self.logger.debug("The following metadata have been read:")
        [self.logger.debug("%20s : %-4s" % (key, val)) for key, val in self.MetaObj.overview.items()]
314

315
        # set some object attributes directly linked to metadata
316
        self.subsystem = self.MetaObj.Subsystem
317
        self.arr.nodata = self.MetaObj.spec_vals['fill']
318

319
320
321
        # update acquisition date to a complete datetime object incl. date, time and timezone
        if self.MetaObj.AcqDateTime:
            self.acq_datetime = self.MetaObj.AcqDateTime
322

323
        # set arr_desc
324
325
326
        self.arr_desc = \
            'DN' if self.MetaObj.PhysUnit == 'DN' else \
            'Rad' if self.MetaObj.PhysUnit == "W * m-2 * sr-1 * micrometer-1" else \
327
328
329
            'TOA_Ref' if re.search(r'TOA_Reflectance', self.MetaObj.PhysUnit, re.I) else \
            'BOA_Ref' if re.search(r'BOA_Reflectance', self.MetaObj.PhysUnit, re.I) else \
            'Temp' if re.search(r'Degrees Celsius', self.MetaObj.PhysUnit, re.I) else None
330
331

        assert self.arr_desc, 'GMS_obj contains an unexpected physical unit: %s' % self.MetaObj.PhysUnit
332
333

    def calc_TOARadRefTemp(self, subset=None):
334
        """Convert DN, Rad or TOA_Ref data to TOA Reflectance, to Radiance or to Surface Temperature
335
        (depending on CFG.target_radunit_optical and target_radunit_thermal).
336
337
        The function can be executed by a L1A_object representing a full scene or a tile. To process a file from disk
        in tiles, provide an item of self.tile_pos as the 'subset' argument."""
Daniel Scheffler's avatar
Daniel Scheffler committed
338

339
        proc_opt_therm = sorted(list(set(self.dict_LayerOptTherm.values())))
340
        assert proc_opt_therm in [['optical', 'thermal'], ['optical'], ['thermal']]
341

342
343
344
        inFill, inZero, inSaturated = \
            self.MetaObj.spec_vals['fill'], self.MetaObj.spec_vals['zero'], self.MetaObj.spec_vals['saturated']
        (rS, rE), (cS, cE) = self.arr_pos or ((0, self.shape_fullArr[0]), (0, self.shape_fullArr[1]))
345
346
        #        in_mem = hasattr(self, 'arr') and isinstance(self.arr, np.ndarray)
        #        if in_mem:
347
348
        #            (rS, rE), (cS, cE) =
        #                self.arr_pos if self.arr_pos else ((0,self.shape_fullArr[0]),(0,self.shape_fullArr[1]))
349
350
351
        #            bands = true_bands = self.arr.shape[2] if len(self.arr.shape) == 3 else 1
        #        else:
        #            subset = subset if subset else ['block', self.arr_pos] if self.arr_pos else ['cube', None]
352
353
        #            bands, rS, rE, cS, cE =
        #                list(GEOP.get_subsetProps_from_subsetArg(self.shape_fullArr, subset).values())[2:7]
354
355
356
357
        #            ds = gdal.Open(self.MetaObj.Dataname); true_bands = ds.RasterCount; ds = None
        assert len(self.LayerBandsAssignment) == self.arr.bands, \
            "DN2RadRef-Input data have %s bands although %s bands are specified in self.LayerBandsAssignment." \
            % (self.arr.bands, len(self.LayerBandsAssignment))
358

359
360
361
362
        ################################
        # perform conversion if needed #
        ################################

363
364
        data_optical, data_thermal, optical_bandsList, thermal_bandsList = None, None, [], []
        for optical_thermal in ['optical', 'thermal']:
365
366
            if optical_thermal not in self.dict_LayerOptTherm.values():
                continue
367
            conv = getattr(CFG, 'target_radunit_%s' % optical_thermal)
368
369
            conv = conv if conv != 'BOA_Ref' else 'TOA_Ref'
            assert conv in ['Rad', 'TOA_Ref', 'Temp'], 'Unsupported conversion type: %s' % conv
370
            arr_desc = self.arr_desc.split('/')[0] if optical_thermal == 'optical' else self.arr_desc.split('/')[-1]
371
            assert arr_desc in ['DN', 'Rad', 'TOA_Ref', 'Temp'], 'Unsupported array description: %s' % arr_desc
372

373
374
            bList = [i for i, lr in enumerate(self.LayerBandsAssignment) if
                     self.dict_LayerOptTherm[lr] == optical_thermal]
375

376
            # custom_subsetProps = [[rS,rE],[cS,cE],bList]
377
378

            inArray = np.take(self.arr, bList, axis=2)
379
            # inArray = np.take(self.arr,bList,axis=2) if in_mem else \
380
            #    INP_R.read_ENVI_image_data_as_array(self.MetaObj.Dataname,'custom',custom_subsetProps,self.logger,q=1)
381
            inArray = inArray[:, :, 0] if len(inArray.shape) == 3 and inArray.shape[2] == 1 else inArray  # 3D to 2D
382

383
384
            if arr_desc == 'DN':
                self.log_for_fullArr_or_firstTile('Calculating %s...' % conv, subset)
385
386
387
388
389
390
391
392
393
394
395

                # get input parameters #
                ########################

                off = [self.MetaObj.Offsets[b] for b in self.LayerBandsAssignment]
                gai = [self.MetaObj.Gains[b] for b in self.LayerBandsAssignment]
                irr = [self.MetaObj.SolIrradiance[b] for b in self.LayerBandsAssignment]
                zen, esd = 90 - float(self.MetaObj.SunElevation), self.MetaObj.EarthSunDist
                k1 = [self.MetaObj.ThermalConstK1[b] for b in self.LayerBandsAssignment]
                k2 = [self.MetaObj.ThermalConstK2[b] for b in self.LayerBandsAssignment]

396
                OFF, GAI, IRR, K1, K2 = [list(np.array(par)[bList]) for par in [off, gai, irr, k1, k2]]
397
398
399

                # perform conversion #
                ######################
400
401
402
403
404
405
                res = \
                    GEOP.DN2Rad(inArray, OFF, GAI, inFill, inZero, inSaturated) if conv == 'Rad' else \
                    GEOP.DN2TOARef(inArray, OFF, GAI, IRR, zen, esd, inFill, inZero,
                                   inSaturated) if conv == 'TOA_Ref' else \
                    GEOP.DN2DegreesCelsius_fastforward(inArray, OFF, GAI, K1, K2, 0.95, inFill, inZero, inSaturated)
                if conv == 'TOA_Ref':
406
                    self.MetaObj.ScaleFactor = CFG.scale_factor_TOARef
407
408

            elif arr_desc == 'Rad':
409
                raise NotImplementedError("Conversion Rad to %s is currently not supported." % conv)
410

411
            elif arr_desc == 'TOA_Ref':
412
                if conv == 'Rad':
413
414
                    """http://s2tbx.telespazio-vega.de/sen2three/html/r2rusage.html?highlight=quantification182
                    rToa = (float)(DN_L1C_band / QUANTIFICATION_VALUE);
415
                    L = (rToa * e0__SOLAR_IRRADIANCE_For_band * cos(Z__Sun_Angles_Grid_Zenith_Values)) /
416
                        (PI * U__earth_sun_distance_correction_factor);
417
418
                    L = (U__earth_sun_distance_correction_factor * rToa * e0__SOLAR_IRRADIANCE_For_band * cos(
                        Z__Sun_Angles_Grid_Zenith_Values)) / PI;"""
419
                    if re.search(r'Sentinel-2', self.satellite, re.I):
420
421
422
                        warnings.warn('Physical gain values unclear for Sentinel-2! This may cause errors when '
                                      'calculating radiance from TOA Reflectance. ESA provides only 12 gain values for '
                                      '13 bands and it not clear for which bands the gains are provided.')
423
424
                    raise NotImplementedError("Conversion TOA_Ref to %s is currently not supported." % conv)
                else:  # conv=='TOA_Ref'
425
426
427
                    if self.MetaObj.ScaleFactor != CFG.scale_factor_TOARef:
                        res = self.rescale_array(inArray, CFG.scale_factor_TOARef, self.MetaObj.ScaleFactor)
                        self.MetaObj.ScaleFactor = CFG.scale_factor_TOARef
428
                        self.log_for_fullArr_or_firstTile(
429
                            'Rescaling Ref data to scaling factor %d.' % CFG.scale_factor_TOARef)
430
431
                    else:
                        res = inArray
432
                        self.log_for_fullArr_or_firstTile('The input data already represents TOA '
433
                                                          'reflectance with the desired scale factor of %d.'
434
                                                          % CFG.scale_factor_TOARef)
435

436
437
            else:  # arr_desc == 'Temp'
                raise NotImplementedError("Conversion Temp to %s is currently not supported." % conv)
438

439
            # ensure int16 as output data type (also relevant for auto-setting of nodata value)
440
            if isinstance(res, (np.ndarray, GeoArray)):
441
                # change data type to int16 and update nodata values within array
442
443
                res = res if res.dtype == np.int16 else res.astype(np.int16)
                res[res == inFill] = get_outFillZeroSaturated(np.int16)[0]
444

445
446
447
448
            if optical_thermal == 'optical':
                data_optical, optical_bandsList = res, bList
            else:
                data_thermal, thermal_bandsList = res, bList
449

450
451
452
453
        #################################################
        # combine results from optical and thermal data #
        #################################################

454
455
456
457
        if data_optical is not None and data_thermal is not None:
            bands_opt, bands_therm = [1 if len(d.shape) == 2 else d.shape[2] for d in [data_optical, data_thermal]]
            r, c, b = data_optical.shape[0], data_optical.shape[1], bands_opt + bands_therm
            dataOut = np.empty((r, c, b), data_optical.dtype)
458

459
460
461
462
463
464
465
            for idx_out, idx_in in zip(optical_bandsList, range(bands_opt)):
                dataOut[:, :, idx_out] = data_optical[:, :, idx_in]
            for idx_out, idx_in in zip(thermal_bandsList, range(bands_therm)):
                dataOut[:, :, idx_out] = data_thermal[:, :, idx_in]
        else:
            dataOut = data_optical if data_thermal is None else data_thermal
        assert dataOut is not None
466

Daniel Scheffler's avatar
Daniel Scheffler committed
467
        self.update_spec_vals_according_to_dtype('int16')
468
        tiles_desc = '_'.join([desc for op_th, desc in zip(['optical', 'thermal'],
469
470
                                                           [CFG.target_radunit_optical,
                                                            CFG.target_radunit_thermal])
471
                               if desc in self.dict_LayerOptTherm.values()])
Daniel Scheffler's avatar
Daniel Scheffler committed
472

473
        self.arr = dataOut
474
475
476
        self.arr_desc = tiles_desc

        return {'desc': tiles_desc, 'row_start': rS, 'row_end': rE, 'col_start': cS, 'col_end': cE, 'data': dataOut}
Daniel Scheffler's avatar
Daniel Scheffler committed
477
478
479
480
481
482

    def validate_GeoTransProj_GeoAlign(self):
        project_dir = os.path.abspath(os.curdir)
        if self.MetaObj.Dataname.startswith('/vsi'):
            os.chdir(os.path.dirname(self.path_archive))
        rasObj = GEOP.GEOPROCESSING(self.MetaObj.Dataname, self.logger)
483
        if rasObj.geotransform == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) and rasObj.projection == '':
484
            if re.search(r'ALOS', self.satellite) and self.MetaObj.ProcLCode == '1B2':
485
                self.GeoTransProj_ok, self.GeoAlign_ok = False, True
486
            else:
487
                self.GeoTransProj_ok, self.GeoAlign_ok = False, False
Daniel Scheffler's avatar
Daniel Scheffler committed
488
        os.chdir(project_dir)
489

Daniel Scheffler's avatar
Daniel Scheffler committed
490
491
492
493
494
    def reference_data(self, out_CS='UTM'):
        """Perform georeferencing of self.arr or the corresponding data on disk respectively.
        Method is skipped if self.GeoAlign_ok and self.GeoTransProj_ok evaluate to 'True'. All attributes connected
        with the georeference of self.arr are automatically updated."""

495
496
        from ..io.output_writer import add_attributes_to_ENVIhdr

Daniel Scheffler's avatar
Daniel Scheffler committed
497
498
        if False in [self.GeoAlign_ok, self.GeoTransProj_ok]:
            previous_dataname = self.MetaObj.Dataname
499
            if hasattr(self, 'arr') and isinstance(self.arr, (GeoArray, np.ndarray)) and \
500
               self.MetaObj.Dataname.startswith('/vsi'):
501
502
503
504
505
506
                outP = os.path.join(self.ExtractedFolder, self.baseN + '__' + self.arr_desc)
                # FIXME ineffective but needed as long as georeference_by_TieP_or_inherent_GCPs does not support
                # FIXME direct array inputs
                GEOP.ndarray2gdal(self.arr, outPath=outP, geotransform=mapinfo2geotransform(self.MetaObj.map_info),
                                  projection=self.MetaObj.projection,
                                  direction=3)
Daniel Scheffler's avatar
Daniel Scheffler committed
507
508
                assert os.path.isfile(outP), 'Writing tempfile failed.'
                self.MetaObj.Dataname = outP
509
            rasObj = GEOP.GEOPROCESSING(self.MetaObj.Dataname, self.logger)
Daniel Scheffler's avatar
Daniel Scheffler committed
510
511
512

            # --Georeference if neccessary
            if self.GeoAlign_ok and not self.GeoTransProj_ok:
513
                self.logger.info('Adding georeference from metadata to image...')
Daniel Scheffler's avatar
Daniel Scheffler committed
514
                rasObj.add_GeoTransform_Projection_using_MetaData(self.MetaObj.CornerTieP_LonLat,
515
516
517
518
519
                                                                  CS_TYPE=self.MetaObj.CS_TYPE,
                                                                  CS_DATUM=self.MetaObj.CS_DATUM,
                                                                  CS_UTM_ZONE=self.MetaObj.CS_UTM_ZONE,
                                                                  gResolution=self.MetaObj.gResolution,
                                                                  shape_fullArr=self.shape_fullArr)
Daniel Scheffler's avatar
Daniel Scheffler committed
520
                self.add_rasterInfo_to_MetaObj(rasObj)
521
                add_attributes_to_ENVIhdr(
522
523
524
                    {'map info': self.MetaObj.map_info, 'coordinate system string': self.MetaObj.projection},
                    os.path.splitext(self.MetaObj.Dataname)[0] + '.hdr')
                self.arr = self.MetaObj.Dataname
Daniel Scheffler's avatar
Daniel Scheffler committed
525
526
527
                self.GeoTransProj_ok = True

            if not self.GeoAlign_ok:
528
                self.logger.info('Georeferencing image...')
Daniel Scheffler's avatar
Daniel Scheffler committed
529
                rasObj.georeference_by_TieP_or_inherent_GCPs(TieP=self.MetaObj.CornerTieP_LonLat, dst_CS=out_CS,
530
531
                                                             dst_CS_datum='WGS84', mode='GDAL', use_workspace=True,
                                                             inFill=self.MetaObj.spec_vals['fill'])
Daniel Scheffler's avatar
Daniel Scheffler committed
532

533
                if not CFG.inmem_serialization:
534
                    path_warped = os.path.join(self.ExtractedFolder, self.baseN + '__' + self.arr_desc)
Daniel Scheffler's avatar
Daniel Scheffler committed
535
536
                    GEOP.ndarray2gdal(rasObj.tondarray(direction=3), path_warped, importFile=rasObj.desc, direction=3)
                    self.MetaObj.Dataname = path_warped
537
                    self.arr = path_warped
538
                else:  # CFG.inmem_serialization is True
Daniel Scheffler's avatar
Daniel Scheffler committed
539
540
                    self.arr = rasObj.tondarray(direction=3)

541
542
                self.shape_fullArr = [rasObj.rows, rasObj.cols, rasObj.bands]
                self.add_rasterInfo_to_MetaObj()  # uses self.MetaObj.Dataname as inFile
Daniel Scheffler's avatar
Daniel Scheffler committed
543
                self.update_spec_vals_according_to_dtype()
544
                self.calc_mask_nodata()  # uses nodata value from self.MetaObj.spec_vals
Daniel Scheffler's avatar
Daniel Scheffler committed
545
546
547
548
549
550
551
552
                self.GeoTransProj_ok, self.GeoAlign_ok = True, True

            if rasObj.get_projection_type() == 'LonLat':
                self.MetaObj.CornerTieP_LonLat = rasObj.get_corner_coordinates('LonLat')

            if rasObj.get_projection_type() == 'UTM':
                self.MetaObj.CornerTieP_UTM = rasObj.get_corner_coordinates('UTM')

553
            if CFG.inmem_serialization:
554
555
                self.delete_tempFiles()  # these files are needed later in Python execution mode
                self.MetaObj.Dataname = previous_dataname  # /vsi.. pointing directly to raw data archive (which exists)
556

557
    def calc_corner_positions(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
558
        """Get true corner positions in the form
559
        [UL, UR, LL, LR] as [(ULrow,ULcol),(URrow,URcol),...],[(ULLon,ULLat),(URLon,URLat),..]"""
560

561
        # set lonlat corner positions for outer image edges
562
563
564
        rows, cols = self.shape_fullArr[:2]
        CorPosXY = [(0, 0), (cols, 0), (0, rows), (cols, rows)]
        gt = self.mask_nodata.geotransform
565
        # using EPSG code ensures that exactly the same WKT code is used in case of in-mem and disk serialization
566
567
        prj = EPSG2WKT(self.mask_nodata.epsg)
        CorLatLon = pixelToLatLon(CorPosXY, geotransform=gt, projection=prj)
568
569
570
        self.corner_lonlat = [tuple(reversed(i)) for i in CorLatLon]

        # set true data corner positions (image coordinates)
571
572
573
574
        assert self.arr_shape == 'cube', 'The given GMS object must represent the full image cube for calculating ,' \
                                         'true corner posistions. Received %s.' % self.arr_shape
        assert hasattr(self,
                       'mask_nodata') and self.mask_nodata is not None, "The L1A object needs to have a nodata mask."
575
        self.logger.info('Calculating true data corner positions (image and world coordinates)...')
576

577
        # if re.search(r'ETM+', self.sensor) and self.acq_datetime > datetime.datetime(year=2003, month=5, day=31,
578
579
        #                                                                             tzinfo=datetime.timezone.utc):
        if is_dataset_provided_as_fullScene(self.GMS_identifier):
580
            kw = dict(algorithm='numpy', assert_four_corners=True)
581
        else:
582
583
            kw = dict(algorithm='shapely', assert_four_corners=False)
        self.trueDataCornerPos = calc_FullDataset_corner_positions(self.mask_nodata, **kw)
584

585
        # set true data corner positions (lon/lat coordinates)
586
587
        trueCorPosXY = [tuple(reversed(i)) for i in self.trueDataCornerPos]
        trueCorLatLon = pixelToLatLon(trueCorPosXY, geotransform=gt, projection=prj)
Daniel Scheffler's avatar
Daniel Scheffler committed
588
        self.trueDataCornerLonLat = [tuple(reversed(i)) for i in trueCorLatLon]
589

590
591
592
593
594
595
        # set true data corner positions (UTM coordinates)
        # calculate trueDataCornerUTM
        if isProjectedOrGeographic(prj) == 'projected':
            data_corners_utmYX = [imYX2mapYX(imYX, gt=gt) for imYX in self.trueDataCornerPos]
            self.trueDataCornerUTM = [(YX[1], YX[0]) for YX in data_corners_utmYX]

596
597
        # set full scene corner positions (image and lonlat coordinates)
        if is_dataset_provided_as_fullScene(self.GMS_identifier):
598
599
600
601
602
            assert len(self.trueDataCornerLonLat) == 4, \
                "Dataset %s (entity ID %s) is provided as full scene. Thus exactly 4 image corners must be present " \
                "within the dataset. Found %s corners instead."\
                % (self.scene_ID, self.entity_ID, len(self.trueDataCornerLonLat))
            self.fullSceneCornerPos = self.trueDataCornerPos
603
            self.fullSceneCornerLonLat = self.trueDataCornerLonLat
604

605
        else:
606

607
            if re.search(r'AVNIR', self.sensor):
608
609
610
611
612
613
                self.fullSceneCornerPos = calc_FullDataset_corner_positions(self.mask_nodata, algorithm='numpy',
                                                                            assert_four_corners=False)
                # set true data corner positions (lon/lat coordinates)
                trueCorPosXY = [tuple(reversed(i)) for i in self.trueDataCornerPos]
                trueCorLatLon = pixelToLatLon(trueCorPosXY, geotransform=gt, projection=prj)
                self.fullSceneCornerLonLat = [tuple(reversed(i)) for i in trueCorLatLon]
614
615
616
617

            else:
                # RapidEye or Sentinel-2 data

618
                if re.search(r'Sentinel-2', self.satellite):
619
                    # get fullScene corner coordinates by database query
620
621
622
                    # -> calculate footprints for all granules of the same S2 datatake
                    # -> merge them and calculate overall corner positions
                    self.fullSceneCornerPos = [tuple(reversed(i)) for i in CorPosXY]  # outer corner positions
623
                    self.fullSceneCornerLonLat = self.corner_lonlat  # outer corner positions
624
625
                else:
                    # RapidEye
626
                    raise NotImplementedError  # FIXME
627

628
    def calc_center_AcqTime(self):
629
        """Calculate scene acquisition time if not given in provider metadata."""
630

631
        if self.MetaObj.AcqTime == '':
632
633
634
635
            self.MetaObj.calc_center_acquisition_time(self.fullSceneCornerLonLat, self.logger)

        # update timestamp
        self.acq_datetime = self.MetaObj.AcqDateTime
636
637

    def calc_orbit_overpassParams(self):
638
        """Calculate orbit parameters."""
639

640
641
642
643
        if is_dataset_provided_as_fullScene(self.GMS_identifier):
            self.MetaObj.overpassDurationSec, self.MetaObj.scene_length = \
                self.MetaObj.get_overpassDuration_SceneLength(
                    self.fullSceneCornerLonLat, self.fullSceneCornerPos, self.shape_fullArr, self.logger)
644
645

    def add_rasterInfo_to_MetaObj(self, custom_rasObj=None):
646
647
648
        """
        Add the attributes 'rows','cols','bands','map_info','projection' and 'physUnit' to self.MetaObj
        """
649
        # TODO is this info still needed in MetaObj?
Daniel Scheffler's avatar
Daniel Scheffler committed
650
        project_dir = os.path.abspath(os.curdir)
651
652
653
        if self.MetaObj.Dataname.startswith('/vsi'):
            os.chdir(os.path.dirname(self.path_archive))

654
        rasObj = custom_rasObj if custom_rasObj else GEOP.GEOPROCESSING(self.MetaObj.Dataname, self.logger)
655
        self.MetaObj.add_rasObj_dims_projection_physUnit(rasObj, self.dict_LayerOptTherm, self.logger)
656

Daniel Scheffler's avatar
Daniel Scheffler committed
657
        prj_type = rasObj.get_projection_type()
658
        assert prj_type, "Projections other than LonLat or UTM are currently not supported. Got. %s." % prj_type
Daniel Scheffler's avatar
Daniel Scheffler committed
659
660
        if prj_type == 'LonLat':
            self.MetaObj.CornerTieP_LonLat = rasObj.get_corner_coordinates('LonLat')
661
662
        else:  # UTM
            self.MetaObj.CornerTieP_UTM = rasObj.get_corner_coordinates('UTM')
663

664
        if self.MetaObj.Dataname.startswith('/vsi'):  # only if everthing is kept in memory
Daniel Scheffler's avatar
Daniel Scheffler committed
665
            os.chdir(project_dir)
666
            self.MetaObj.bands = 1 if len(self.arr.shape) == 2 else self.arr.shape[2]
Daniel Scheffler's avatar
Daniel Scheffler committed
667

668
        self.arr.gt = mapinfo2geotransform(self.MetaObj.map_info)
669
670
671
        if not self.arr.prj:
            self.arr.prj = self.MetaObj.projection

672
673
        # must be set here because nodata mask has been computed from self.arr without geoinfos:
        self.mask_nodata.gt = self.arr.gt
674
        self.mask_nodata.prj = self.arr.prj
675

676
    def update_spec_vals_according_to_dtype(self, dtype=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
677
        """Update self.MetaObj.spec_vals.
678

Daniel Scheffler's avatar
Daniel Scheffler committed
679
        :param dtype:   <str> or <numpy data type> The data type to be used for updating.
680
681
                        If not specified the data type of self.arr is used.
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
682
        if dtype is None:
683
            if hasattr(self, 'arr') and isinstance(self.arr, np.ndarray):
Daniel Scheffler's avatar
Daniel Scheffler committed
684
685
                dtype = self.arr.dtype
            else:
686
687
688
                arr = gdalnumeric.LoadFile(self.arr, 0, 0, 1, 1) if hasattr(self, 'arr') and isinstance(self.arr,
                                                                                                        str) else \
                    gdalnumeric.LoadFile(self.MetaObj.Dataname, 0, 0, 1, 1)
Daniel Scheffler's avatar
Daniel Scheffler committed
689
690
                assert arr is not None
                dtype = arr.dtype
691

692
693
        self.MetaObj.spec_vals['fill'], self.MetaObj.spec_vals['zero'], self.MetaObj.spec_vals['saturated'] = \
            get_outFillZeroSaturated(dtype)
694
        self.arr.nodata = self.MetaObj.spec_vals['fill']
695
696

    def calc_mean_VAA(self):
697
        """Calculates mean viewing azimuth angle using sensor flight line derived from full scene corner coordinates."""
698

699
700
701
702
703
        if is_dataset_provided_as_fullScene(self.GMS_identifier):
            self.VAA_mean = \
                GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
        else:
            # e.g. Sentinel-2 / RapidEye
704
705
706
707
            self.logger.debug('No precise calculation of mean viewing azimuth angle possible because orbit track '
                              'cannot be reconstructed from dataset since full scene corner positions are unknown. '
                              'Mean VAA angle is filled with the mean value of the viewing azimuth array provided '
                              'in metadata.')
708
            self.VAA_mean = self.MetaObj.IncidenceAngle
709

710
        self.logger.info('Calculation of mean VAA...: %s' % round(self.VAA_mean, 2))