metadata.py 119 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
"""Module 'metadata' for handling any type of metadata of GeoMultiSens compatible sensor systems."""
Daniel Scheffler's avatar
Daniel Scheffler committed
3
4

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

6
7
import collections
import datetime
8
import glob
9
10
11
12
import math
import os
import re
import sys
13
import warnings
14
import iso8601
15
import xml.etree.ElementTree as ET
16
from typing import List, TYPE_CHECKING  # noqa F401  # flake8 issue
17
18
19

import numpy as np
import pyproj
20
21
from matplotlib import dates as mdates
from pyorbital import astronomy
22

Daniel Scheffler's avatar
Daniel Scheffler committed
23
24
from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.geo.projection import WKT2EPSG
25
from sicor.options import get_options as get_ac_options
26

27
28
29
30
31
32
33
34
35
from ..options.config import GMS_config as CFG
from ..io.input_reader import open_specific_file_within_archive, Solar_Irradiance_reader, SRF_reader
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
36

37
38
if TYPE_CHECKING:
    from ..model.gms_object import GMS_identifier  # noqa F401  # flake8 issue
39

40
41
__author__ = 'Daniel Scheffler', 'Robert Behling'

42

43
class METADATA(object):
44
    def __init__(self, GMS_id):
45
46
47
        # private attributes
        self._AcqDateTime = None

48
        # unpack GMS_identifier
49
50
51
52
53
54
55
        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
56
57
58

        self.Dataname = ''
        self.FolderOrArchive = ''
Daniel Scheffler's avatar
Daniel Scheffler committed
59
60
61
        self.Metafile = ""  # File containing image metadata (automatically found)
        self.EntityID = ""  # ID to identify the original scene
        self.SceneID = ''  # postgreSQL-database identifier
62
        self.Sensormode = ""
Daniel Scheffler's avatar
Daniel Scheffler committed
63
64
65
        self.gResolution = -99.  # resolution [m]
        self.AcqDate = ""  # YYYY-MM-DD
        self.AcqTime = ""  # HH:MM:SS
66
67
68
69
70
71
        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
72
73
74
75
76
        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
77
        self.SolIrradiance = []
78
79
        self.ThermalConstK1 = {}
        self.ThermalConstK2 = {}
80
        self.EarthSunDist = 1.0
Daniel Scheffler's avatar
Daniel Scheffler committed
81
82
83
        # viewing zenith angle of the Sensor (offNadir) [Deg] (usually) in Case of RapidEye "+"
        # being East and “-” being West
        self.ViewingAngle = -99.0
84
        self.ViewingAngle_arrProv = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
85
86
87
        # 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
88
        self.IncidenceAngle_arrProv = {}
Daniel Scheffler's avatar
Daniel Scheffler committed
89
90
91
        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 = []
92
        self.PhysUnit = "DN"
Daniel Scheffler's avatar
Daniel Scheffler committed
93
94
95
96
        # 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
97
98
99
100
101
        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
102
103
        # List of Corner Coordinates in order of Lon/Lat/DATA_X/Data_Y for all given image corners
        self.CornerTieP_LonLat = []
104
        self.CornerTieP_UTM = []
Daniel Scheffler's avatar
Daniel Scheffler committed
105
        self.LayerBandsAssignment = []  # List of spectral bands in the same order as they are stored in the rasterfile.
106
107
108
109
        self.additional = []
        self.image_type = 'RSD'
        self.orbitParams = {}
        self.overpassDurationSec = -99.
Daniel Scheffler's avatar
Daniel Scheffler committed
110
        self.scene_length = -99.
111
112
113
        self.rows = -99.
        self.cols = -99.
        self.bands = -99.
114
        self.nBands = -99.
115
116
117
        self.map_info = []
        self.projection = ""
        self.wvlUnit = ""
Daniel Scheffler's avatar
Daniel Scheffler committed
118
        self.spec_vals = {'fill': None, 'zero': None, 'saturated': None}
119

120
121
122
        self.version_gms_preprocessing = CFG.version
        self.versionalias_gms_preprocessing = CFG.versionalias

123
124
125
126
127
128
129
130
131
132
133
134
135
136
    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 []

        if re.search("SPOT", self.Satellite, re.I):
            self.Read_SPOT_dimap2()
        elif re.search("Terra", self.Satellite, re.I):
            self.Read_ASTER_hdffile(self.Subsystem)
137
138
        elif re.search("Sentinel-2A", self.Satellite, re.I) or re.search("Sentinel-2B", self.Satellite, re.I):
            self.Read_Sentinel2_xmls()
139
140
141
142
143
144
145
146
        elif re.search("LANDSAT", self.Satellite, re.I):
            self.Read_LANDSAT_mtltxt(self.LayerBandsAssignment)
        elif re.search("RapidEye", self.Satellite, re.I):
            self.Read_RE_metaxml()
        elif re.search("ALOS", self.Satellite, re.I):
            self.Read_ALOS_summary()
            self.Read_ALOS_LEADER()  # for gains and offsets
        else:
147
            raise NotImplementedError("%s is not a supported sensor or sensorname is misspelled." % self.Satellite)
148
149

        return self
150

151
152
153
    @property
    def AcqDateTime(self):
        """Returns a datetime.datetime object containing date, time and timezone (UTC time)."""
154
155
156
        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
157

158
159
160
161
162
        return self._AcqDateTime

    @AcqDateTime.setter
    def AcqDateTime(self, DateTime):
        # type: (datetime.datetime) -> None
163
164
165
166
167
168
169
        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')
170

171
172
173
    @property
    def overview(self):
        attr2include = \
Daniel Scheffler's avatar
Daniel Scheffler committed
174
175
            ['Dataname', 'Metafile', 'EntityID', 'Satellite', 'Sensor', 'Subsystem', 'gResolution', 'AcqDate',
             'AcqTime',
176
177
             'CWL', 'FWHM', 'Offsets', 'Gains', 'ProcLCode', 'ProcLStr', 'SunElevation', 'SunAzimuth',
             'ViewingAngle', 'IncidenceAngle', 'FOV', 'SolIrradiance', 'ThermalConstK1', 'ThermalConstK2',
Daniel Scheffler's avatar
Daniel Scheffler committed
178
179
             'EarthSunDist', 'Quality', 'PhysUnit', 'additional', 'GainsRef', 'OffsetsRef', 'CornerTieP_LonLat',
             'CS_EPSG',
180
181
182
             'CS_TYPE', 'CS_DATUM', 'CS_UTM_ZONE', 'LayerBandsAssignment']
        return collections.OrderedDict((attr, getattr(self, attr)) for attr in attr2include)

183
184
185
186
187
188
189
190
191
    @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)

192
193
194
195
    def Read_SPOT_dimap2(self):
        """----METHOD_2------------------------------------------------------------
        # read metadata from spot dimap file
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
196

Daniel Scheffler's avatar
Daniel Scheffler committed
197
        # self.default_attr()
198
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
199
200
            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
201
202
            self.Metafile = glob_res[0]
            dim_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
203
        else:  # archive
204
            dim_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/scene01/metadata.dim')
205

Daniel Scheffler's avatar
Daniel Scheffler committed
206
207
208
        # special values
        h1 = re.findall("<SPECIAL_VALUE_INDEX>([a-zA-Z0-9]*)</SPECIAL_VALUE_INDEX>", dim_, re.I)
        h2 = re.findall("<SPECIAL_VALUE_TEXT>([a-zA-Z0-9]*)</SPECIAL_VALUE_TEXT>", dim_, re.I)
209
        self.additional.append(["SpecialValues:"])
Daniel Scheffler's avatar
Daniel Scheffler committed
210
211
        for ii, ind in enumerate(h1):
            self.additional[0].append(["%s:%s" % (ind, h2[ii])])
212

Daniel Scheffler's avatar
Daniel Scheffler committed
213
214
        # EntityID
        h3 = re.search("<SOURCE_ID>([a-zA-Z0-9]*)</SOURCE_ID>", dim_, re.I)
215
216
        self.EntityID = h3.group(1)

Daniel Scheffler's avatar
Daniel Scheffler committed
217
218
        # AcqDate
        h4 = re.search("<IMAGING_DATE>([0-9-]*)</IMAGING_DATE>", dim_, re.I)
219
        AcqDate = h4.group(1)
220

Daniel Scheffler's avatar
Daniel Scheffler committed
221
        # Earth sun distance
222
        self.EarthSunDist = self.get_EarthSunDistance(AcqDate)
223

Daniel Scheffler's avatar
Daniel Scheffler committed
224
225
        # AcqTime
        h5 = re.search("<IMAGING_TIME>([0-9:]*)</IMAGING_TIME>", dim_, re.I)
226
227
        AcqTime = h5.group(1)

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

Daniel Scheffler's avatar
Daniel Scheffler committed
232
        # Satellite + Sensor
233
234
        h6 = re.search("<MISSION>([a-zA-Z]*)</MISSION>[a-zA-Z0-9\s]*<MISSION_INDEX>([0-9]*)</MISSION_INDEX>"
                       "[a-zA-Z0-9\s]*<INSTRUMENT>([a-zA-Z]*)</INSTRUMENT>[a-zA-Z0-9\s]*<INSTRUMENT_INDEX>([0-9]*)"
Daniel Scheffler's avatar
Daniel Scheffler committed
235
236
237
                       "</INSTRUMENT_INDEX>[a-zA-Z0-9\s]*<SENSOR_CODE>([a-zA-Z0-9]*)</SENSOR_CODE>", dim_, re.I)
        self.Satellite = "-".join([h6.group(1), h6.group(2)])
        self.Sensor = "".join([h6.group(3), h6.group(4), h6.group(5)])
238

Daniel Scheffler's avatar
Daniel Scheffler committed
239
        # Angles: incidence angle, sunAzimuth, sunElevation
240
        h7 = re.search("<INCIDENCE_ANGLE>(.*)</INCIDENCE_ANGLE>[\s\S]*<SUN_AZIMUTH>(.*)</SUN_AZIMUTH>[\s\S]"
Daniel Scheffler's avatar
Daniel Scheffler committed
241
                       "*<SUN_ELEVATION>(.*)</SUN_ELEVATION>", dim_, re.I)
242
243
244
245
        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
246
247
        # Viewing Angle: not always included in the metadata.dim file
        h8 = re.search("<VIEWING_ANGLE>(.*)</VIEWING_ANGLE>", dim_, re.I)
248
        if type(h8).__name__ == 'NoneType':
Daniel Scheffler's avatar
Daniel Scheffler committed
249
            self.ViewingAngle = 0
250
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
251
            self.ViewingAngle = float(h8.group(1))
252

Daniel Scheffler's avatar
Daniel Scheffler committed
253
        # Field of View
254
        self.FOV = get_FieldOfView(self.GMS_identifier)
255

Daniel Scheffler's avatar
Daniel Scheffler committed
256
        # Units
257
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
258
        self.PhysUnit = "DN"
259

260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
        # ProcLevel
        h11 = re.search("<PROCESSING_LEVEL>([a-zA-Z0-9]*)</PROCESSING_LEVEL>", dim_, re.I)
        self.ProcLCode = h11.group(1)

        # Quality
        h12 = re.findall("<Bad_Pixel>[\s]*<PIXEL_INDEX>([0-9]*)</PIXEL_INDEX>[\s]*<BAD_PIXEL_STATUS>([^<]*)"
                         "</BAD_PIXEL_STATUS></Bad_Pixel>", dim_, re.I)
        self.Quality = h12

        # Coordinate Reference System
        re_CS_EPSG = re.search('<HORIZONTAL_CS_CODE>epsg:([0-9]*)</HORIZONTAL_CS_CODE>', dim_, re.I)
        re_CS_TYPE = re.search('<HORIZONTAL_CS_TYPE>([a-zA-Z0-9]*)</HORIZONTAL_CS_TYPE>', dim_, re.I)
        re_CS_DATUM = re.search('<HORIZONTAL_CS_NAME>([\w+\s]*)</HORIZONTAL_CS_NAME>', dim_, re.I)
        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
        h121 = re.findall("<TIE_POINT_CRS_X>(.*)</TIE_POINT_CRS_X>", dim_, re.I)
        h122 = re.findall("<TIE_POINT_CRS_Y>(.*)</TIE_POINT_CRS_Y>", dim_, re.I)
        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 #
        ##########################

        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

Daniel Scheffler's avatar
Daniel Scheffler committed
296
297
        # Gains and Offsets
        h9 = re.search("<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I)
298
299
300
301
302
        h9_ = h9.group(0)
        h91 = re.findall("<PHYSICAL_UNIT>([^<]*)</PHYSICAL_UNIT>", h9_, re.I)
        h92 = re.findall("<PHYSICAL_BIAS>([^<]*)</PHYSICAL_BIAS>", h9_, re.I)
        h93 = re.findall("<PHYSICAL_GAIN>([^<]*)</PHYSICAL_GAIN>", h9_, re.I)
        self.additional.append(["Physical Units per band:"])
303
        for ii, ind in enumerate(h91):  # FIXME does not respect sort_bands_by_cwl
304
            self.additional[-1].append(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
305
        # Offsets
306
307
        for ii, (ind, band) in enumerate(zip(h92, LBA_full_sorted)):
            self.Offsets[band] = float(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
308
        # Gains
309
        for ii, (ind, band) in enumerate(zip(h93, LBA_full_sorted)):
Daniel Scheffler's avatar
Daniel Scheffler committed
310
            # gains in dimap file represent reciprocal 1/gain
311
            self.Gains[band] = 1 / float(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
312
313
314
315

        # 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
316
        # noinspection PyBroadException
Daniel Scheffler's avatar
Daniel Scheffler committed
317
318
319
320
        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)
321
        except Exception:
Daniel Scheffler's avatar
Daniel Scheffler committed
322
323
            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)
324
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)
325

Daniel Scheffler's avatar
Daniel Scheffler committed
326
        # Exact values calculated based in SRFs
327
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
328
        # Provider values
329
        if not self.SolIrradiance:
Daniel Scheffler's avatar
Daniel Scheffler committed
330
            h10 = re.search("<Solar_Irradiance>[\s\S]*</Solar_Irradiance>", dim_, re.I)
331
332
333
            h10_ = h10.group(0)
            h101 = re.findall("<SOLAR_IRRADIANCE_VALUE>([^<]*)</SOLAR_IRRADIANCE_VALUE>", h10_, re.I)
            if h101:
334
                self.SolIrradiance = dict(zip(LBA_full_sorted, h101))
Daniel Scheffler's avatar
Daniel Scheffler committed
335
336
337
338
                #    self.additional.append(["Solar Irradiance per band:"])
                #    for ii,ind in enumerate(h101):
                #       self.additional[-1].append(ind)
            else:  # Preconfigured Irradiation values
339
                spot_irr_dic = {
340
341
342
343
344
345
346
347
348
349
                    '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.]))}
350
                self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)]
Daniel Scheffler's avatar
Daniel Scheffler committed
351
            # Preconfigured CWLs --   # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined
352
353
            # 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
354
355
356
357
358
359
360
            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.]))
361

362
363
364
365
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

366
367
368
369
370
    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
371

372
373
        # self.default_attr()
        self.LayerBandsAssignment = LayerBandsAssignment
Daniel Scheffler's avatar
Daniel Scheffler committed
374
        self.nBands = len(LayerBandsAssignment)
375
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
376
377
            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
378
379
            self.Metafile = glob_res[0]
            mtl_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
380
        else:  # archive
381
            mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt')
382

383
        # Processing Level
Daniel Scheffler's avatar
Daniel Scheffler committed
384
        h1 = re.search('PRODUCT_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
385
        if h1 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
386
            h1 = re.search('DATA_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
387
388
389
        self.ProcLCode = h1.group(1)

        # Satellite + Sensor + Sensor Mode
Daniel Scheffler's avatar
Daniel Scheffler committed
390
391
        h2 = re.search('SPACECRAFT_ID = "([a-zA-Z0-9_]*)"[\s]*SENSOR_ID = "([a-zA-Z0-9+]*)"[\s]*'
                       'SENSOR_MODE = "([\S]*)"', mtl_, re.I)
392
        if h2:
Daniel Scheffler's avatar
Daniel Scheffler committed
393
394
            self.Satellite = 'Landsat-%s' % re.search('LANDSAT[\D]*([0-9])', h2.group(1), re.I).group(1)
            self.Sensor = h2.group(2)
395
396
            self.Sensormode = h2.group(3)
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
397
398
399
400
401
402
403
            h2a = re.search('SPACECRAFT_ID = "([a-zA-Z0-9_]*)"', mtl_, re.I)
            h2b = re.search('SENSOR_ID = "([a-zA-Z0-9_+]*)"', mtl_, re.I)
            h2c = re.search('SENSOR_MODE = "([a-zA-Z0-9_+]*)"', mtl_, re.I)
            self.Satellite = 'Landsat-%s' % re.search('LANDSAT[\D]*([0-9])', h2a.group(1), re.I).group(1)
            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
404
405

        # EntityID
Daniel Scheffler's avatar
Daniel Scheffler committed
406
        h2c = re.search('LANDSAT_SCENE_ID = "([A-Z0-9]*)"', mtl_, re.I)
407
        if h2c:
408
409
            self.EntityID = h2c.group(1)

410
        # Acquisition Date + Time
Daniel Scheffler's avatar
Daniel Scheffler committed
411
        h3 = re.search('ACQUISITION_DATE = ([0-9-]*)[\s]*SCENE_CENTER_SCAN_TIME = "?([0-9:]*)"?', mtl_, re.I)
412
        if h3 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
413
            h3 = re.search('DATE_ACQUIRED = ([0-9-]*)[\s]*SCENE_CENTER_TIME = "?([0-9:]*)"?', mtl_, re.I)
414
415
416
417
418
419
        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')
420
421
422
423

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

424
425
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
426
        self.PhysUnit = "DN"
427

428
        # Angles: incidence angle, sunAzimuth, sunElevation, field of view
Daniel Scheffler's avatar
Daniel Scheffler committed
429
430
        h5 = re.search("SUN_AZIMUTH = ([\S]*)[\s]*SUN_ELEVATION = ([\S]*)", mtl_, re.I)
        self.SunAzimuth = float(h5.group(1))
431
        self.SunElevation = float(h5.group(2))
Daniel Scheffler's avatar
Daniel Scheffler committed
432
        self.FOV = get_FieldOfView(self.GMS_identifier)
433
434

        # Quality
Daniel Scheffler's avatar
Daniel Scheffler committed
435
        h6 = re.search("GROUP = CORRECTIONS_APPLIED[\s\S]*END_GROUP = CORRECTIONS_APPLIED", mtl_, re.I)
436
        if h6 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
437
            h6 = re.search("GROUP = IMAGE_ATTRIBUTES[\s\S]*END_GROUP = IMAGE_ATTRIBUTES", mtl_, re.I)
438

Daniel Scheffler's avatar
Daniel Scheffler committed
439
        h6_ = h6.group(0)
440
        h61 = (h6_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
441
        x = 0
442
        for i in h61:
Daniel Scheffler's avatar
Daniel Scheffler committed
443
            if x == 0 or x + 1 == len(h61):
444
445
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
446
                i_ = i.strip().replace('"', "")
447
                self.Quality.append(i_.split(" = "))
Daniel Scheffler's avatar
Daniel Scheffler committed
448
            x += 1
449

450
        # Additional: coordinate system, geom. Resolution
Daniel Scheffler's avatar
Daniel Scheffler committed
451
452
        h7 = re.search("GROUP = PROJECTION_PARAMETERS[\s\S]*END_GROUP = L1_METADATA_FILE", mtl_, re.I)
        h7_ = h7.group(0)
453
        h71 = (h7_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
454
455
        for x, i in enumerate(h71):
            if re.search("Group", i, re.I):
456
457
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
458
                i_ = i.strip().replace('"', "")
459
                self.additional.append(i_.split(" = "))
Daniel Scheffler's avatar
Daniel Scheffler committed
460
461
462
463
464
465
466
467
        re_CS_TYPE = re.search('MAP_PROJECTION = "([\w+\s]*)"', h7_, re.I)
        re_CS_DATUM = re.search('DATUM = "([\w+\s]*)"', h7_, re.I)
        re_CS_UTM_ZONE = re.search('ZONE_NUMBER = ([0-9]*)\n', mtl_, re.I)
        re_CS_UTM_ZONE = re.search('UTM_ZONE = ([0-9]*)\n', h7_,
                                   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
468
469
470
        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
471
472
        # Landsat8
        h8 = re.search("ROLL_ANGLE = ([\S]*)", mtl_, re.I)
473
        if h8:
474
475
476
477
478
479
480
            self.ViewingAngle = float(h8.group(1))

        # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples)
        h10_UL = re.findall("PRODUCT_UL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_UR = re.findall("PRODUCT_UR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_LL = re.findall("PRODUCT_LL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
        h10_LR = re.findall("PRODUCT_LR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
481
        if not h10_UL:  # Landsat8
482
483
484
485
            h10_UL = re.findall("CORNER_UL_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_UR = re.findall("CORNER_UR_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_LL = re.findall("CORNER_LL_[\S]+ = (.*)[\S]*", mtl_, re.I)
            h10_LR = re.findall("CORNER_LR_[\S]+ = (.*)[\S]*", mtl_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
486
487
488
489
490
491
492
493
494
        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]]))
495
496

        # WRS path/row
Daniel Scheffler's avatar
Daniel Scheffler committed
497
        h11_p = re.search('WRS_PATH = ([0-9]*)', mtl_, re.I)
498
        if h11_p:
499
            self.WRS_path = h11_p.group(1)
Daniel Scheffler's avatar
Daniel Scheffler committed
500
501
502
503
        h11_r1 = re.search('STARTING_ROW = ([0-9]*)', mtl_, re.I)
        h11_r2 = re.search('ENDING_ROW = ([0-9]*)', mtl_, re.I)
        if h11_r1 is None:  # Landsat-8
            h11_r = re.search('WRS_ROW = ([0-9]*)', mtl_, re.I)
504
            self.WRS_row = int(h11_r.group(1))
505
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
506
            self.WRS_row = int(h11_r1.group(1)) if h11_r1.group(1) == h11_r2.group(1) else self.WRS_row
507
508

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

Daniel Scheffler's avatar
Daniel Scheffler committed
515
            if len(result) == 1:  # e.g. [('LE71950282003121EDC00',)]
516
517
                self.EntityID = result[0][0]
            elif len(result) == 0:
518
519
                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
520
            else:  # e.g. [('LE71950282003121EDC00',), ('LE71950282003105ASN00',)]
521
                self.logger.warning("Scene ID could not be assigned because %s datasets matching to the query "
Daniel Scheffler's avatar
Daniel Scheffler committed
522
                                    "parameters were found in metadata database." % len(result))
523
        # if self.EntityID=='':
Daniel Scheffler's avatar
Daniel Scheffler committed
524
525
        #     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 ''
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
        #     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]
        #             path_res = re.search('WRS_PATH = ([0-9]+)',mtl_, re.I)
        #             path_str = "%03d"%int(path_res.group(1)) if path_res!=None else '000'
        #             row_res  = re.search('STARTING_ROW = ([0-9]+)',mtl_, re.I)
        #             if row_res is None: row_res = re.search('WRS_ROW = ([0-9]+)',mtl_, re.I)
        #             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)
        #             station_res  = re.search('GROUND_STATION = "([\S]+)"',mtl_, re.I)
        #             if station_res is None: station_res  = re.search('STATION_ID = "([\S]+)"',mtl_, re.I)
        #             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']

550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
        ##########################
        # band specific metadata #
        ##########################
        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

        # Gains and Offsets
        h4 = re.search("GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
        h4_ = h4.group(0)
        h4max_rad = re.findall(" LMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4min_rad = re.findall(" LMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4max_pix = re.findall("QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall("QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        if not h4max_rad:
            h4max_rad = re.findall(" RADIANCE_MAXIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4min_rad = re.findall(" RADIANCE_MINIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4max_pix = re.findall("QUANTIZE_CAL_MAX_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
            h4min_pix = re.findall("QUANTIZE_CAL_MIN_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
        # 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
            if re.search("[\D]*5", self.Satellite):
                # Landsat 5; resource for center wavelength (6 = thermal)
                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("[\D]*7", self.Satellite):
                # Landsat 7; resource for center wavelength:
                # http://opticks.org/confluence/display/opticksDev/Sensor+Wavelength+Definitions
                # 6L(thermal),6H(thermal),B8(pan)
                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("[\D]*8", self.Satellite):  # Landsat 8
                # no irradiation values available
                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:

        # Thermal Constants (K1 = [W*m-2*um-1]; K1 = [K])
        if re.search("[\D]*5", self.Satellite):
            # 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.]
        if re.search("[\D]*7", self.Satellite):
            # 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.]
        if re.search("[\D]*8", self.Satellite):  # Landsat 8
            K1_res = re.findall("K1_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
            K2_res = re.findall("K2_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
            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)
        h9 = re.search("GROUP = RADIOMETRIC_RESCALING[\s\S]*END_GROUP = RADIOMETRIC_RESCALING", mtl_, re.I)
        if h9:
            h9_ = h9.group(0)
            h9gain_ref = re.findall(" REFLECTANCE_MULT_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9gain_ref_bandNr = re.findall(" REFLECTANCE_MULT_BAND_([0-9]+) = [\S]*", h9_, re.I)
            h9offs_ref = re.findall(" REFLECTANCE_ADD_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9offs_ref_bandNr = re.findall(" REFLECTANCE_ADD_BAND_([0-9]+) = [\S]*", h9_, re.I)
            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)}

649
650
651
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)
652
653
654
655
656
657
658

        # mtl.close()

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

660
661
        # self.default_attr()
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
662
663
            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
664
665
            self.Metafile = glob_res[0]
            mxml_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
666
        else:  # archive
667
            mxml_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*_metadata.xml')
668

669
        # EntityID
Daniel Scheffler's avatar
Daniel Scheffler committed
670
        h1 = re.search("<[a-z]*:identifier>([\S]*)</[a-z]*:identifier>", mxml_, re.I)
671
        self.EntityID = h1.group(1) if h1 else "Not found"
672

673
        # Processing Level
Daniel Scheffler's avatar
Daniel Scheffler committed
674
        h2 = re.search("<[a-z]*:productType>([a-zA-Z0-9]*)</[a-z]*:productType>", mxml_, re.I)
675
        self.ProcLCode = h2.group(1) if h2 else "Not found"
676

677
        # Satellite
Daniel Scheffler's avatar
Daniel Scheffler committed
678
679
        h3 = re.search("<[a-z]*:serialIdentifier>([a-zA-Z0-9-]*)</[a-z]*:serialIdentifier>", mxml_, re.I)
        self.Satellite = 'RapidEye-%s' % re.search('[\D]*([0-9])', h3.group(1), re.I).group(1) if h3 else "Not found"
680

681
        # Sensor (Instrument shortname)
Daniel Scheffler's avatar
Daniel Scheffler committed
682
        h4 = re.search("<[a-z]*:Instrument>[\s]*<eop:shortName>([\S]*)</[a-z]*:shortName>", mxml_, re.I)
683
        self.Sensor = h4.group(1) if h4 else "Not found"
684

685
        # Acquisition Parameters: Incidence Angle, SunAzimuth, SunElevation, ViewingAngle, FieldOfView, AcqDate, AcqTime
686
687
688
689
690
691
        try:
            h5 = re.search('<[a-z]*:incidenceAngle uom="deg">([\S]*)</[a-z]*:incidenceAngle>[\s]*<opt:illumination'
                           'AzimuthAngle uom="deg">([\S]*)</opt:illuminationAzimuthAngle>[\s]*'
                           '<opt:illuminationElevationAngle uom="deg">([\S]*)</opt:illuminationElevationAngle>[\s]*'
                           '<re:azimuthAngle uom="deg">([\S]*)</re:azimuthAngle>[\s]*'
                           '<re:spaceCraftViewAngle uom="deg">([\S]*)</re:spaceCraftViewAngle>[\s]*<re:acquisitionDate'
Daniel Scheffler's avatar
Daniel Scheffler committed
692
                           'Time>([0-9-]*)T([0-9:]*)[\S]*</re:acquisitionDateTime>', mxml_, re.I)
693
            self.IncidenceAngle = float(h5.group(1))
Daniel Scheffler's avatar
Daniel Scheffler committed
694
            self.SunAzimuth = float(h5.group(2))
695
            self.SunElevation = float(h5.group(3))
Daniel Scheffler's avatar
Daniel Scheffler committed
696
697
698
            # 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))])
699
            self.ViewingAngle = float(h5.group(5))
Daniel Scheffler's avatar
Daniel Scheffler committed
700
            self.FOV = get_FieldOfView(self.GMS_identifier)
701
702
            AcqDate = h5.group(6)
            AcqTime = h5.group(7)
703
704

        except AttributeError:
Daniel Scheffler's avatar
Daniel Scheffler committed
705
706
707
708
709
            h5 = re.search('<hma:acrossTrackIncidenceAngle uom="deg">([\S]*)</hma:acrossTrackIncidenceAngle>[\s]*'
                           '<ohr:illuminationAzimuthAngle uom="deg">([\S]*)</ohr:illuminationAzimuthAngle>[\s]*'
                           '<ohr:illuminationElevationAngle uom="deg">([\S]*)</ohr:illuminationElevationAngle>[\s]*'
                           '<re:azimuthAngle uom="deg">([\S]*)</re:azimuthAngle>[\s]*<re:acquisitionDateTime>([0-9-]*)'
                           'T([0-9:]*)[\S]*</re:acquisitionDateTime>', mxml_, re.I)  #
710
711
712
            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
713
            self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))])
714
            self.ViewingAngle = 9999.99
Daniel Scheffler's avatar
Daniel Scheffler committed
715
            self.FOV = get_FieldOfView(self.GMS_identifier)
716
717
718
719
720
            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')
721

722
        # Earth sun distance
723
724
725
        self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate)

        # Additional: Projection
Daniel Scheffler's avatar
Daniel Scheffler committed
726
727
        h6 = re.search("<re:geodeticDatum>([\S]*)</re:geodeticDatum>[\s]*<re:projection>([\S^<\s]*)</re:projection>",
                       mxml_, re.I)
728
        try:
Daniel Scheffler's avatar
Daniel Scheffler committed
729
730
            self.additional.append(["SpatialReference", ["geodeticDatum", h6.group(1)], ["projection", h6.group(2)]])
        except AttributeError:  # Level1-B data has no geodaticDatum
731
732
            pass

733
        # Corrections applied + Quality
Daniel Scheffler's avatar
Daniel Scheffler committed
734
735
736
737
738
739
        h7 = re.search("<re:radiometricCorrectionApplied>([\w]*)</re:radiometricCorrectionApplied>[\s]*"
                       "<re:radiometricCalibrationVersion>([\S]*)</re:radiometricCalibrationVersion>[\s]*"
                       "<re:geoCorrectionLevel>([\S\s^<]*)</re:geoCorrectionLevel>[\s]*"
                       "<re:elevationCorrectionApplied>([\S]*)</re:elevationCorrectionApplied>[\s]*"
                       "<re:atmosphericCorrectionApplied>([\w]*)</re:atmosphericCorrectionApplied>[\s]*"
                       "<re:productAccuracy>([\S\s^<]*)</re:productAccuracy>", mxml_, re.I)
740
741
        # fuer IRIS ihre Daten
        if h7 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
742
743
744
745
746
747
            h7 = re.search("<re:radiometricCorrectionApplied>([\w]*)</re:radiometricCorrectionApplied>[\s]*"
                           "<re:geoCorrectionLevel>([\S\s^<]*)</re:geoCorrectionLevel>[\s]*"
                           "<re:elevationCorrectionApplied>([\S]*)</re:elevationCorrectionApplied>[\s]*"
                           "<re:atmosphericCorrectionApplied>([\w]*)</re:atmosphericCorrectionApplied>", mxml_, re.I)
            self.additional.append(["Corrections", ["radiometricCorrectionApplied", h7.group(1)],
                                    ["geoCorrectionLevel", h7.group(2)], ["elevationCorrectionApplied", h7.group(3)],
748
                                    ["atmosphericCorrectionApplied", h7.group(4)]])
749
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
750
751
            self.additional.append(["Corrections", ["radiometricCorrectionApplied", h7.group(1)],
                                    ["radiometricCalibrationVersion", h7.group(2)], ["geoCorrectionLevel", h7.group(3)],
752
753
                                    ["elevationCorrectionApplied", h7.group(4)],
                                    ["atmosphericCorrectionApplied", h7.group(5)]])
Daniel Scheffler's avatar
Daniel Scheffler committed
754
755
            self.Quality.append(["geomProductAccuracy[m]:", str(
                round(float(h7.group(6)), 1))])  # Estimated product horizontal CE90 uncertainty [m]
756

757
        # additional. missing lines, suspectlines, binning, shifting, masking
758
759
760
761
762
763
        h81 = re.findall("<re:percentMissingLines>([^<]*)</re:percentMissingLines>", mxml_, re.I)
        h82 = re.findall("<re:percentSuspectLines>([^<]*)</re:percentSuspectLines>", mxml_, re.I)
        h83 = re.findall("<re:binning>([^<]*)</re:binning>", mxml_, re.I)
        h84 = re.findall("<re:shifting>([^<]*)</re:shifting>", mxml_, re.I)
        h85 = re.findall("<re:masking>([^<]*)</re:masking>", mxml_, re.I)

Daniel Scheffler's avatar
Daniel Scheffler committed
764
765
766
767
768
769
770
771
772
        self.Quality.append(
            ["MissingLines[%]perBand", h81])  # Percentage of missing lines in the source data of this band
        # Percentage of suspect lines (lines that contained downlink errors) in the source data for the band
        self.Quality.append(["SuspectLines[%]perBand", h82])
        self.additional.append(
            ["binning_perBand", h83])  # Indicates the binning used (across track x along track) [1x1,2x2,3x3,1x2,2x1]
        self.additional.append(
            ["shifting_perBand", h84])  # Indicates the sensor applied right shifting [none, 1bit, 2bits, 3bits, 4bits]
        self.additional.append(["masking_perBand", h85])  # Indicates the sensor applied masking [111, 110, 100, 000]
773

774
775
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
776
        self.PhysUnit = "DN"
777

778
        # Coordinate Reference System
Daniel Scheffler's avatar
Daniel Scheffler committed
779
780
781
782
783
784
785
786
787
788
789
790
        re_CS_EPSG = re.search('<re:ProductInformation>[\s\S]*<re:epsgCode>([0-9]*)</re:epsgCode>[\s\S]*'
                               '</re:ProductInformation>', mxml_, re.I)
        re_CS_TYPE = re.search('<re:ProductInformation>[\s\S]*<re:projection>([\s\S]*)</re:projection>[\s\S]*'
                               '</re:ProductInformation>', mxml_, re.I)
        re_CS_DATUM = re.search('<re:ProductInformation>[\s\S]*<re:geodeticDatum>([\w+\s]*)</re:geodeticDatum>[\s\S]*'
                                '</re:ProductInformation>', mxml_, re.I)
        re_CS_UTM_ZONE = re.search('<re:ProductInformation>[\s\S]*<re:projectionZone>([0-9]*)'
                                   '</re:projectionZone>[\s\S]*</re:ProductInformation>', mxml_, re.I)
        self.CS_EPSG = int(re_CS_EPSG.group(1)) if re_CS_EPSG else self.CS_EPSG
        self.CS_TYPE = 'LonLat' if re_CS_TYPE and not re.search('UTM', re_CS_TYPE.group(1)) else 'UTM' \
            if re_CS_TYPE and re.search('UTM', re_CS_TYPE.group(1)) else self.CS_TYPE
        self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS_1984' else self.CS_DATUM
791
792
793
        self.CS_UTM_ZONE = int(re_CS_UTM_ZONE.group(1)) if re_CS_UTM_ZONE is not None else self.CS_UTM_ZONE

        # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples)
Daniel Scheffler's avatar
Daniel Scheffler committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
        h10_UL = re.findall('<re:topLeft><re:latitude>(.*)</re:latitude><re:longitude>(.*)</re:longitude></re:topLeft>',
                            mxml_, re.I)
        h10_UR = re.findall(
            '<re:topRight><re:latitude>(.*)</re:latitude><re:longitude>(.*)</re:longitude></re:topRight>', mxml_, re.I)
        h10_LL = re.findall(
            '<re:bottomLeft><re:latitude>(.*)</re:latitude><re:longitude>(.*)</re:longitude></re:bottomLeft>', mxml_,
            re.I)
        h10_LR = re.findall(
            '<re:bottomRight><re:latitude>(.*)</re:latitude><re:longitude>(.*)</re:longitude></re:bottomRight>', mxml_,
            re.I)
        if h10_UL:  # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR)
            self.CornerTieP_LonLat.append(tuple([float(h10_UL[0][1]), float(h10_UL[0][0])]))
            self.CornerTieP_LonLat.append(tuple([float(h10_UR[0][1]), float(h10_UR[0][0])]))
            self.CornerTieP_LonLat.append(tuple([float(h10_LL[0][1]), float(h10_LL[0][0])]))
            self.CornerTieP_LonLat.append(tuple([float(h10_LR[0][1]), float(h10_LR[0][0])]))
809

810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
        ##########################
        # band specific metadata #
        ##########################

        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

        # Gains + Offsets
        h9 = re.findall("<re:radiometricScaleFactor>([^<]*)</re:radiometricScaleFactor>", mxml_, re.I)
        self.Gains = dict(zip(LBA_full_sorted, [float(gain) for gain in h9]))
        self.Offsets = dict(zip(LBA_full_sorted, [0, 0, 0, 0, 0]))

        # 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
        # noinspection PyBroadException
        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)
        except Exception:
            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)
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)

        # Exact values calculated based in SRFs
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
        # Provider values
        if not self.SolIrradiance:
            # Preconfigured Irradiation values
            self.SolIrradiance = dict(zip(LBA_full_sorted, [1997.8, 1863.5, 1560.4, 1395.0, 1124.4]))
            # Preconfigured CWLs
            self.CWL = dict(zip(LBA_full_sorted, [475., 555., 657., 710., 805.]))

843
844
845
846
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

847
848
849
850
851
852
853
854
855
    def Read_ASTER_hdffile(self, subsystem):
        """#----METHOD_5------------------------------------------------------------
        read metadata from ASTER hdf
        input:
            hdffile:
            subsystem:
        output:
        :param subsystem:
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
856

857
        #    self.default_attr()
Daniel Scheffler's avatar
Daniel Scheffler committed
858
        dat_ = open(self.FolderOrArchive, "r").read() if sys.version_info[0] < 3 else \
859
860
861
            open(self.FolderOrArchive, "rb").read().decode('latin-1')

        # Split meta from raster data
Daniel Scheffler's avatar
Daniel Scheffler committed
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
        meta = re.search("GROUP[\s]*=[\s]ASTERGENERICMETADATA[\s\S]*?END_GROUP[\s]*=[\s]INVENTORYMETADATA", dat_, re.I)
        meta_ = meta.group(0) if meta else ''
        genericmeta = re.search("GROUP[\s]*=[\s]ASTERGENERICMETADATA[\s\S]*?END_GROUP[\s]*=[\s]ASTERGENERICMETADATA",
                                meta_, re.I)
        inventorymeta = re.search("GROUP[\s]*=[\s]INVENTORYMETADATA[\s\S]*?END_GROUP[\s]*=[\s]INVENTORYMETADATA", meta_,
                                  re.I)
        gcsgenericmeta = re.search("GROUP[\s]*=[\s]GDSGENERICMETADATA[\s\S]*?END_GROUP[\s]*=[\s]GDSGENERICMETADATA",
                                   meta_, re.I)
        vnirmeta = re.search(
            "GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATAVNIR[\s\S]*?END_GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATAVNIR", meta_,
            re.I)
        swirmeta = re.search(
            "GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATASWIR[\s\S]*?END_GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATASWIR", meta_,
            re.I)
        tirmeta = re.search(
            "GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATATIR[\s\S]*?END_GROUP[\s]*=[\s]PRODUCTSPECIFICMETADATATIR", meta_,
            re.I)
        genericmeta_ = genericmeta.group(0) if genericmeta else ''
        inventorymeta_ = inventorymeta.group(0) if inventorymeta else ''
        gcsgenericmeta_ = gcsgenericmeta.group(0) if gcsgenericmeta else ''
        vnirmeta_ = vnirmeta.group(0) if vnirmeta else ''
        swirmeta_ = swirmeta.group(0) if swirmeta else ''
        tirmeta_ = tirmeta.group(0) if tirmeta else ''
        h_ = '\n\n\n'.join([genericmeta_, inventorymeta_, gcsgenericmeta_, vnirmeta_, swirmeta_, tirmeta_])
886
887
888
889
890

        # with open("./testing/out/ASTER_HDFmeta__h_.txt", "w") as allMetaOut:
        #     allMetaOut.write(h_)

        # EntityID
Daniel Scheffler's avatar
Daniel Scheffler committed
891
892
893
894
        h1 = re.search("OBJECT[\s]*=[\s]*IDOFASTERGDSDATAGRANULE[\s\S]*END_OBJECT[\s]*=[\s]*IDOFASTERGDSDATAGRANULE",
                       h_, re.I)
        h11 = re.search('VALUE[\s]*=[\s]*"([\s\S]*)"', h1.group(0), re.I)
        self.EntityID = [h11.group(1), re.search("\"(ASTL1A[\s0-9]*)\"", h_).group(1)]
895
896

        # Viewing Angle
Daniel Scheffler's avatar
Daniel Scheffler committed
897
898
899
900
        h2 = re.search("GROUP[\s]*=[\s]*POINTINGANGLES[\s\S]*END_GROUP[\s]*=[\s]*POINTINGANGLES", h_, re.I)
        h21 = re.findall("VALUE[\s]*=[\s]*([-]*[0-9.0-9]+)", h2.group(0), re.I)
        self.additional.append(["ViewingAngles", "VNIR", float(h21[0]), "SWIR", float(h21[1]), "TIR", float(h21[2])])
        if max(float(h21[0]), float(h21[1]), float(h21[2])) - min(float(h21[0]), float(h21[1]), float(h21[2])) < 1:
901
902
903
904
            self.ViewingAngle = float(h21[0])
        else:
            self.ViewingAngle = -99.0

905
        # additional GainMode
Daniel Scheffler's avatar
Daniel Scheffler committed
906
907
908
909
910
911
912
913
914
        h3 = re.search("GROUP[\s]*=[\s]*GAININFORMATION[\s\S]*END_GROUP[\s]*=[\s]*GAININFORMATION", genericmeta_, re.I)
        h31 = re.findall('VALUE[\s]*=[\s]*[\S]?\"[0-9A-Z]*\", \"([A-Z]*)\"', h3.group(0), re.I)
        gains = {'HGH': [170.8, 179.0, 106.8, 27.5, 8.8, 7.9, 7.55, 5.27, 4.02, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'],
                 'NOR': [427.0, 358.0, 218.0, 55.0, 17.6, 15.8, 15.1, 10.55, 8.04, 28.17, 27.75, 26.97, 23.30, 21.38],
                 'LOW': [569.0, 477.0, 290.0, 73.3, 23.4, 21.0, 20.1, 14.06, 10.72, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'],
                 'LO1': [569.0, 477.0, 290.0, 73.3, 23.4, 21.0, 20.1, 14.06, 10.72, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'],
                 'LO2': ['N/A', 'N/A', 'N/A', 73.3, 103.5, 98.7, 83.8, 62.0, 67.0, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'],
                 'OFF': "OFF"}
        self.additional.append([["GainMode:"], ["Max_radiance:"]])
915
916
917
918
        for x, i in enumerate(h31[:15]):
            self.additional[-1][-2].append(i)
            self.additional[-1][-1].append(gains[i][x])

919
920
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
921
        self.PhysUnit = "DN"
922
923

        # Satellite
Daniel Scheffler's avatar
Daniel Scheffler committed
924
        self.Satellite = 'Terra'
925
926

        # Sensor
Daniel Scheffler's avatar
Daniel Scheffler committed
927
928
        h10 = re.search("OBJECT[\s]*=[\s]*INSTRUMENTSHORTNAME[\s\S]*END_OBJECT[\s]*=[\s]*INSTRUMENTSHORTNAME", h_, re.I)
        self.Sensor = re.search("VALUE[\s]*=[\s]*[\"]([A-Za-z]*)\"", h10.group(0), re.I).group(1)
929
930

        # Subsystem
Daniel Scheffler's avatar
Daniel Scheffler committed
931
932
933
        h5 = re.search("GROUP[\s]*=[\s]*OBSERVATIONMODE[\s\S]*END_GROUP[\s]*=[\s]*OBSERVATIONMODE", h_, re.I)
        avail_subsystems = re.findall('VALUE[\s]*=[\s]*[(]\"([a-zA-Z0-9]*)\", \"([ONF]*)\"', h5.group(0), re.I)
        if subsystem not in ['VNIR1', 'VNIR2', 'SWIR', 'TIR']:
934
935
            raise ValueError('Unexpected subsystem >%s<. Unable to specify the correct ASTER subsystem to be '
                             'processed.' % subsystem)
936
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
937
938
939
940
941
942
943
944
            if subsystem == 'VNIR1' and avail_subsystems[0][1] == 'ON':
                self.nBands = 3
            if subsystem == 'VNIR2' and avail_subsystems[1][1] == 'ON':
                self.nBands = 1
            if subsystem == 'SWIR' and avail_subsystems[2][1] == 'ON':
                self.nBands = 6
            if subsystem == 'TIR' and avail_subsystems[3][1] == 'ON':
                self.nBands = 5
945
946
947
        self.Subsystem = subsystem

        # Field of view (requires Satellite, Sensor, Subsystem)
948
        self.FOV = get_FieldOfView(self.GMS_identifier)
949
950

        # Ground resolution
Daniel Scheffler's avatar
Daniel Scheffler committed
951
952
953
954
955
956
        re_res_GSD = re.findall('OBJECT[\s]*=[\s]*SPATIALRESOLUTION[\s\S]*VALUE[\s]*=[\s]*[\S]{1}([1-9][0-9]), '
                                '([1-9][0-9]), ([1-9][0-9])[\s\S]*END_OBJECT[\s]*=[\s]*SPATIALRESOLUTION', h_, re.I)
        idx_subsystem = \
            0 if subsystem[:4] == 'VNIR' else \
            1 if subsystem[:4] == 'SWIR' else \
            2 if subsystem[:4] == 'TIR' else None
957
        self.gResolution = float(re_res_GSD[0][idx_subsystem]) \
958
            if re_res_GSD and idx_subsystem is not None else self.gResolution
959
960

        # Flight direction
Daniel Scheffler's avatar
Daniel Scheffler committed
961
962
        h6 = re.search("OBJECT[\s]*=[\s]*FLYINGDIRECTION[\s\S]*\"([ASDE]*?)\"", h_, re.I)
        Fdir = {'AS': "Ascending", 'DE': "Descending"}
963
964
965
        self.additional.append(["Flight Direction", Fdir[h6.group(1)]])

        # SunAzimuth SunElevation
Daniel Scheffler's avatar
Daniel Scheffler committed
966
967
968
        h7 = re.search("OBJECT[\s]*=[\s]*SOLARDIRECTION[\s\S]*END_OBJECT[\s]*=[\s]*SOLARDIRECTION", h_, re.I)
        h71 = re.search("VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)", h7.group(0), re.I)
        self.SunAzimuth = float(h71.group(1))
969
970
971
        self.SunElevation = float(h71.group(2))

        # data Quality
Daniel Scheffler's avatar
Daniel Scheffler committed
972
973
974
        h8 = re.findall(r"GROUP[\s]*=[\s]*DATAQUALITY[1-9][0-4BN]?[^.]*END_GROUP[\s]*=[\s]*DATAQUALITY[1-9][0-4BN]?",
                        h_, re.I)
        h81 = re.findall('VALUE[\s]*=[\s]*([^\n]*)', "\n".join(h8), re.I)
975
976
977
978
        self.Quality.append(["NoOfBadPixel:"])
        for i in h81:
            self.Quality[-1].append(i)