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