metadata.py 128 KB
Newer Older
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
"""Module 'metadata' for handling any type of metadata of GeoMultiSens compatible sensor systems."""
Daniel Scheffler's avatar
Daniel Scheffler committed
28 29

from __future__ import (division, print_function, unicode_literals, absolute_import)
30

31 32
import collections
import datetime
33
import glob
34 35 36 37
import math
import os
import re
import sys
38
import warnings
39
import iso8601
40
import xml.etree.ElementTree as ET
41
from typing import List, TYPE_CHECKING  # noqa F401  # flake8 issue
42 43 44

import numpy as np
import pyproj
45 46
from matplotlib import dates as mdates
from pyorbital import astronomy
47
from natsort import natsorted
48

Daniel Scheffler's avatar
Daniel Scheffler committed
49 50
from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.geo.projection import WKT2EPSG
Daniel Scheffler's avatar
Daniel Scheffler committed
51
from pyrsr import RSR
52
from sicor.options import get_options as get_ac_options
53

54
from ..options.config import GMS_config as CFG
55
from ..io.input_reader import open_specific_file_within_archive, Solar_Irradiance_reader, SRF_Reader
56 57 58 59 60 61 62
from ..io.output_writer import enviHdr_keyOrder
from ..algorithms import geoprocessing as GEOP
from ..misc import helper_functions as HLP_F
from ..misc import database_tools as DB_T
from ..misc.path_generator import path_generator, get_path_ac_options
from ..misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene, datasetid_to_sat_sen
from ..misc.exceptions import ACNotSupportedError
63

64 65
if TYPE_CHECKING:
    from ..model.gms_object import GMS_identifier  # noqa F401  # flake8 issue
66

67 68
__author__ = 'Daniel Scheffler', 'Robert Behling'

69

70
class METADATA(object):
71
    def __init__(self, GMS_id):
72 73 74
        # private attributes
        self._AcqDateTime = None

75
        # unpack GMS_identifier
76 77 78 79 80 81 82
        self.GMS_identifier = GMS_id
        self.image_type = GMS_id.image_type
        self.Satellite = GMS_id.satellite
        self.Sensor = GMS_id.sensor
        self.Subsystem = GMS_id.subsystem
        self.proc_level = GMS_id.proc_level
        self.logger = GMS_id.logger
83 84 85

        self.Dataname = ''
        self.FolderOrArchive = ''
Daniel Scheffler's avatar
Daniel Scheffler committed
86 87 88
        self.Metafile = ""  # File containing image metadata (automatically found)
        self.EntityID = ""  # ID to identify the original scene
        self.SceneID = ''  # postgreSQL-database identifier
89
        self.Sensormode = ""
Daniel Scheffler's avatar
Daniel Scheffler committed
90 91 92
        self.gResolution = -99.  # resolution [m]
        self.AcqDate = ""  # YYYY-MM-DD
        self.AcqTime = ""  # HH:MM:SS
93 94 95 96 97 98
        self.Offsets = {}  # Dict with offset for each band (radiance)
        self.Gains = {}  # Dict with gain for each band (radiance)
        self.OffsetsRef = {}  # Dict with offset for each band for conversion to Reflectance (Landsat8)
        self.GainsRef = {}  # Dict with gain for each band for conversion to Reflectance (Landsat8)
        self.CWL = {}
        self.FWHM = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
99 100 101 102 103
        self.ProcLCode = ""  # processing Level: Sensor specific Code
        self.ProcLStr = ""  # processing Level: Sensor independent String (raw,rad cal, rad+geom cal, ortho)
        self.SunElevation = -99.0  # Sun elevation angle at center of product [Deg]
        # Sun azimuth angle at center of product, in degrees from North (clockwise) at the time of the first image line
        self.SunAzimuth = -99.0
104
        self.SolIrradiance = []
105 106
        self.ThermalConstK1 = {}
        self.ThermalConstK2 = {}
107
        self.EarthSunDist = 1.0
Daniel Scheffler's avatar
Daniel Scheffler committed
108 109 110
        # viewing zenith angle of the Sensor (offNadir) [Deg] (usually) in Case of RapidEye "+"
        # being East and “-” being West
        self.ViewingAngle = -99.0
111
        self.ViewingAngle_arrProv = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
112 113 114
        # viewing azimuth angle. The angle between the view direction of the satellite and a line perpendicular to
        # the image or tile center.[Deg]
        self.IncidenceAngle = -9999.99
115
        self.IncidenceAngle_arrProv = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
116 117 118
        self.FOV = 9999.99  # field of view of the sensor [Deg]
        # Sensor specific quality code: See ro2/behling/Satelliten/doc_allg/ReadOutMetadata/SatMetadata.xls
        self.Quality = []
119
        self.PhysUnit = "DN"
Daniel Scheffler's avatar
Daniel Scheffler committed
120 121 122 123
        # Factor to get reflectance values in [0-1]. Sentinel2A provides scaling factors for the Level1C
        # TOA-reflectance products
        self.ScaleFactor = -99
        self.CS_EPSG = -99.  # EPSG-Code of coordinate system
124 125 126 127 128
        self.CS_TYPE = ""
        self.CS_DATUM = ""
        self.CS_UTM_ZONE = -99.
        self.WRS_path = -99.
        self.WRS_row = -99.
Daniel Scheffler's avatar
Daniel Scheffler committed
129 130
        # List of Corner Coordinates in order of Lon/Lat/DATA_X/Data_Y for all given image corners
        self.CornerTieP_LonLat = []
131
        self.CornerTieP_UTM = []
Daniel Scheffler's avatar
Daniel Scheffler committed
132
        self.LayerBandsAssignment = []  # List of spectral bands in the same order as they are stored in the rasterfile.
133 134 135 136
        self.additional = []
        self.image_type = 'RSD'
        self.orbitParams = {}
        self.overpassDurationSec = -99.
Daniel Scheffler's avatar
Daniel Scheffler committed
137
        self.scene_length = -99.
138 139 140
        self.rows = -99.
        self.cols = -99.
        self.bands = -99.
141
        self.nBands = -99.
142 143 144
        self.map_info = []
        self.projection = ""
        self.wvlUnit = ""
Daniel Scheffler's avatar
Daniel Scheffler committed
145
        self.spec_vals = {'fill': None, 'zero': None, 'saturated': None}
146

147 148 149
        self.version_gms_preprocessing = CFG.version
        self.versionalias_gms_preprocessing = CFG.versionalias

150 151 152 153 154 155 156 157 158 159
    def read_meta(self, scene_ID, stacked_image, data_folderOrArchive, LayerBandsAssignment=None):
        """
        Read metadata.
        """

        self.SceneID = scene_ID
        self.Dataname = stacked_image
        self.FolderOrArchive = data_folderOrArchive
        self.LayerBandsAssignment = LayerBandsAssignment if LayerBandsAssignment else []

160
        if re.search(r"SPOT", self.Satellite, re.I):
161
            self.Read_SPOT_dimap2()
162
        elif re.search(r"Terra", self.Satellite, re.I):
163
            self.Read_ASTER_hdffile(self.Subsystem)
164
        elif re.search(r"Sentinel-2A", self.Satellite, re.I) or re.search(r"Sentinel-2B", self.Satellite, re.I):
165
            self.Read_Sentinel2_xmls()
166
        elif re.search(r"LANDSAT", self.Satellite, re.I):
167
            self.Read_LANDSAT_mtltxt(self.LayerBandsAssignment)
168
        elif re.search(r"RapidEye", self.Satellite, re.I):
169
            self.Read_RE_metaxml()
170
        elif re.search(r"ALOS", self.Satellite, re.I):
171 172 173
            self.Read_ALOS_summary()
            self.Read_ALOS_LEADER()  # for gains and offsets
        else:
174
            raise NotImplementedError("%s is not a supported sensor or sensorname is misspelled." % self.Satellite)
175 176

        return self
177

178 179 180 181 182 183 184
    def __getstate__(self):
        """Defines how the attributes of MetaObj instances are pickled."""
        if self.logger:
            self.logger.close()

        return self.__dict__

185 186 187
    @property
    def AcqDateTime(self):
        """Returns a datetime.datetime object containing date, time and timezone (UTC time)."""
188 189 190
        if not self._AcqDateTime and self.AcqDate and self.AcqTime:
            self._AcqDateTime = datetime.datetime.strptime('%s %s%s' % (self.AcqDate, self.AcqTime, '.000000+0000'),
                                                           '%Y-%m-%d %H:%M:%S.%f%z')
Daniel Scheffler's avatar
Daniel Scheffler committed
191

192 193 194 195 196
        return self._AcqDateTime

    @AcqDateTime.setter
    def AcqDateTime(self, DateTime):
        # type: (datetime.datetime) -> None
197 198 199 200 201 202 203
        if isinstance(DateTime, str):
            self._AcqDateTime = datetime.datetime.strptime(DateTime, '%Y-%m-%d %H:%M:%S.%f%z')
        elif isinstance(DateTime, datetime.datetime):
            self._AcqDateTime = DateTime

        self.AcqDate = DateTime.strftime('%Y-%m-%d')
        self.AcqTime = DateTime.strftime('%H:%M:%S')
204

205 206 207
    @property
    def overview(self):
        attr2include = \
Daniel Scheffler's avatar
Daniel Scheffler committed
208
            ['Dataname', 'Metafile', 'EntityID', 'Satellite', 'Sensor', 'Subsystem', 'gResolution', 'AcqDate',
209
             'AcqTime', 'CWL', 'FWHM', 'Offsets', 'Gains', 'ProcLCode', 'ProcLStr', 'SunElevation', 'SunAzimuth',
210
             'ViewingAngle', 'IncidenceAngle', 'FOV', 'SolIrradiance', 'ThermalConstK1', 'ThermalConstK2',
Daniel Scheffler's avatar
Daniel Scheffler committed
211
             'EarthSunDist', 'Quality', 'PhysUnit', 'additional', 'GainsRef', 'OffsetsRef', 'CornerTieP_LonLat',
212
             'CS_EPSG', 'CS_TYPE', 'CS_DATUM', 'CS_UTM_ZONE', 'LayerBandsAssignment']
213 214
        return collections.OrderedDict((attr, getattr(self, attr)) for attr in attr2include)

215 216 217 218 219 220 221 222 223
    @property
    def LayerBandsAssignment_full(self):
        # type: () -> list
        """Return complete LayerBandsAssignment without excluding thermal or panchromatic bands.

        NOTE: CFG.sort_bands_by_cwl is respected, so returned list may be sorted by central wavelength
        """
        return get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False, return_fullLBA=True)

224 225 226 227
    @property
    def bandnames(self):
        return [('Band %s' % i) for i in self.LayerBandsAssignment]

228 229 230 231
    def Read_SPOT_dimap2(self):
        """----METHOD_2------------------------------------------------------------
        # read metadata from spot dimap file
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
232

Daniel Scheffler's avatar
Daniel Scheffler committed
233
        # self.default_attr()
234
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
235 236
            glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/scene01/metadata.dim'))
            assert len(glob_res) > 0, 'No metadata.dim file can be found in %s!' % self.FolderOrArchive
237 238
            self.Metafile = glob_res[0]
            dim_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
239
        else:  # archive
240
            dim_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/scene01/metadata.dim')
241

Daniel Scheffler's avatar
Daniel Scheffler committed
242
        # special values
243 244
        h1 = re.findall(r"<SPECIAL_VALUE_INDEX>([a-zA-Z0-9]*)</SPECIAL_VALUE_INDEX>", dim_, re.I)
        h2 = re.findall(r"<SPECIAL_VALUE_TEXT>([a-zA-Z0-9]*)</SPECIAL_VALUE_TEXT>", dim_, re.I)
245
        self.additional.append(["SpecialValues:"])
Daniel Scheffler's avatar
Daniel Scheffler committed
246 247
        for ii, ind in enumerate(h1):
            self.additional[0].append(["%s:%s" % (ind, h2[ii])])
248

Daniel Scheffler's avatar
Daniel Scheffler committed
249
        # EntityID
250
        h3 = re.search(r"<SOURCE_ID>([a-zA-Z0-9]*)</SOURCE_ID>", dim_, re.I)
251 252
        self.EntityID = h3.group(1)

Daniel Scheffler's avatar
Daniel Scheffler committed
253
        # AcqDate
254
        h4 = re.search(r"<IMAGING_DATE>([0-9-]*)</IMAGING_DATE>", dim_, re.I)
255
        AcqDate = h4.group(1)
256

Daniel Scheffler's avatar
Daniel Scheffler committed
257
        # Earth sun distance
258
        self.EarthSunDist = self.get_EarthSunDistance(AcqDate)
259

Daniel Scheffler's avatar
Daniel Scheffler committed
260
        # AcqTime
261
        h5 = re.search(r"<IMAGING_TIME>([0-9:]*)</IMAGING_TIME>", dim_, re.I)
262 263
        AcqTime = h5.group(1)

Daniel Scheffler's avatar
Daniel Scheffler committed
264
        # set self.AcqDateTime as well as self.AcqDate and self.AcqTime
265
        self.AcqDateTime = datetime.datetime.strptime(
Daniel Scheffler's avatar
Daniel Scheffler committed
266
            '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z')
267

Daniel Scheffler's avatar
Daniel Scheffler committed
268
        # Satellite + Sensor
269 270 271 272 273 274
        h6 = re.search(r"<MISSION>([a-zA-Z]*)</MISSION>[a-zA-Z0-9\s]*"
                       r"<MISSION_INDEX>([0-9]*)</MISSION_INDEX>[a-zA-Z0-9\s]*"
                       r"<INSTRUMENT>([a-zA-Z]*)</INSTRUMENT>[a-zA-Z0-9\s]*"
                       r"<INSTRUMENT_INDEX>([0-9]*)</INSTRUMENT_INDEX>[a-zA-Z0-9\s]*"
                       r"<SENSOR_CODE>([a-zA-Z0-9]*)</SENSOR_CODE>",
                       dim_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
275 276
        self.Satellite = "-".join([h6.group(1), h6.group(2)])
        self.Sensor = "".join([h6.group(3), h6.group(4), h6.group(5)])
277

Daniel Scheffler's avatar
Daniel Scheffler committed
278
        # Angles: incidence angle, sunAzimuth, sunElevation
279 280 281
        h7 = re.search(r"<INCIDENCE_ANGLE>(.*)</INCIDENCE_ANGLE>[\s\S]*"
                       r"<SUN_AZIMUTH>(.*)</SUN_AZIMUTH>[\s\S]"
                       r"*<SUN_ELEVATION>(.*)</SUN_ELEVATION>", dim_, re.I)
282 283 284 285
        self.IncidenceAngle = float(h7.group(1))
        self.SunAzimuth = float(h7.group(2))
        self.SunElevation = float(h7.group(3))

Daniel Scheffler's avatar
Daniel Scheffler committed
286
        # Viewing Angle: not always included in the metadata.dim file
287
        h8 = re.search(r"<VIEWING_ANGLE>(.*)</VIEWING_ANGLE>", dim_, re.I)
288
        if type(h8).__name__ == 'NoneType':
Daniel Scheffler's avatar
Daniel Scheffler committed
289
            self.ViewingAngle = 0
290
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
291
            self.ViewingAngle = float(h8.group(1))
292

Daniel Scheffler's avatar
Daniel Scheffler committed
293
        # Field of View
294
        self.FOV = get_FieldOfView(self.GMS_identifier)
295

Daniel Scheffler's avatar
Daniel Scheffler committed
296
        # Units
297
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
298
        self.PhysUnit = "DN"
299

300
        # ProcLevel
301
        h11 = re.search(r"<PROCESSING_LEVEL>([a-zA-Z0-9]*)</PROCESSING_LEVEL>", dim_, re.I)
302 303 304
        self.ProcLCode = h11.group(1)

        # Quality
305 306 307 308 309
        h12 = re.findall(r"<Bad_Pixel>[\s]*"
                         r"<PIXEL_INDEX>([0-9]*)</PIXEL_INDEX>[\s]*"
                         r"<BAD_PIXEL_STATUS>([^<]*)</BAD_PIXEL_STATUS>"
                         r"</Bad_Pixel>", dim_,
                         re.I)
310 311 312
        self.Quality = h12

        # Coordinate Reference System
313 314 315
        re_CS_EPSG = re.search(r'<HORIZONTAL_CS_CODE>epsg:([0-9]*)</HORIZONTAL_CS_CODE>', dim_, re.I)
        re_CS_TYPE = re.search(r'<HORIZONTAL_CS_TYPE>([a-zA-Z0-9]*)</HORIZONTAL_CS_TYPE>', dim_, re.I)
        re_CS_DATUM = re.search(r'<HORIZONTAL_CS_NAME>([\w+\s]*)</HORIZONTAL_CS_NAME>', dim_, re.I)
316 317 318 319 320 321
        self.CS_EPSG = int(re_CS_EPSG.group(1)) if re_CS_EPSG is not None else self.CS_EPSG
        self.CS_TYPE = 'LonLat' if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'GEOGRAPHIC' else 'UTM' \
            if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'UTM' else self.CS_TYPE
        self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS 84' else self.CS_DATUM

        # Corner Coordinates
322 323
        h121 = re.findall(r"<TIE_POINT_CRS_X>(.*)</TIE_POINT_CRS_X>", dim_, re.I)
        h122 = re.findall(r"<TIE_POINT_CRS_Y>(.*)</TIE_POINT_CRS_Y>", dim_, re.I)
324 325 326 327 328 329 330 331 332 333 334 335 336
        if len(h121) == 4 and self.CS_TYPE == 'LonLat':
            # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR)
            self.CornerTieP_LonLat.append(tuple([float(h121[0]), float(h122[0])]))  # UL
            self.CornerTieP_LonLat.append(tuple([float(h121[1]), float(h122[1])]))  # UR
            self.CornerTieP_LonLat.append(tuple([float(h121[3]), float(h122[3])]))  # LL
            self.CornerTieP_LonLat.append(tuple([float(h121[2]), float(h122[2])]))  # LR
            #    ul_lon,ul_lat = self.inDs.GetGCPs()[0].GCPX,self.inDs.GetGCPs()[0].GCPY # funktioniert nur bei SPOT
            #    lr_lon,lr_lat = self.inDs.GetGCPs()[2].GCPX,self.inDs.GetGCPs()[2].GCPY

        ##########################
        # band specific metadata #
        ##########################

337
        LBA_full_sorted = natsorted(self.LayerBandsAssignment_full)
338

Daniel Scheffler's avatar
Daniel Scheffler committed
339
        # Gains and Offsets
340
        h9 = re.search(r"<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I)
341
        h9_ = h9.group(0)
342 343 344
        h91 = re.findall(r"<PHYSICAL_UNIT>([^<]*)</PHYSICAL_UNIT>", h9_, re.I)
        h92 = re.findall(r"<PHYSICAL_BIAS>([^<]*)</PHYSICAL_BIAS>", h9_, re.I)
        h93 = re.findall(r"<PHYSICAL_GAIN>([^<]*)</PHYSICAL_GAIN>", h9_, re.I)
345
        self.additional.append(["Physical Units per band:"])
346
        for ii, ind in enumerate(h91):  # FIXME does not respect sort_bands_by_cwl
347
            self.additional[-1].append(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
348
        # Offsets
349 350
        for ii, (ind, band) in enumerate(zip(h92, LBA_full_sorted)):
            self.Offsets[band] = float(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
351
        # Gains
352
        for ii, (ind, band) in enumerate(zip(h93, LBA_full_sorted)):
Daniel Scheffler's avatar
Daniel Scheffler committed
353
            # gains in dimap file represent reciprocal 1/gain
354
            self.Gains[band] = 1 / float(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
355 356 357 358

        # Solar irradiance, central wavelengths , full width half maximum per band
        self.wvlUnit = 'Nanometers'
        # derive number of bands from number of given gains/offsets in header file or from raster data
359
        # noinspection PyBroadException
Daniel Scheffler's avatar
Daniel Scheffler committed
360 361 362 363
        try:
            self.nBands = (np.mean([len(self.Gains), len(self.Offsets)])
                           if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else
                           GEOP.GEOPROCESSING(self.Dataname, self.logger).bands)
364
        except Exception:
Daniel Scheffler's avatar
Daniel Scheffler committed
365 366
            self.logger.warning('Unable to get number of bands for the dataset %s Provider values are used for '
                                'solar irradiation, CWL and FWHM!.' % self.Dataname)
367
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)
368

Daniel Scheffler's avatar
Daniel Scheffler committed
369
        # Exact values calculated based in SRFs
370
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
371
        # Provider values
372
        if not self.SolIrradiance:
373 374
            h10 = re.search(r"<Solar_Irradiance>[\s\S]*"
                            r"</Solar_Irradiance>", dim_, re.I)
375
            h10_ = h10.group(0)
376 377
            h101 = re.findall(r"<SOLAR_IRRADIANCE_VALUE>([^<]*)"
                              r"</SOLAR_IRRADIANCE_VALUE>", h10_, re.I)
378
            if h101:
379
                self.SolIrradiance = dict(zip(LBA_full_sorted, h101))
Daniel Scheffler's avatar
Daniel Scheffler committed
380 381 382 383
                #    self.additional.append(["Solar Irradiance per band:"])
                #    for ii,ind in enumerate(h101):
                #       self.additional[-1].append(ind)
            else:  # Preconfigured Irradiation values
384
                spot_irr_dic = {
385 386 387 388 389 390 391 392 393 394
                    'SPOT1a': dict(zip(LBA_full_sorted, [1855., 1615., 1090., 1680.])),
                    'SPOT1b': dict(zip(LBA_full_sorted, [1845., 1575., 1040., 1690.])),
                    'SPOT2a': dict(zip(LBA_full_sorted, [1865., 1620., 1085., 1705.])),
                    'SPOT2b': dict(zip(LBA_full_sorted, [1865., 1615., 1090., 1670.])),
                    'SPOT3a': dict(zip(LBA_full_sorted, [1854., 1580., 1065., 1668.])),
                    'SPOT3b': dict(zip(LBA_full_sorted, [1855., 1597., 1067., 1667.])),
                    'SPOT4a': dict(zip(LBA_full_sorted, [1843., 1568., 1052., 233., 1568.])),
                    'SPOT4b': dict(zip(LBA_full_sorted, [1851., 1586., 1054., 240., 1586.])),
                    'SPOT5a': dict(zip(LBA_full_sorted, [1858., 1573., 1043., 236., 1762.])),
                    'SPOT5b': dict(zip(LBA_full_sorted, [1858., 1575., 1047., 234., 1773.]))}
395
                self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)]
Daniel Scheffler's avatar
Daniel Scheffler committed
396
            # Preconfigured CWLs --   # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined
397 398
            # from SPOT-HRV and Landsat-TM Data'; Hill 1990 Comparative Analysis of Landsat-5 TM and SPOT HRV-1 Data for
            # Use in Multiple Sensor Approaches ; # resource: SPOT techical report: Resolutions and spectral modes
399 400 401 402 403 404 405
            sensorcode = get_GMS_sensorcode(self.GMS_identifier)[:2]
            if sensorcode in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b', 'SPOT3a', 'SPOT3b']:
                self.CWL = dict(zip(LBA_full_sorted, [545., 638., 819., 615.]))
            elif sensorcode in ['SPOT4a', 'SPOT4b']:
                self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 645.]))
            elif sensorcode == ['SPOT5a', 'SPOT5b']:
                self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 595.]))
406

407 408 409 410
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

411 412 413 414 415
    def Read_LANDSAT_mtltxt(self, LayerBandsAssignment):
        """----METHOD_3------------------------------------------------------------
        read metadata from LANDSAT metafile: <dataname>.MTL.txt. Metadatafile of LPGS processing chain
        :param LayerBandsAssignment:
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
416

417 418
        # self.default_attr()
        self.LayerBandsAssignment = LayerBandsAssignment
Daniel Scheffler's avatar
Daniel Scheffler committed
419
        self.nBands = len(LayerBandsAssignment)
420
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
421 422
            glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*MTL.txt'))
            assert len(glob_res) > 0, 'No *.MTL metadata file can be found in %s!' % self.FolderOrArchive
423 424
            self.Metafile = glob_res[0]
            mtl_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
425
        else:  # archive
426
            mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt')
427

428
        # Processing Level
429
        h1 = re.search(r'PRODUCT_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
430
        if h1 is None:
431
            h1 = re.search(r'DATA_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
432 433 434
        self.ProcLCode = h1.group(1)

        # Satellite + Sensor + Sensor Mode
435 436 437 438
        h2 = re.search(r'SPACECRAFT_ID = "([a-zA-Z0-9_]*)"[\s]*'
                       r'SENSOR_ID = "([a-zA-Z0-9+]*)"[\s]*'
                       r'SENSOR_MODE = "([\S]*)"',
                       mtl_, re.I)
439
        if h2:
440
            self.Satellite = 'Landsat-%s' % re.search(r'LANDSAT[\D]*([0-9])', h2.group(1), re.I).group(1)
Daniel Scheffler's avatar
Daniel Scheffler committed
441
            self.Sensor = h2.group(2)
442 443
            self.Sensormode = h2.group(3)
        else:
444 445 446 447
            h2a = re.search(r'SPACECRAFT_ID = "([a-zA-Z0-9_]*)"', mtl_, re.I)
            h2b = re.search(r'SENSOR_ID = "([a-zA-Z0-9_+]*)"', mtl_, re.I)
            h2c = re.search(r'SENSOR_MODE = "([a-zA-Z0-9_+]*)"', mtl_, re.I)
            self.Satellite = 'Landsat-%s' % re.search(r'LANDSAT[\D]*([0-9])', h2a.group(1), re.I).group(1)
Daniel Scheffler's avatar
Daniel Scheffler committed
448 449 450
            self.Sensor = h2b.group(1)
            self.Sensormode = h2c.group(1) if h2c is not None else self.Sensormode  # Landsat-8
        self.Sensor = 'ETM+' if self.Sensor == 'ETM' else self.Sensor
451 452

        # EntityID
453
        h2c = re.search(r'LANDSAT_SCENE_ID = "([A-Z0-9]*)"', mtl_, re.I)
454
        if h2c:
455 456
            self.EntityID = h2c.group(1)

457
        # Acquisition Date + Time
458 459 460
        h3 = re.search(r'ACQUISITION_DATE = ([0-9-]*)[\s]*'
                       r'SCENE_CENTER_SCAN_TIME = "?([0-9:]*)"?',
                       mtl_, re.I)
461
        if h3 is None:
462 463 464
            h3 = re.search(r'DATE_ACQUIRED = ([0-9-]*)[\s]*'
                           r'SCENE_CENTER_TIME = "?([0-9:]*)"?',
                           mtl_, re.I)
465 466 467 468 469 470
        AcqDate = h3.group(1)
        AcqTime = h3.group(2)

        # set self.AcqDateTime as well as self.AcqDate and self.AcqTime
        self.AcqDateTime = datetime.datetime.strptime(
            '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z')
471 472 473 474

        # Earth sun distance
        self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate)

475 476
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
477
        self.PhysUnit = "DN"
478

479
        # Angles: incidence angle, sunAzimuth, sunElevation, field of view
480 481 482
        h5 = re.search(r"SUN_AZIMUTH = ([\S]*)[\s]*"
                       r"SUN_ELEVATION = ([\S]*)",
                       mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
483
        self.SunAzimuth = float(h5.group(1))
484
        self.SunElevation = float(h5.group(2))
Daniel Scheffler's avatar
Daniel Scheffler committed
485
        self.FOV = get_FieldOfView(self.GMS_identifier)
486 487

        # Quality
488 489 490
        h6 = re.search(r"GROUP = CORRECTIONS_APPLIED[\s\S]*"
                       r"END_GROUP = CORRECTIONS_APPLIED",
                       mtl_, re.I)
491
        if h6 is None:
492 493 494
            h6 = re.search(r"GROUP = IMAGE_ATTRIBUTES[\s\S]*"
                           r"END_GROUP = IMAGE_ATTRIBUTES",
                           mtl_, re.I)
495

Daniel Scheffler's avatar
Daniel Scheffler committed
496
        h6_ = h6.group(0)
497
        h61 = (h6_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
498
        x = 0
499
        for i in h61:
Daniel Scheffler's avatar
Daniel Scheffler committed
500
            if x == 0 or x + 1 == len(h61):
501 502
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
503
                i_ = i.strip().replace('"', "")
504
                self.Quality.append(i_.split(" = "))
Daniel Scheffler's avatar
Daniel Scheffler committed
505
            x += 1
506

507
        # Additional: coordinate system, geom. Resolution
508
        h7 = re.search(r"GROUP = PROJECTION_PARAMETERS[\s\S]*END_GROUP = L1_METADATA_FILE", mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
509
        h7_ = h7.group(0)
510
        h71 = (h7_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
511
        for x, i in enumerate(h71):
512
            if re.search(r"Group", i, re.I):
513 514
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
515
                i_ = i.strip().replace('"', "")
516
                self.additional.append(i_.split(" = "))
517 518 519 520
        re_CS_TYPE = re.search(r'MAP_PROJECTION = "([\w+\s]*)"', h7_, re.I)
        re_CS_DATUM = re.search(r'DATUM = "([\w+\s]*)"', h7_, re.I)
        re_CS_UTM_ZONE = re.search(r'ZONE_NUMBER = ([0-9]*)\n', mtl_, re.I)
        re_CS_UTM_ZONE = re.search(r'UTM_ZONE = ([0-9]*)\n', h7_,
Daniel Scheffler's avatar
Daniel Scheffler committed
521 522 523 524
                                   re.I) if re_CS_UTM_ZONE is None else re_CS_UTM_ZONE  # Landsat8
        self.CS_TYPE = 'LonLat' if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'GEOGRAPHIC' else 'UTM' \
            if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'UTM' else self.CS_TYPE
        self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS84' else self.CS_DATUM
525 526 527
        self.CS_UTM_ZONE = int(re_CS_UTM_ZONE.group(1)) if re_CS_UTM_ZONE is not None else self.CS_UTM_ZONE
        # viewing Angle
        self.ViewingAngle = 0
Daniel Scheffler's avatar
Daniel Scheffler committed
528
        # Landsat8
529
        h8 = re.search(r"ROLL_ANGLE = ([\S]*)", mtl_, re.I)
530
        if h8:
531 532 533
            self.ViewingAngle = float(h8.group(1))

        # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples)
534 535 536 537
        h10_UL = re.findall(r"PRODUCT_UL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_UR = re.findall(r"PRODUCT_UR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_LL = re.findall(r"PRODUCT_LL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_LR = re.findall(r"PRODUCT_LR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
538
        if not h10_UL:  # Landsat8
539 540 541 542
            h10_UL = re.findall(r"CORNER_UL_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_UR = re.findall(r"CORNER_UR_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_LL = re.findall(r"CORNER_LL_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_LR = re.findall(r"CORNER_LR_[\S]+ = (.*)[\S]*", mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
543 544 545 546 547 548 549 550 551
        if h10_UL:  # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR)
            self.CornerTieP_LonLat.append(tuple([float(h10_UL[i]) for i in [1, 0]]))
            self.CornerTieP_LonLat.append(tuple([float(h10_UR[i]) for i in [1, 0]]))
            self.CornerTieP_LonLat.append(tuple([float(h10_LL[i]) for i in [1, 0]]))
            self.CornerTieP_LonLat.append(tuple([float(h10_LR[i]) for i in [1, 0]]))
            self.CornerTieP_UTM.append(tuple([float(h10_UL[i]) for i in [2, 3]]))
            self.CornerTieP_UTM.append(tuple([float(h10_UR[i]) for i in [2, 3]]))
            self.CornerTieP_UTM.append(tuple([float(h10_LL[i]) for i in [2, 3]]))
            self.CornerTieP_UTM.append(tuple([float(h10_LR[i]) for i in [2, 3]]))
552 553

        # WRS path/row
554
        h11_p = re.search(r'WRS_PATH = ([0-9]*)', mtl_, re.I)
555
        if h11_p:
556
            self.WRS_path = h11_p.group(1)
557 558
        h11_r1 = re.search(r'STARTING_ROW = ([0-9]*)', mtl_, re.I)
        h11_r2 = re.search(r'ENDING_ROW = ([0-9]*)', mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
559
        if h11_r1 is None:  # Landsat-8
560
            h11_r = re.search(r'WRS_ROW = ([0-9]*)', mtl_, re.I)
561
            self.WRS_row = int(h11_r.group(1))
562
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
563
            self.WRS_row = int(h11_r1.group(1)) if h11_r1.group(1) == h11_r2.group(1) else self.WRS_row
564 565

        # Fill missing values
Daniel Scheffler's avatar
Daniel Scheffler committed
566 567 568
        if self.EntityID == '':
            self.logger.info('Scene-ID could not be extracted and has to be retrieved from %s metadata database...'
                             % self.Satellite)
569 570
            result = DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['entityID'],
                                                     {'id': self.SceneID})
571

Daniel Scheffler's avatar
Daniel Scheffler committed
572
            if len(result) == 1:  # e.g. [('LE71950282003121EDC00',)]
573 574
                self.EntityID = result[0][0]
            elif len(result) == 0:
575 576
                self.logger.warning("Scene ID could not be assigned because no dataset matching to the query "
                                    "parameters could be found in metadata database.")
Daniel Scheffler's avatar
Daniel Scheffler committed
577
            else:  # e.g. [('LE71950282003121EDC00',), ('LE71950282003105ASN00',)]
578
                self.logger.warning("Scene ID could not be assigned because %s datasets matching to the query "
Daniel Scheffler's avatar
Daniel Scheffler committed
579
                                    "parameters were found in metadata database." % len(result))
580
        # if self.EntityID=='':
Daniel Scheffler's avatar
Daniel Scheffler committed
581 582
        #     dataset = 'LANDSAT_TM' if self.Satellite=='Landsat-5' else \
        #          'LANDSAT_ETM' if self.Satellite=='Landsat-7' else 'LANDSAT_8' if self.Satellite=='Landsat-8' else ''
583 584 585 586 587 588 589 590 591
        #     if dataset != '':
        #         webmeta = list(usgsapi.search(dataset,'EE',distance=0,ll={'longitude':self.CornerTieP_LonLat[2][0], \
        #                     'latitude':self.CornerTieP_LonLat[2][1]},ur={'longitude':self.CornerTieP_LonLat[1][0], \
        #                     'latitude':self.CornerTieP_LonLat[1][1]},start_date=self.AcqDate,end_date=self.AcqDate))
        #         if len(webmeta)==1:
        #             self.EntityID=webmeta[0]['entityId']
        #         else:
        #             sen  = {'MSS':'M','TM':'T','ETM+':'E','OLI_TIRS':'C','OLI':'O'}[self.Sensor]
        #             sat  = self.Satellite.split('-')[1]
592
        #             path_res = re.search(r'WRS_PATH = ([0-9]+)',mtl_, re.I)
593
        #             path_str = "%03d"%int(path_res.group(1)) if path_res!=None else '000'
594 595
        #             row_res  = re.search(r'STARTING_ROW = ([0-9]+)',mtl_, re.I)
        #             if row_res is None: row_res = re.search(r'WRS_ROW = ([0-9]+)',mtl_, re.I)
596 597 598
        #             row_str  = "%03d"%int(row_res.group(1)) if row_res!=None else '000'
        #             tt       = datetime.datetime.strptime(self.AcqDate, '%Y-%m-%d').timetuple()
        #             julianD  = '%d%03d'%(tt.tm_year,tt.tm_yday)
599 600
        #             station_res  = re.search(r'GROUND_STATION = "([\S]+)"',mtl_, re.I)
        #             if station_res is None: station_res  = re.search(r'STATION_ID = "([\S]+)"',mtl_, re.I)
601 602 603 604 605 606
        #             station_str = station_res.group(1) if station_res is not None else 'XXX'
        #             idWithoutVersion = 'L%s%s%s%s%s%s'%(sen,sat,path_str,row_str,julianD,station_str)
        #             filt_webmeta     = [i for i in webmeta if i['entityId'].startswith(idWithoutVersion)]
        #             if len(filt_webmeta) == 1:
        #                 self.EntityID = filt_webmeta[0]['entityId']

607 608 609
        ##########################
        # band specific metadata #
        ##########################
610
        LBA_full_sorted = natsorted(self.LayerBandsAssignment_full)
611 612

        # Gains and Offsets
613
        h4 = re.search(r"GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
614
        h4_ = h4.group(0)
615 616 617 618
        h4max_rad = re.findall(r" LMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4min_rad = re.findall(r" LMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4max_pix = re.findall(r"QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall(r"QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
619
        if not h4max_rad:
620
            h4max_rad = re.findall(r" RADIANCE_MAXIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
621
                                   re.I)  # space in front IS needed
622
            h4min_rad = re.findall(r" RADIANCE_MINIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
623
                                   re.I)  # space in front IS needed
624 625
            h4max_pix = re.findall(r"QUANTIZE_CAL_MAX_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
            h4min_pix = re.findall(r"QUANTIZE_CAL_MIN_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646
        # additional: LMAX, LMIN, QCALMAX, QCALMIN
        self.additional.append([["LMAX"], ["LMIN"], ["QCALMAX"], ["QCALMIN"]])
        # Offsets + Gains
        Gains, Offsets = [], []
        for ii, ind in enumerate(h4min_rad):
            Gains.append(
                (float(h4max_rad[ii]) - float(h4min_rad[ii])) / (float(h4max_pix[ii]) - float(h4min_pix[ii])))
            Offsets.append(float(ind) - float(h4min_pix[ii]) * Gains[-1])
            self.additional[-1][-4].append(h4max_rad[ii])
            self.additional[-1][-3].append(h4min_rad[ii])
            self.additional[-1][-2].append(h4max_pix[ii])
            self.additional[-1][-1].append(h4min_pix[ii])
        self.Gains = {bN: Gains[i] for i, bN in enumerate(LBA_full_sorted)}
        self.Offsets = {bN: Offsets[i] for i, bN in enumerate(LBA_full_sorted)}

        # Solar irradiance, central wavelengths , full width half maximum per band
        self.wvlUnit = 'Nanometers'
        # Exact values calculated based in SRFs
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()  # 3x dict
        # Provider values
        if not self.SolIrradiance:  # Preconfigured Irradiation and CWL values
647
            if re.search(r"[\D]*5", self.Satellite):
648
                # Landsat 5; resource for center wavelength (6 = thermal)
649 650 651 652 653 654 655 656 657 658 659 660 661 662 663
                self.SolIrradiance = {'1': 1957.,
                                      '2': 1826.,
                                      '3': 1554.,
                                      '4': 1036.,
                                      '5': 215.,
                                      '6': 0.0,
                                      '7': 80.67}
                self.CWL = {'1': 485.,
                            '2': 560.,
                            '3': 660.,
                            '4': 830.,
                            '5': 1650.,
                            '6': 11450.,
                            '7': 2215.}
            if re.search(r"[\D]*7", self.Satellite):
664 665 666
                # Landsat 7; resource for center wavelength:
                # http://opticks.org/confluence/display/opticksDev/Sensor+Wavelength+Definitions
                # 6L(thermal),6H(thermal),B8(pan)
667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685
                self.SolIrradiance = {'1': 1969.,
                                      '2': 1840.,
                                      '3': 1551.,
                                      '4': 1044.,
                                      '5': 225.7,
                                      '6L': 0.0,
                                      '6H': 0.0,
                                      '7': 82.07,
                                      '8': 1368.}
                self.CWL = {'1': 482.5,
                            '2': 565.,
                            '3': 660.,
                            '4': 825.,
                            '5': 1650.,
                            '6L': 11450.,
                            '6H': 11450.,
                            '7': 2215.,
                            '8': 710.}
            if re.search(r"[\D]*8", self.Satellite):  # Landsat 8
686
                # no irradiation values available
687 688 689 690 691 692 693 694 695 696 697 698
                self.CWL = {'1': 443.,
                            '2': 482.6,
                            '3': 561.3,
                            '4': 654.6,
                            '5': 864.6,
                            '6': 1609.1,
                            '7': 2201.2,
                            '8': 591.7,
                            '9': 1373.5,
                            '10': 10903.6,
                            '11': 12003.}
        # if None in SolIrradiance:
699 700

        # Thermal Constants (K1 = [W*m-2*um-1]; K1 = [K])
701
        if re.search(r"[\D]*5", self.Satellite):
702 703 704
            # resource: http://www.yale.edu/ceo/Documentation/Landsat_DN_to_Kelvin.pdf
            ThermalConstK1 = [0., 0., 0., 0., 0., 607.76, 0.]
            ThermalConstK2 = [0., 0., 0., 0., 0., 1260.56, 0.]
705
        if re.search(r"[\D]*7", self.Satellite):
706 707 708
            # resource: http://www.yale.edu/ceo/Documentation/Landsat_DN_to_Kelvin.pdf
            ThermalConstK1 = [0., 0., 0., 0., 0., 666.09, 666.09, 0., 0.]
            ThermalConstK2 = [0., 0., 0., 0., 0., 1282.71, 1282.71, 0., 0.]
709 710 711
        if re.search(r"[\D]*8", self.Satellite):  # Landsat 8
            K1_res = re.findall(r"K1_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
            K2_res = re.findall(r"K2_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
712 713 714 715 716 717 718 719 720 721
            if len(K1_res) == 2:
                ThermalConstK1 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K1_res[0]), float(K1_res[1])]
                ThermalConstK2 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K2_res[0]), float(K2_res[1])]
            else:
                self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. '
                                  'Found %s' % len(K1_res))
        self.ThermalConstK1 = {bN: ThermalConstK1[i] for i, bN in enumerate(LBA_full_sorted)}
        self.ThermalConstK2 = {bN: ThermalConstK2[i] for i, bN in enumerate(LBA_full_sorted)}

        # reflectance coefficients (Landsat8)
722
        h9 = re.search(r"GROUP = RADIOMETRIC_RESCALING[\s\S]*END_GROUP = RADIOMETRIC_RESCALING", mtl_, re.I)
723 724
        if h9:
            h9_ = h9.group(0)
725 726 727 728
            h9gain_ref = re.findall(r" REFLECTANCE_MULT_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9gain_ref_bandNr = re.findall(r" REFLECTANCE_MULT_BAND_([0-9]+) = [\S]*", h9_, re.I)
            h9offs_ref = re.findall(r" REFLECTANCE_ADD_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9offs_ref_bandNr = re.findall(r" REFLECTANCE_ADD_BAND_([0-9]+) = [\S]*", h9_, re.I)
729 730 731 732 733 734 735 736 737 738 739 740
            if h9gain_ref:
                GainsRef = [None] * len(self.Gains)
                OffsetsRef = [None] * len(self.Offsets)

                for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref):
                    GainsRef[LBA_full_sorted.index(ii)] = ind
                for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref):
                    OffsetsRef[LBA_full_sorted.index(ii)] = ind

                self.GainsRef = {bN: GainsRef[i] for i, bN in enumerate(LBA_full_sorted)}
                self.OffsetsRef = {bN: OffsetsRef[i] for i, bN in enumerate(LBA_full_sorted)}

741 742 743
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)
744 745 746 747 748 749 750

        # mtl.close()

    def Read_RE_metaxml(self):
        """----METHOD_4------------------------------------------------------------
        read metadata from RapidEye metafile: <dataname>metadata.xml
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
751

752 753
        # self.default_attr()
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
754 755
            glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/*_metadata.xml'))
            assert len(glob_res) > 0, 'No *metadata.xml file can be found in %s/*/!' % self.FolderOrArchive
756 757
            self.Metafile = glob_res[0]
            mxml_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
758
        else:  # archive
759
            mxml_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*_metadata.xml')
760

761
        # EntityID
762
        h1 = re.search(r"<[a-z]*:identifier>([\S]*)</[a-z]*:identifier>", mxml_, re.I)
763
        self.EntityID = h1.group(1) if h1 else "Not found"
764

765
        # Processing Level
766
        h2 = re.search(r"<[a-z]*:productType>([a-zA-Z0-9]*)</[a-z]*:productType>", mxml_, re.I)
767
        self.ProcLCode = h2.group(1) if h2 else "Not found"
768

769