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
        # Satellite
770
771
        h3 = re.search(r"<[a-z]*:serialIdentifier>([a-zA-Z0-9-]*)</[a-z]*:serialIdentifier>", mxml_, re.I)
        self.Satellite = 'RapidEye-%s' % re.search(r'[\D]*([0-9])', h3.group(1), re.I).group(1) if h3 else "Not found"
772

773
        # Sensor (Instrument shortname)
774
        h4 = re.search(r"<[a-z]*:Instrument>[\s]*<eop:shortName>([\S]*)</[a-z]*:shortName>", mxml_, re.I)
775
        self.Sensor = h4.group(1) if h4 else "Not found"
776

777
        # Acquisition Parameters: Incidence Angle, SunAzimuth, SunElevation, ViewingAngle, FieldOfView, AcqDate, AcqTime
778
        try:
779
780
781
782
783
784
785
786
787
788
789
790
791
            h5 = re.search(r'<[a-z]*:incidenceAngle uom="deg">([\S]*)'
                           r'</[a-z]*:incidenceAngle>[\s]*'
                           r'<opt:illuminationAzimuthAngle uom="deg">([\S]*)'
                           r'</opt:illuminationAzimuthAngle>[\s]*'
                           r'<opt:illuminationElevationAngle uom="deg">([\S]*)'
                           r'</opt:illuminationElevationAngle>[\s]*'
                           r'<re:azimuthAngle uom="deg">([\S]*)'
                           r'</re:azimuthAngle>[\s]*'
                           r'<re:spaceCraftViewAngle uom="deg">([\S]*)'
                           r'</re:spaceCraftViewAngle>[\s]*'
                           r'<re:acquisitionDateTime>([0-9-]*)T([0-9:]*)[\S]*'
                           r'</re:acquisitionDateTime>',
                           mxml_, re.I)
792
            self.IncidenceAngle = float(h5.group(1))
Daniel Scheffler's avatar
Daniel Scheffler committed
793
            self.SunAzimuth = float(h5.group(2))
794
            self.SunElevation = float(h5.group(3))
Daniel Scheffler's avatar
Daniel Scheffler committed
795
796
797
            # angle from true north at the image or tile center to the scan (line) direction at image center,
            # in clockwise positive degrees.
            self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))])
798
            self.ViewingAngle = float(h5.group(5))
Daniel Scheffler's avatar
Daniel Scheffler committed
799
            self.FOV = get_FieldOfView(self.GMS_identifier)
800
801
            AcqDate = h5.group(6)
            AcqTime = h5.group(7)
802
803

        except AttributeError:
804
805
806
807
808
809
810
811
812
813
814
815
            h5 = re.search(r'<hma:acrossTrackIncidenceAngle uom="deg">([\S]*)'
                           r'</hma:acrossTrackIncidenceAngle>[\s]*'
                           r'<ohr:illuminationAzimuthAngle uom="deg">([\S]*)'
                           r'</ohr:illuminationAzimuthAngle>[\s]*'
                           r'<ohr:illuminationElevationAngle uom="deg">([\S]*)'
                           r'</ohr:illuminationElevationAngle>[\s]*'
                           r'<re:azimuthAngle uom="deg">([\S]*)'
                           r'</re:azimuthAngle>[\s]*'
                           r'<re:acquisitionDateTime>([0-9-]*)'
                           r'T([0-9:]*)[\S]*'
                           r'</re:acquisitionDateTime>',
                           mxml_, re.I)  #
816
817
818
            self.IncidenceAngle = float(h5.group(1))
            self.SunAzimuth = float(h5.group(2))
            self.SunElevation = float(h5.group(3))
Daniel Scheffler's avatar
Daniel Scheffler committed
819
            self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))])
820
            self.ViewingAngle = 9999.99
Daniel Scheffler's avatar
Daniel Scheffler committed
821
            self.FOV = get_FieldOfView(self.GMS_identifier)
822
823
824
825
826
            AcqDate = h5.group(5)
            AcqTime = h5.group(6)
        # 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')
827

828
        # Earth sun distance
829
830
831
        self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate)

        # Additional: Projection
832
833
834
835
        h6 = re.search(r"<re:geodeticDatum>([\S]*)"
                       r"</re:geodeticDatum>[\s]*"
                       r"<re:projection>([\S^<\s]*)"
                       r"</re:projection>",
Daniel Scheffler's avatar
Daniel Scheffler committed
836
                       mxml_, re.I)
837
        try:
838
839
840
841
            self.additional.append(["SpatialReference",
                                    ["geodeticDatum", h6.group(1)],
                                    ["projection", h6.group(2)]
                                    ])
Daniel Scheffler's avatar
Daniel Scheffler committed
842
        except AttributeError:  # Level1-B data has no geodaticDatum
843
844
            pass

845
        # Corrections applied + Quality
846
847
848
849
850
851
852
853
854
855
856
857
        h7 = re.search(r"<re:radiometricCorrectionApplied>([\w]*)"
                       r"</re:radiometricCorrectionApplied>[\s]*"
                       r"<re:radiometricCalibrationVersion>([\S]*)"
                       r"</re:radiometricCalibrationVersion>[\s]*"
                       r"<re:geoCorrectionLevel>([\S\s^<]*)"
                       r"</re:geoCorrectionLevel>[\s]*"
                       r"<re:elevationCorrectionApplied>([\S]*)"
                       r"</re:elevationCorrectionApplied>[\s]*"
                       r"<re:atmosphericCorrectionApplied>([\w]*)"
                       r"</re:atmosphericCorrectionApplied>[\s]*"
                       r"<re:productAccuracy>([\S\s^<]*)"
                       r"</re:productAccuracy>", mxml_, re.I)
858
859
        # fuer IRIS ihre Daten
        if h7 is None:
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
            h7 = re.search(r"<re:radiometricCorrectionApplied>([\w]*)"
                           r"</re:radiometricCorrectionApplied>[\s]*"
                           r"<re:geoCorrectionLevel>([\S\s^<]*)"
                           r"</re:geoCorrectionLevel>[\s]*"
                           r"<re:elevationCorrectionApplied>([\S]*)"
                           r"</re:elevationCorrectionApplied>[\s]*"
                           r"<re:atmosphericCorrectionApplied>([\w]*)"
                           r"</re:atmosphericCorrectionApplied>",
                           mxml_, re.I)
            self.additional.append(
                ["Corrections",
                 ["radiometricCorrectionApplied", h7.group(1)],
                 ["geoCorrectionLevel", h7.group(2)],
                 ["elevationCorrectionApplied", h7.group(3)],
                 ["atmosphericCorrectionApplied", h7.group(4)]
                 ])
876
        else:
877
878
879
880
881
882
883
884
            self.additional.append(
                ["Corrections",
                 ["radiometricCorrectionApplied", h7.group(1)],
                 ["radiometricCalibrationVersion", h7.group(2)],
                 ["geoCorrectionLevel", h7.group(3)],
                 ["elevationCorrectionApplied", h7.group(4)],
                 ["atmosphericCorrectionApplied", h7.group(5)]
                 ])
Daniel Scheffler's avatar
Daniel Scheffler committed
885
886
            self.Quality.append(["geomProductAccuracy[m]:", str(
                round(float(h7.group(6)), 1))])  # Estimated product horizontal CE90 uncertainty [m]
887

Daniel Scheffler's avatar