metadata.py 112 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
from gms_preprocessing.misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene
33

34
from sicor.options import get_options as get_ac_options
35

36
37
__author__ = 'Daniel Scheffler', 'Robert Behling'

38

39
class METADATA(object):
40
    def __init__(self, GMS_identifier):
41
42
43
        # private attributes
        self._AcqDateTime = None

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

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

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

        return self
143

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

148
149
150
151
152
153
        return self._AcqDateTime

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

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

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

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

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

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
198
        # Earth sun distance
199
        self.EarthSunDist = self.get_EarthSunDistance(AcqDate)
200

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

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

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
230
        # Field of View
231
        self.FOV = get_FieldOfView(self.GMS_identifier)
232

Daniel Scheffler's avatar
Daniel Scheffler committed
233
        # Units
234
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
235
        self.PhysUnit = "DN"
236

Daniel Scheffler's avatar
Daniel Scheffler committed
237
238
        # Gains and Offsets
        h9 = re.search("<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I)
239
240
241
242
243
        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
244
        for ii, ind in enumerate(h91):
245
            self.additional[-1].append(ind)
Daniel Scheffler's avatar
Daniel Scheffler committed
246
247
        # Offsets
        for ii, ind in enumerate(h92):
248
            self.Offsets.append(float(ind))
Daniel Scheffler's avatar
Daniel Scheffler committed
249
250
251
252
253
254
255
256
        # 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
257
        # noinspection PyBroadException
Daniel Scheffler's avatar
Daniel Scheffler committed
258
259
260
261
        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)
262
        except Exception:
Daniel Scheffler's avatar
Daniel Scheffler committed
263
264
            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)
265
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)
266

Daniel Scheffler's avatar
Daniel Scheffler committed
267
        # Exact values calculated based in SRFs
268
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
269
        # Provider values
270
        if not self.SolIrradiance:
Daniel Scheffler's avatar
Daniel Scheffler committed
271
            h10 = re.search("<Solar_Irradiance>[\s\S]*</Solar_Irradiance>", dim_, re.I)
272
273
274
275
            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
276
277
278
279
                #    self.additional.append(["Solar Irradiance per band:"])
                #    for ii,ind in enumerate(h101):
                #       self.additional[-1].append(ind)
            else:  # Preconfigured Irradiation values
280
                spot_irr_dic = {
Daniel Scheffler's avatar
Daniel Scheffler committed
281
282
283
284
285
                    '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.]}
286
                self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)]
Daniel Scheffler's avatar
Daniel Scheffler committed
287
            # Preconfigured CWLs --   # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined
288
289
            # 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
290
291
292
            if get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b',
                                                               'SPOT3a', 'SPOT3b']:
                self.CWL = [545., 638., 819., 615.]
293
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT4a', 'SPOT4b']:
Daniel Scheffler's avatar
Daniel Scheffler committed
294
                self.CWL = [540., 650., 835., 1630., 645.]
295
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] == ['SPOT5a', 'SPOT5b']:
Daniel Scheffler's avatar
Daniel Scheffler committed
296
                self.CWL = [540., 650., 835., 1630., 595.]
297

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
307
308
309
310
311
312
313
314
315
316
317
318
        # 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)
319
320
        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
321
322
323
324
325
326
            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
327

328
329
330
331
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

332
333
334
335
336
    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
337

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

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

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

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

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

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

390
391
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
392
        self.PhysUnit = "DN"
393

394
        # Gains and Offsets
Daniel Scheffler's avatar
Daniel Scheffler committed
395
        h4 = re.search("GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
396
        h4_ = h4.group(0)
Daniel Scheffler's avatar
Daniel Scheffler committed
397
398
        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
399
400
        h4max_pix = re.findall("QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall("QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
401
        if not h4max_rad:
Daniel Scheffler's avatar
Daniel Scheffler committed
402
403
404
405
            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
406
407
408
            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
409
410
411
412
413
414
        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])
415
416
417
418
419
420
            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
421
422
        self.wvlUnit = 'Nanometers'
        # Exact values calculated based in SRFs
423
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
Daniel Scheffler's avatar
Daniel Scheffler committed
424
425
426
427
428
429
430
431
432
433
434
435
        # 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
436
                # no irradiation values available
Daniel Scheffler's avatar
Daniel Scheffler committed
437
438
439
440
441
442
443
444
445
446
447
448
449
                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
450
451
452
            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
453
454
                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])]
455
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
456
457
                self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. '
                                  'Found %s' % len(K1_res))
458
459

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

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

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

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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
520
                for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref):
521
                    self.GainsRef[FullLBA.index(ii)] = ind
Daniel Scheffler's avatar
Daniel Scheffler committed
522
                for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref):
523
524
525
526
527
528
529
                    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
530
        if not h10_UL:  # Landsat8
531
532
533
534
            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
535
536
537
538
539
540
541
542
543
        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]]))
544
545

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

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

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

598
599
600
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)
601
602
603
604
605
606
607

        # mtl.close()

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

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

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

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

626
        # Satellite
Daniel Scheffler's avatar
Daniel Scheffler committed
627
628
        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"
629

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

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

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

671
        # Earth sun distance
672
673
674
        self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate)

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

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

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

723
724
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
725
        self.PhysUnit = "DN"
726

727
728
729
730
        # 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
731
        self.Offsets = [0, 0, 0, 0, 0]
732
733

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

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

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

787
788
789
790
        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

791
792
793
794
795
796
797
798
799
    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
800

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

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

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

        # EntityID
Daniel Scheffler's avatar
Daniel Scheffler committed
835
836
837
838
        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)]
839
840

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

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

863
864
        # Units
        self.ScaleFactor = 1
Daniel Scheffler's avatar
Daniel Scheffler committed
865
        self.PhysUnit = "DN"
866
867

        # Gains/Offsets
Daniel Scheffler's avatar
Daniel Scheffler committed
868
869
870
871
        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)
872
873
874
875
        for x, i in enumerate(h41):
            self.Gains.append(float(i))
            self.Offsets.append(float(h42[x]))