metadata.py 114 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
17
18

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

Daniel Scheffler's avatar
Daniel Scheffler committed
22
23
from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.geo.projection import WKT2EPSG
24

25
from gms_preprocessing.options.config import GMS_config as CFG
26
27
28
from gms_preprocessing.io.input_reader import open_specific_file_within_archive, Solar_Irradiance_reader, SRF_reader
from gms_preprocessing.io.output_writer import enviHdr_keyOrder
from gms_preprocessing.algorithms import geoprocessing as GEOP
29
from gms_preprocessing.misc import helper_functions as HLP_F
Daniel Scheffler's avatar
Daniel Scheffler committed
30
from gms_preprocessing.misc import database_tools as DB_T
31
from gms_preprocessing.misc.path_generator import path_generator, get_path_ac_options
32
33
34
35
from gms_preprocessing.misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene, \
    datasetid_to_sat_sen
from gms_preprocessing.misc.exceptions import ACNotSupportedError

36

37
from sicor.options import get_options as get_ac_options
38

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

41

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

47
48
49
        # unpack GMS_identifier
        self.GMS_identifier = GMS_identifier
        self.image_type = GMS_identifier['image_type']
Daniel Scheffler's avatar
Daniel Scheffler committed
50
51
52
        self.Satellite = GMS_identifier['Satellite']
        self.Sensor = GMS_identifier['Sensor']
        self.Subsystem = GMS_identifier['Subsystem']
53
        self.proc_level = GMS_identifier['proc_level']
Daniel Scheffler's avatar
Daniel Scheffler committed
54
        self.logger = GMS_identifier['logger']
55
56
57

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

119
120
121
122
123
124
125
126
127
128
129
130
131
132
    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)
133
134
        elif re.search("Sentinel-2A", self.Satellite, re.I) or re.search("Sentinel-2B", self.Satellite, re.I):
            self.Read_Sentinel2_xmls()
135
136
137
138
139
140
141
142
        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:
143
            raise NotImplementedError("%s is not a supported sensor or sensorname is misspelled." % self.Satellite)
144
145

        return self
146

147
148
149
    @property
    def AcqDateTime(self):
        """Returns a datetime.datetime object containing date, time and timezone (UTC time)."""
Daniel Scheffler's avatar
Daniel Scheffler committed
150

151
152
153
154
155
156
        return self._AcqDateTime

    @AcqDateTime.setter
    def AcqDateTime(self, DateTime):
        # type: (datetime.datetime) -> None
        self._AcqDateTime = DateTime
Daniel Scheffler's avatar
Daniel Scheffler committed
157
158
        self.AcqDate = datetime.datetime.strftime(DateTime, format='%Y-%m-%d')
        self.AcqTime = datetime.datetime.strftime(DateTime, format='%H:%M:%S')
159

160
161
162
    @property
    def overview(self):
        attr2include = \
Daniel Scheffler's avatar
Daniel Scheffler committed
163
164
            ['Dataname', 'Metafile', 'EntityID', 'Satellite', 'Sensor', 'Subsystem', 'gResolution', 'AcqDate',
             'AcqTime',
165
166
             'CWL', 'FWHM', 'Offsets', 'Gains', 'ProcLCode', 'ProcLStr', 'SunElevation', 'SunAzimuth',
             'ViewingAngle', 'IncidenceAngle', 'FOV', 'SolIrradiance', 'ThermalConstK1', 'ThermalConstK2',
Daniel Scheffler's avatar
Daniel Scheffler committed
167
168
             'EarthSunDist', 'Quality', 'PhysUnit', 'additional', 'GainsRef', 'OffsetsRef', 'CornerTieP_LonLat',
             'CS_EPSG',
169
170
171
             'CS_TYPE', 'CS_DATUM', 'CS_UTM_ZONE', 'LayerBandsAssignment']
        return collections.OrderedDict((attr, getattr(self, attr)) for attr in attr2include)

172
173
174
175
    def Read_SPOT_dimap2(self):
        """----METHOD_2------------------------------------------------------------
        # read metadata from spot dimap file
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
176

Daniel Scheffler's avatar
Daniel Scheffler committed
177
        # self.default_attr()
178
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
179
180
            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
181
182
            self.Metafile = glob_res[0]
            dim_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
183
        else:  # archive
184
            dim_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/scene01/metadata.dim')
185

Daniel Scheffler's avatar
Daniel Scheffler committed
186
187
188
        # 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)
189
        self.additional.append(["SpecialValues:"])
Daniel Scheffler's avatar
Daniel Scheffler committed
190
191
        for ii, ind in enumerate(h1):
            self.additional[0].append(["%s:%s" % (ind, h2[ii])])
192

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

Daniel Scheffler's avatar
Daniel Scheffler committed
197
198
        # AcqDate
        h4 = re.search("<IMAGING_DATE>([0-9-]*)</IMAGING_DATE>", dim_, re.I)
199
        AcqDate = h4.group(1)
200

Daniel Scheffler's avatar
Daniel Scheffler committed
201
        # Earth sun distance
202
        self.EarthSunDist = self.get_EarthSunDistance(AcqDate)
203

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
212
        # Satellite + Sensor
213
214
        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
215
216
217
                       "</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)])
218

Daniel Scheffler's avatar
Daniel Scheffler committed
219
        # Angles: incidence angle, sunAzimuth, sunElevation
220
        h7 = re.search("<INCIDENCE_ANGLE>(.*)</INCIDENCE_ANGLE>[\s\S]*<SUN_AZIMUTH>(.*)</SUN_AZIMUTH>[\s\S]"
Daniel Scheffler's avatar
Daniel Scheffler committed
221
                       "*<SUN_ELEVATION>(.*)</SUN_ELEVATION>", dim_, re.I)
222
223
224
225
        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
226
227
        # Viewing Angle: not always included in the metadata.dim file
        h8 = re.search("<VIEWING_ANGLE>(.*)</VIEWING_ANGLE>", dim_, re.I)
228
        if type(h8).__name__ == 'NoneType':
Daniel Scheffler's avatar
Daniel Scheffler committed
229
            self.ViewingAngle = 0
230
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
231
            self.ViewingAngle = float(h8.group(1))
232

Daniel Scheffler's avatar
Daniel Scheffler committed
233
        # Field of View
234
        self.FOV = get_FieldOfView(self.GMS_identifier)
235

Daniel Scheffler's avatar
Daniel Scheffler committed
236
        # Units
237
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
238
        self.PhysUnit = "DN"
239

Daniel Scheffler's avatar
Daniel Scheffler committed
240
241
        # Gains and Offsets
        h9 = re.search("<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I)
242
243
244
245
246
        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:"])
Daniel Scheffler's avatar
Daniel Scheffler committed
247
        for ii, ind in enumerate(h91):
248
            self.additional[-1].append(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
249
250
        # Offsets
        for ii, ind in enumerate(h92):
251
            self.Offsets.append(float(ind))
Daniel Scheffler's avatar
Daniel Scheffler committed
252
253
254
255
256
257
258
259
        # Gains
        for ii, ind in enumerate(h93):
            # gains in dimap file represent reciprocal 1/gain
            self.Gains.append(1 / float(ind))

        # 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
260
        # noinspection PyBroadException
Daniel Scheffler's avatar
Daniel Scheffler committed
261
262
263
264
        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)
265
        except Exception:
Daniel Scheffler's avatar
Daniel Scheffler committed
266
267
            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)
268
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)
269

Daniel Scheffler's avatar
Daniel Scheffler committed
270
        # Exact values calculated based in SRFs
271
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
272
        # Provider values
273
        if not self.SolIrradiance:
Daniel Scheffler's avatar
Daniel Scheffler committed
274
            h10 = re.search("<Solar_Irradiance>[\s\S]*</Solar_Irradiance>", dim_, re.I)
275
276
277
278
            h10_ = h10.group(0)
            h101 = re.findall("<SOLAR_IRRADIANCE_VALUE>([^<]*)</SOLAR_IRRADIANCE_VALUE>", h10_, re.I)
            if h101:
                self.SolIrradiance = h101
Daniel Scheffler's avatar
Daniel Scheffler committed
279
280
281
282
                #    self.additional.append(["Solar Irradiance per band:"])
                #    for ii,ind in enumerate(h101):
                #       self.additional[-1].append(ind)
            else:  # Preconfigured Irradiation values
283
                spot_irr_dic = {
Daniel Scheffler's avatar
Daniel Scheffler committed
284
285
286
287
288
                    'SPOT1a': [1855., 1615., 1090., 1680.], 'SPOT1b': [1845., 1575., 1040., 1690.],
                    'SPOT2a': [1865., 1620., 1085., 1705.], 'SPOT2b': [1865., 1615., 1090., 1670.],
                    'SPOT3a': [1854., 1580., 1065., 1668.], 'SPOT3b': [1855., 1597., 1067., 1667.],
                    'SPOT4a': [1843., 1568., 1052., 233., 1568.], 'SPOT4b': [1851., 1586., 1054., 240., 1586.],
                    'SPOT5a': [1858., 1573., 1043., 236., 1762.], 'SPOT5b': [1858., 1575., 1047., 234., 1773.]}
289
                self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)]
Daniel Scheffler's avatar
Daniel Scheffler committed
290
            # Preconfigured CWLs --   # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined
291
292
            # 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
Daniel Scheffler's avatar
Daniel Scheffler committed
293
294
295
            if get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b',
                                                               'SPOT3a', 'SPOT3b']:
                self.CWL = [545., 638., 819., 615.]
296
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT4a', 'SPOT4b']:
Daniel Scheffler's avatar
Daniel Scheffler committed
297
                self.CWL = [540., 650., 835., 1630., 645.]
298
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] == ['SPOT5a', 'SPOT5b']:
Daniel Scheffler's avatar
Daniel Scheffler committed
299
                self.CWL = [540., 650., 835., 1630., 595.]
300

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

Daniel Scheffler's avatar
Daniel Scheffler committed
305
        # Quality
306
        h12 = re.findall("<Bad_Pixel>[\s]*<PIXEL_INDEX>([0-9]*)</PIXEL_INDEX>[\s]*<BAD_PIXEL_STATUS>([^<]*)"
Daniel Scheffler's avatar
Daniel Scheffler committed
307
                         "</BAD_PIXEL_STATUS></Bad_Pixel>", dim_, re.I)
308
309
        self.Quality = h12

Daniel Scheffler's avatar
Daniel Scheffler committed
310
311
312
313
314
315
316
317
318
319
320
321
        # 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)
322
323
        if len(h121) == 4 and self.CS_TYPE == 'LonLat':
            # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR)
Daniel Scheffler's avatar
Daniel Scheffler committed
324
325
326
327
328
329
            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
330

331
332
333
334
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

335
336
337
338
339
    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
340

341
342
        # self.default_attr()
        self.LayerBandsAssignment = LayerBandsAssignment
Daniel Scheffler's avatar
Daniel Scheffler committed
343
        self.nBands = len(LayerBandsAssignment)
344
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
345
346
            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
347
348
            self.Metafile = glob_res[0]
            mtl_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
349
        else:  # archive
350
            mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt')
351
352

            # Processing Level
Daniel Scheffler's avatar
Daniel Scheffler committed
353
        h1 = re.search('PRODUCT_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
354
        if h1 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
355
            h1 = re.search('DATA_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I)
356
357
358
        self.ProcLCode = h1.group(1)

        # Satellite + Sensor + Sensor Mode
Daniel Scheffler's avatar
Daniel Scheffler committed
359
360
        h2 = re.search('SPACECRAFT_ID = "([a-zA-Z0-9_]*)"[\s]*SENSOR_ID = "([a-zA-Z0-9+]*)"[\s]*'
                       'SENSOR_MODE = "([\S]*)"', mtl_, re.I)
361
        if h2:
Daniel Scheffler's avatar
Daniel Scheffler committed
362
363
            self.Satellite = 'Landsat-%s' % re.search('LANDSAT[\D]*([0-9])', h2.group(1), re.I).group(1)
            self.Sensor = h2.group(2)
364
365
            self.Sensormode = h2.group(3)
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
366
367
368
369
370
371
372
            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
373
374

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

379
        # Acquisition Date + Time
Daniel Scheffler's avatar
Daniel Scheffler committed
380
        h3 = re.search('ACQUISITION_DATE = ([0-9-]*)[\s]*SCENE_CENTER_SCAN_TIME = "?([0-9:]*)"?', mtl_, re.I)
381
        if h3 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
382
            h3 = re.search('DATE_ACQUIRED = ([0-9-]*)[\s]*SCENE_CENTER_TIME = "?([0-9:]*)"?', mtl_, re.I)
383
384
385
386
387
388
        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')
389
390
391
392

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

393
394
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
395
        self.PhysUnit = "DN"
396

397
        # Gains and Offsets
Daniel Scheffler's avatar
Daniel Scheffler committed
398
        h4 = re.search("GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
399
        h4_ = h4.group(0)
Daniel Scheffler's avatar
Daniel Scheffler committed
400
401
        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
402
403
        h4max_pix = re.findall("QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall("QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
404
        if not h4max_rad:
Daniel Scheffler's avatar
Daniel Scheffler committed
405
406
407
408
            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
409
410
411
            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
Daniel Scheffler's avatar
Daniel Scheffler committed
412
413
414
415
416
417
        self.additional.append([["LMAX"], ["LMIN"], ["QCALMAX"], ["QCALMIN"]])
        # Offsets + Gains
        for ii, ind in enumerate(h4min_rad):
            self.Gains.append(
                (float(h4max_rad[ii]) - float(h4min_rad[ii])) / (float(h4max_pix[ii]) - float(h4min_pix[ii])))
            self.Offsets.append(float(ind) - float(h4min_pix[ii]) * self.Gains[-1])
418
419
420
421
422
423
            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])

            # Solar irradiance, central wavelengths , full width half maximum per band
Daniel Scheffler's avatar
Daniel Scheffler committed
424
425
        self.wvlUnit = 'Nanometers'
        # Exact values calculated based in SRFs
426
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
427
428
429
430
431
432
433
434
435
436
437
438
        # Provider values
        if not self.SolIrradiance:  # Preconfigured Irradiation and CWL values
            if re.search("[\D]*5", self.Satellite):  # Landsat 5; resource for center wavelength:
                self.SolIrradiance = [1957., 1826., 1554., 1036., 215., 0.0, 80.67]  # B1,B2,B3,B4,B5,B6(thermal),B7
                self.CWL = [485., 560., 660., 830., 1650., 11450., 2215.]
            if re.search("[\D]*7", self.Satellite):
                # Landsat 7; resource for center wavelength:
                # http://opticks.org/confluence/display/opticksDev/Sensor+Wavelength+Definitions
                self.SolIrradiance = [1969., 1840., 1551., 1044., 225.7, 0.0, 0.0, 82.07,
                                      1368.]  # B1,B2,B3,B4,B5,B61(thermal),B62(thermal),B7,B8(pan)
                self.CWL = [482.5, 565., 660., 825., 1650., 11450., 11450., 2215., 710.]
            if re.search("[\D]*8", self.Satellite):  # Landsat 8
439
                # no irradiation values available
Daniel Scheffler's avatar
Daniel Scheffler committed
440
441
442
443
444
445
446
447
448
449
450
451
452
                self.CWL = [443., 482.6, 561.3, 654., 864.6, 1609., 2201., 1307.5, 10895., 12005.]
                #    if None in self.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
            self.ThermalConstK1 = [0., 0., 0., 0., 0., 607.76, 0.]
            self.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
            self.ThermalConstK1 = [0., 0., 0., 0., 0., 666.09, 666.09, 0., 0.]
            self.ThermalConstK2 = [0., 0., 0., 0., 0., 1282.71, 1282.71, 0., 0.]
        if re.search("[\D]*8", self.Satellite):  # Landsat 8
453
454
455
            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:
Daniel Scheffler's avatar
Daniel Scheffler committed
456
457
                self.ThermalConstK1 = [0., 0., 0., 0., 0., 0., 0., 0., float(K1_res[0]), float(K1_res[1])]
                self.ThermalConstK2 = [0., 0., 0., 0., 0., 0., 0., 0., float(K2_res[0]), float(K2_res[1])]
458
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
459
460
                self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. '
                                  'Found %s' % len(K1_res))
461
462

                # Angles: incidence angle, sunAzimuth, sunElevation, field of view
Daniel Scheffler's avatar
Daniel Scheffler committed
463
464
        h5 = re.search("SUN_AZIMUTH = ([\S]*)[\s]*SUN_ELEVATION = ([\S]*)", mtl_, re.I)
        self.SunAzimuth = float(h5.group(1))
465
        self.SunElevation = float(h5.group(2))
Daniel Scheffler's avatar
Daniel Scheffler committed
466
        self.FOV = get_FieldOfView(self.GMS_identifier)
467
468

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

Daniel Scheffler's avatar
Daniel Scheffler committed
473
        h6_ = h6.group(0)
474
        h61 = (h6_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
475
        x = 0
476
        for i in h61:
Daniel Scheffler's avatar
Daniel Scheffler committed
477
            if x == 0 or x + 1 == len(h61):
478
479
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
480
                i_ = i.strip().replace('"', "")
481
                self.Quality.append(i_.split(" = "))
Daniel Scheffler's avatar
Daniel Scheffler committed
482
            x += 1
483
484

            # Additional: coordinate system, geom. Resolution
Daniel Scheffler's avatar
Daniel Scheffler committed
485
486
        h7 = re.search("GROUP = PROJECTION_PARAMETERS[\s\S]*END_GROUP = L1_METADATA_FILE", mtl_, re.I)
        h7_ = h7.group(0)
487
        h71 = (h7_.split("\n"))
Daniel Scheffler's avatar
Daniel Scheffler committed
488
489
        for x, i in enumerate(h71):
            if re.search("Group", i, re.I):
490
491
                pass
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
492
                i_ = i.strip().replace('"', "")
493
                self.additional.append(i_.split(" = "))
Daniel Scheffler's avatar
Daniel Scheffler committed
494
495
496
497
498
499
500
501
        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
502
503
504
        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
505
506
        # Landsat8
        h8 = re.search("ROLL_ANGLE = ([\S]*)", mtl_, re.I)
507
        if h8:
508
509
510
            self.ViewingAngle = float(h8.group(1))

            # reflectance coefficients (Landsat8)
Daniel Scheffler's avatar
Daniel Scheffler committed
511
        h9 = re.search("GROUP = RADIOMETRIC_RESCALING[\s\S]*END_GROUP = RADIOMETRIC_RESCALING", mtl_, re.I)
512
        if h9:
513
            h9_ = h9.group(0)
Daniel Scheffler's avatar
Daniel Scheffler committed
514
            h9gain_ref = re.findall(" REFLECTANCE_MULT_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
515
            h9gain_ref_bandNr = re.findall(" REFLECTANCE_MULT_BAND_([0-9]+) = [\S]*", h9_, re.I)
Daniel Scheffler's avatar
Daniel Scheffler committed
516
517
            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)
518
            if h9gain_ref:
Daniel Scheffler's avatar
Daniel Scheffler committed
519
520
                self.GainsRef = [None] * len(self.Gains)
                self.OffsetsRef = [None] * len(self.Offsets)
521
                FullLBA = get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False)
522

Daniel Scheffler's avatar
Daniel Scheffler committed
523
                for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref):
524
                    self.GainsRef[FullLBA.index(ii)] = ind
Daniel Scheffler's avatar
Daniel Scheffler committed
525
                for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref):
526
527
528
529
530
531
532
                    self.OffsetsRef[FullLBA.index(ii)] = ind

        # 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
533
        if not h10_UL:  # Landsat8
534
535
536
537
            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
538
539
540
541
542
543
544
545
546
        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]]))
547
548

        # WRS path/row
Daniel Scheffler's avatar
Daniel Scheffler committed
549
        h11_p = re.search('WRS_PATH = ([0-9]*)', mtl_, re.I)
550
        if h11_p:
551
            self.WRS_path = h11_p.group(1)
Daniel Scheffler's avatar
Daniel Scheffler committed
552
553
554
555
        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)
556
            self.WRS_row = int(h11_r.group(1))
557
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
558
            self.WRS_row = int(h11_r1.group(1)) if h11_r1.group(1) == h11_r2.group(1) else self.WRS_row
559
560

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

Daniel Scheffler's avatar
Daniel Scheffler committed
566
            if len(result) == 1:  # e.g. [('LE71950282003121EDC00',)]
567
568
                self.EntityID = result[0][0]
            elif len(result) == 0:
569
570
                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
571
            else:  # e.g. [('LE71950282003121EDC00',), ('LE71950282003105ASN00',)]
572
                self.logger.warning("Scene ID could not be assigned because %s datasets matching to the query "
Daniel Scheffler's avatar
Daniel Scheffler committed
573
                                    "parameters were found in metadata database." % len(result))
574
        # if self.EntityID=='':
Daniel Scheffler's avatar
Daniel Scheffler committed
575
576
        #     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 ''
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
        #     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']

601
602
603
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)
604
605
606
607
608
609
610

        # mtl.close()

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

612
613
        # self.default_attr()
        if os.path.isdir(self.FolderOrArchive):
Daniel Scheffler's avatar
Daniel Scheffler committed
614
615
            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
616
617
            self.Metafile = glob_res[0]
            mxml_ = open(self.Metafile, "r").read()
Daniel Scheffler's avatar
Daniel Scheffler committed
618
        else:  # archive
619
            mxml_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*_metadata.xml')
620

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

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

629
        # Satellite
Daniel Scheffler's avatar
Daniel Scheffler committed
630
631
        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"
632

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

637
        # Acquisition Parameters: Incidence Angle, SunAzimuth, SunElevation, ViewingAngle, FieldOfView, AcqDate, AcqTime
638
639
640
641
642
643
        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
644
                           'Time>([0-9-]*)T([0-9:]*)[\S]*</re:acquisitionDateTime>', mxml_, re.I)
645
            self.IncidenceAngle = float(h5.group(1))
Daniel Scheffler's avatar
Daniel Scheffler committed
646
            self.SunAzimuth = float(h5.group(2))
647
            self.SunElevation = float(h5.group(3))
Daniel Scheffler's avatar
Daniel Scheffler committed
648
649
650
            # 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))])
651
            self.ViewingAngle = float(h5.group(5))
Daniel Scheffler's avatar
Daniel Scheffler committed
652
            self.FOV = get_FieldOfView(self.GMS_identifier)
653
654
            AcqDate = h5.group(6)
            AcqTime = h5.group(7)
655
656

        except AttributeError:
Daniel Scheffler's avatar
Daniel Scheffler committed
657
658
659
660
661
            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)  #
662
663
664
            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
665
            self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))])
666
            self.ViewingAngle = 9999.99
Daniel Scheffler's avatar
Daniel Scheffler committed
667
            self.FOV = get_FieldOfView(self.GMS_identifier)
668
669
670
671
672
            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')
673

674
        # Earth sun distance
675
676
677
        self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate)

        # Additional: Projection
Daniel Scheffler's avatar
Daniel Scheffler committed
678
679
        h6 = re.search("<re:geodeticDatum>([\S]*)</re:geodeticDatum>[\s]*<re:projection>([\S^<\s]*)</re:projection>",
                       mxml_, re.I)
680
        try:
Daniel Scheffler's avatar
Daniel Scheffler committed
681
682
            self.additional.append(["SpatialReference", ["geodeticDatum", h6.group(1)], ["projection", h6.group(2)]])
        except AttributeError:  # Level1-B data has no geodaticDatum
683
684
            pass

685
        # Corrections applied + Quality
Daniel Scheffler's avatar
Daniel Scheffler committed
686
687
688
689
690
691
        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)
692
693
        # fuer IRIS ihre Daten
        if h7 is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
694
695
696
697
698
699
            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)],
700
                                    ["atmosphericCorrectionApplied", h7.group(4)]])
701
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
702
703
            self.additional.append(["Corrections", ["radiometricCorrectionApplied", h7.group(1)],
                                    ["radiometricCalibrationVersion", h7.group(2)], ["geoCorrectionLevel", h7.group(3)],
704
705
                                    ["elevationCorrectionApplied", h7.group(4)],
                                    ["atmosphericCorrectionApplied", h7.group(5)]])
Daniel Scheffler's avatar
Daniel Scheffler committed
706
707
            self.Quality.append(["geomProductAccuracy[m]:", str(
                round(float(h7.group(6)), 1))])  # Estimated product horizontal CE90 uncertainty [m]
708

709
        # additional. missing lines, suspectlines, binning, shifting, masking
710
711
712
713
714
715
        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
716
717
718
719
720
721
722
723
724
        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]
725

726
727
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
728
        self.PhysUnit = "DN"
729

730
731
732
733
        # Gains + Offsets
        h9 = re.findall("<re:radiometricScaleFactor>([^<]*)</re:radiometricScaleFactor>", mxml_, re.I)
        for gain in h9:
            self.Gains.append(float(gain))
Daniel Scheffler's avatar
Daniel Scheffler committed
734
        self.Offsets = [0, 0, 0, 0, 0]
735
736

        # Solar irradiance, central wavelengths , full width half maximum per band
Daniel Scheffler's avatar
Daniel Scheffler committed
737
738
        self.wvlUnit = 'Nanometers'
        # derive number of bands from number of given gains/offsets in header file or from raster data
739
        # noinspection PyBroadException
Daniel Scheffler's avatar
Daniel Scheffler committed
740
741
742
743
        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)
744
        except Exception:
Daniel Scheffler's avatar
Daniel Scheffler committed
745
746
            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)
747
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)
748

Daniel Scheffler's avatar
Daniel Scheffler committed
749
        # Exact values calculated based in SRFs
750
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
751
        # Provider values
752
        if not self.SolIrradiance:
Daniel Scheffler's avatar
Daniel Scheffler committed
753
754
755
756
            # Preconfigured Irradiation values
            self.SolIrradiance = [1997.8, 1863.5, 1560.4, 1395.0, 1124.4]
            # Preconfigured CWLs
            self.CWL = [475., 555., 657., 710., 805.]
757
758

            # Coordinate Reference System
Daniel Scheffler's avatar
Daniel Scheffler committed
759
760
761
762
763
764
765
766
767
768
769
770
        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
771
772
773
        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
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
        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])]))
789

790
791
792
793
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

794
795
796
797
798
799
800
801
802
    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
803

804
        #    self.default_attr()
Daniel Scheffler's avatar
Daniel Scheffler committed
805
        dat_ = open(self.FolderOrArchive, "r").read() if sys.version_info[0] < 3 else \
806
807
808
            open(self.FolderOrArchive, "rb").read().decode('latin-1')

        # Split meta from raster data
Daniel Scheffler's avatar
Daniel Scheffler committed
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
        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_])
833
834
835
836
837

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

        # EntityID
Daniel Scheffler's avatar
Daniel Scheffler committed
838
839
840
841
        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)]
842
843

        # Viewing Angle
Daniel Scheffler's avatar
Daniel Scheffler committed
844
845
846
847
        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:
848
849
850
851
            self.ViewingAngle = float(h21[0])
        else:
            self.ViewingAngle = -99.0

852
        # additional GainMode
Daniel Scheffler's avatar
Daniel Scheffler committed
853
854
855
856
857
858
859
860
861
        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:"]])
862
863
864
865
        for x, i in enumerate(h31[:15]):
            self.additional[-1][-2].append(i)
            self.additional[-1][-1].append(gains[i][x])

866
867
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
868
        self.PhysUnit = "DN"
869
870

        # Gains/Offsets
Daniel Scheffler's avatar
Daniel Scheffler committed
871
872
873
874
        h4 = re.findall(r"GROUP[\s]*=[\s]*UNITCONVERSIONCOEFF[1-9][0-4BN]?[^(]*END_GROUP[\s]*=[\s]*"
                        r"UNITCONVERSIONCOEFF[1-9][0-4BN]?", h_, re.I)
        h41 = re.findall('VALUE[\s]*=[\s]*([0-9.0-9]+)', "\n".join(h4), re.I)
        h42 = re.findall('VALUE[\s]*=[\s]*([-]+[0-9.0-9]+)', "\n".join(h4), re.I)
875
876
877
878
        for x, i in