metadata.py

# -*- coding: utf-8 -*-
"""Module 'metadata' for handling any type of metadata of GeoMultiSens compatible sensor systems."""

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

import collections
import datetime
import glob
import math
import os
import re
import sys
import warnings
import iso8601
import xml.etree.ElementTree as ET

import numpy as np
import pyproj
from matplotlib import dates as mdates
from pyorbital import astronomy

from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.geo.projection import WKT2EPSG

from gms_preprocessing.options.config import GMS_config as CFG
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
from gms_preprocessing.misc import helper_functions as HLP_F
from gms_preprocessing.misc import database_tools as DB_T
from gms_preprocessing.misc.path_generator import path_generator, get_path_ac_options
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


from sicor.options import get_options as get_ac_options

__author__ = 'Daniel Scheffler', 'Robert Behling'


class METADATA(object):
    def __init__(self, GMS_identifier):
        # private attributes
        self._AcqDateTime = None

        # unpack GMS_identifier
        self.GMS_identifier = GMS_identifier
        self.image_type = GMS_identifier['image_type']
        self.Satellite = GMS_identifier['Satellite']
        self.Sensor = GMS_identifier['Sensor']
        self.Subsystem = GMS_identifier['Subsystem']
        self.proc_level = GMS_identifier['proc_level']
        self.logger = GMS_identifier['logger']

        self.Dataname = ''
        self.FolderOrArchive = ''
        self.Metafile = ""  # File containing image metadata (automatically found)
        self.EntityID = ""  # ID to identify the original scene
        self.SceneID = ''  # postgreSQL-database identifier
        self.Sensormode = ""
        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)
        self.CWL = []
        self.FWHM = []
        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
        self.SolIrradiance = []
        self.ThermalConstK1 = []
        self.ThermalConstK2 = []
        self.EarthSunDist = 1.0
        # viewing zenith angle of the Sensor (offNadir) [Deg] (usually) in Case of RapidEye "+"
        # being East and “-” being West
        self.ViewingAngle = -99.0
        self.ViewingAngle_arrProv = {}
        # 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
        self.IncidenceAngle_arrProv = {}
        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 = []
        self.PhysUnit = "DN"
        # 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
        self.CS_TYPE = ""
        self.CS_DATUM = ""
        self.CS_UTM_ZONE = -99.
        self.WRS_path = -99.
        self.WRS_row = -99.
        # List of Corner Coordinates in order of Lon/Lat/DATA_X/Data_Y for all given image corners
        self.CornerTieP_LonLat = []
        self.CornerTieP_UTM = []
        self.LayerBandsAssignment = []  # List of spectral bands in the same order as they are stored in the rasterfile.
        self.additional = []
        self.image_type = 'RSD'
        self.orbitParams = {}
        self.overpassDurationSec = -99.
        self.scene_length = -99.
        self.rows = -99.
        self.cols = -99.
        self.bands = -99.
        self.nBands = -99.
        self.map_info = []
        self.projection = ""
        self.wvlUnit = ""
        self.spec_vals = {'fill': None, 'zero': None, 'saturated': None}

    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)
        elif re.search("Sentinel-2A", self.Satellite, re.I) or re.search("Sentinel-2B", self.Satellite, re.I):
            self.Read_Sentinel2_xmls()
        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:
            raise NotImplementedError("%s is not a supported sensor or sensorname is misspelled." % self.Satellite)

        return self

    @property
    def AcqDateTime(self):
        """Returns a datetime.datetime object containing date, time and timezone (UTC time)."""

        return self._AcqDateTime

    @AcqDateTime.setter
    def AcqDateTime(self, DateTime):
        # type: (datetime.datetime) -> None
        self._AcqDateTime = DateTime
        self.AcqDate = datetime.datetime.strftime(DateTime, format='%Y-%m-%d')
        self.AcqTime = datetime.datetime.strftime(DateTime, format='%H:%M:%S')

    @property
    def overview(self):
        attr2include = \
            ['Dataname', 'Metafile', 'EntityID', 'Satellite', 'Sensor', 'Subsystem', 'gResolution', 'AcqDate',
             'AcqTime',
             'CWL', 'FWHM', 'Offsets', 'Gains', 'ProcLCode', 'ProcLStr', 'SunElevation', 'SunAzimuth',
             'ViewingAngle', 'IncidenceAngle', 'FOV', 'SolIrradiance', 'ThermalConstK1', 'ThermalConstK2',
             'EarthSunDist', 'Quality', 'PhysUnit', 'additional', 'GainsRef', 'OffsetsRef', 'CornerTieP_LonLat',
             'CS_EPSG',
             'CS_TYPE', 'CS_DATUM', 'CS_UTM_ZONE', 'LayerBandsAssignment']
        return collections.OrderedDict((attr, getattr(self, attr)) for attr in attr2include)

    def Read_SPOT_dimap2(self):
        """----METHOD_2------------------------------------------------------------
        # read metadata from spot dimap file
        """

        # self.default_attr()
        if os.path.isdir(self.FolderOrArchive):
            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
            self.Metafile = glob_res[0]
            dim_ = open(self.Metafile, "r").read()
        else:  # archive
            dim_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/scene01/metadata.dim')

        # 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)
        self.additional.append(["SpecialValues:"])
        for ii, ind in enumerate(h1):
            self.additional[0].append(["%s:%s" % (ind, h2[ii])])

        # EntityID
        h3 = re.search("<SOURCE_ID>([a-zA-Z0-9]*)</SOURCE_ID>", dim_, re.I)
        self.EntityID = h3.group(1)

        # AcqDate
        h4 = re.search("<IMAGING_DATE>([0-9-]*)</IMAGING_DATE>", dim_, re.I)
        AcqDate = h4.group(1)

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

        # AcqTime
        h5 = re.search("<IMAGING_TIME>([0-9:]*)</IMAGING_TIME>", dim_, re.I)
        AcqTime = h5.group(1)

        # 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')

        # Satellite + Sensor
        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]*)"
                       "</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)])

        # Angles: incidence angle, sunAzimuth, sunElevation
        h7 = re.search("<INCIDENCE_ANGLE>(.*)</INCIDENCE_ANGLE>[\s\S]*<SUN_AZIMUTH>(.*)</SUN_AZIMUTH>[\s\S]"
                       "*<SUN_ELEVATION>(.*)</SUN_ELEVATION>", dim_, re.I)
        self.IncidenceAngle = float(h7.group(1))
        self.SunAzimuth = float(h7.group(2))
        self.SunElevation = float(h7.group(3))

        # Viewing Angle: not always included in the metadata.dim file
        h8 = re.search("<VIEWING_ANGLE>(.*)</VIEWING_ANGLE>", dim_, re.I)
        if type(h8).__name__ == 'NoneType':
            self.ViewingAngle = 0
        else:
            self.ViewingAngle = float(h8.group(1))

        # Field of View
        self.FOV = get_FieldOfView(self.GMS_identifier)

        # Units
        self.ScaleFactor = 1
        self.PhysUnit = "DN"

        # Gains and Offsets
        h9 = re.search("<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I)
        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:"])
        for ii, ind in enumerate(h91):
            self.additional[-1].append(ind)
        # Offsets
        for ii, ind in enumerate(h92):
            self.Offsets.append(float(ind))
        # 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
        # noinspection PyBroadException
        try:
            self.nBands = (np.mean([len(self.Gains), len(self.Offsets)])
                           if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else
                           GEOP.GEOPROCESSING(self.Dataname, self.logger).bands)
        except Exception:
            self.logger.warning('Unable to get number of bands for the dataset %s Provider values are used for '
                                'solar irradiation, CWL and FWHM!.' % self.Dataname)
        self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands)

        # Exact values calculated based in SRFs
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
        # Provider values
        if not self.SolIrradiance:
            h10 = re.search("<Solar_Irradiance>[\s\S]*</Solar_Irradiance>", dim_, re.I)
            h10_ = h10.group(0)
            h101 = re.findall("<SOLAR_IRRADIANCE_VALUE>([^<]*)</SOLAR_IRRADIANCE_VALUE>", h10_, re.I)
            if h101:
                self.SolIrradiance = h101
                #    self.additional.append(["Solar Irradiance per band:"])
                #    for ii,ind in enumerate(h101):
                #       self.additional[-1].append(ind)
            else:  # Preconfigured Irradiation values
                spot_irr_dic = {
                    '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.]}
                self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)]
            # Preconfigured CWLs --   # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined
            # 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
            if get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b',
                                                               'SPOT3a', 'SPOT3b']:
                self.CWL = [545., 638., 819., 615.]
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] in ['SPOT4a', 'SPOT4b']:
                self.CWL = [540., 650., 835., 1630., 645.]
            elif get_GMS_sensorcode(self.GMS_identifier)[:2] == ['SPOT5a', 'SPOT5b']:
                self.CWL = [540., 650., 835., 1630., 595.]

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

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

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

        # Corner Coordinates
        h121 = re.findall("<TIE_POINT_CRS_X>(.*)</TIE_POINT_CRS_X>", dim_, re.I)
        h122 = re.findall("<TIE_POINT_CRS_Y>(.*)</TIE_POINT_CRS_Y>", dim_, re.I)
        if len(h121) == 4 and self.CS_TYPE == 'LonLat':
            # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR)
            self.CornerTieP_LonLat.append(tuple([float(h121[0]), float(h122[0])]))  # UL
            self.CornerTieP_LonLat.append(tuple([float(h121[1]), float(h122[1])]))  # UR
            self.CornerTieP_LonLat.append(tuple([float(h121[3]), float(h122[3])]))  # LL
            self.CornerTieP_LonLat.append(tuple([float(h121[2]), float(h122[2])]))  # LR
            #    ul_lon,ul_lat = self.inDs.GetGCPs()[0].GCPX,self.inDs.GetGCPs()[0].GCPY # funktioniert nur bei SPOT
            #    lr_lon,lr_lat = self.inDs.GetGCPs()[2].GCPX,self.inDs.GetGCPs()[2].GCPY

        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

    def Read_LANDSAT_mtltxt(self, LayerBandsAssignment):
        """----METHOD_3------------------------------------------------------------
        read metadata from LANDSAT metafile: <dataname>.MTL.txt. Metadatafile of LPGS processing chain
        :param LayerBandsAssignment:
        """

        # self.default_attr()
        self.LayerBandsAssignment = LayerBandsAssignment
        self.nBands = len(LayerBandsAssignment)
        if os.path.isdir(self.FolderOrArchive):
            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
            self.Metafile = glob_res[0]
            mtl_ = open(self.Metafile, "r").read()
        else:  # archive
            mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt')

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

        # Satellite + Sensor + Sensor Mode
        h2 = re.search('SPACECRAFT_ID = "([a-zA-Z0-9_]*)"[\s]*SENSOR_ID = "([a-zA-Z0-9+]*)"[\s]*'
                       'SENSOR_MODE = "([\S]*)"', mtl_, re.I)
        if h2:
            self.Satellite = 'Landsat-%s' % re.search('LANDSAT[\D]*([0-9])', h2.group(1), re.I).group(1)
            self.Sensor = h2.group(2)
            self.Sensormode = h2.group(3)
        else:
            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

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

        # Acquisition Date + Time
        h3 = re.search('ACQUISITION_DATE = ([0-9-]*)[\s]*SCENE_CENTER_SCAN_TIME = "?([0-9:]*)"?', mtl_, re.I)
        if h3 is None:
            h3 = re.search('DATE_ACQUIRED = ([0-9-]*)[\s]*SCENE_CENTER_TIME = "?([0-9:]*)"?', mtl_, re.I)
        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')

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

        # Units
        self.ScaleFactor = 1
        self.PhysUnit = "DN"

        # Gains and Offsets
        h4 = re.search("GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
        h4_ = h4.group(0)
        h4max_rad = re.findall(" LMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4min_rad = re.findall(" LMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4max_pix = re.findall("QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall("QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        if not h4max_rad:
            h4max_rad = re.findall(" RADIANCE_MAXIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4min_rad = re.findall(" RADIANCE_MINIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4max_pix = re.findall("QUANTIZE_CAL_MAX_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
            h4min_pix = re.findall("QUANTIZE_CAL_MIN_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
        # additional: LMAX, LMIN, QCALMAX, QCALMIN
        self.additional.append([["LMAX"], ["LMIN"], ["QCALMAX"], ["QCALMIN"]])
        # Offsets + Gains
        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])
            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
        self.wvlUnit = 'Nanometers'
        # Exact values calculated based in SRFs
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()
        # Provider values
        if not self.SolIrradiance:  # Preconfigured Irradiation 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
                # no irradiation values available
                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
            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:
                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])]
            else:
                self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. '
                                  'Found %s' % len(K1_res))

                # Angles: incidence angle, sunAzimuth, sunElevation, field of view
        h5 = re.search("SUN_AZIMUTH = ([\S]*)[\s]*SUN_ELEVATION = ([\S]*)", mtl_, re.I)
        self.SunAzimuth = float(h5.group(1))
        self.SunElevation = float(h5.group(2))
        self.FOV = get_FieldOfView(self.GMS_identifier)

        # Quality
        h6 = re.search("GROUP = CORRECTIONS_APPLIED[\s\S]*END_GROUP = CORRECTIONS_APPLIED", mtl_, re.I)
        if h6 is None:
            h6 = re.search("GROUP = IMAGE_ATTRIBUTES[\s\S]*END_GROUP = IMAGE_ATTRIBUTES", mtl_, re.I)

        h6_ = h6.group(0)
        h61 = (h6_.split("\n"))
        x = 0
        for i in h61:
            if x == 0 or x + 1 == len(h61):
                pass
            else:
                i_ = i.strip().replace('"', "")
                self.Quality.append(i_.split(" = "))
            x += 1

            # Additional: coordinate system, geom. Resolution
        h7 = re.search("GROUP = PROJECTION_PARAMETERS[\s\S]*END_GROUP = L1_METADATA_FILE", mtl_, re.I)
        h7_ = h7.group(0)
        h71 = (h7_.split("\n"))
        for x, i in enumerate(h71):
            if re.search("Group", i, re.I):
                pass
            else:
                i_ = i.strip().replace('"', "")
                self.additional.append(i_.split(" = "))
        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
        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
        # Landsat8
        h8 = re.search("ROLL_ANGLE = ([\S]*)", mtl_, re.I)
        if h8:
            self.ViewingAngle = float(h8.group(1))

            # reflectance coefficients (Landsat8)
        h9 = re.search("GROUP = RADIOMETRIC_RESCALING[\s\S]*END_GROUP = RADIOMETRIC_RESCALING", mtl_, re.I)
        if h9:
            h9_ = h9.group(0)
            h9gain_ref = re.findall(" REFLECTANCE_MULT_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9gain_ref_bandNr = re.findall(" REFLECTANCE_MULT_BAND_([0-9]+) = [\S]*", h9_, re.I)
            h9offs_ref = re.findall(" REFLECTANCE_ADD_BAND_[0-9]+ = ([\S]*)", h9_, re.I)
            h9offs_ref_bandNr = re.findall(" REFLECTANCE_ADD_BAND_([0-9]+) = [\S]*", h9_, re.I)
            if h9gain_ref:
                self.GainsRef = [None] * len(self.Gains)
                self.OffsetsRef = [None] * len(self.Offsets)
                FullLBA = get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False)

                for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref):
                    self.GainsRef[FullLBA.index(ii)] = ind
                for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref):
                    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)
        if not h10_UL:  # Landsat8
            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)
        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]]))

        # WRS path/row
        h11_p = re.search('WRS_PATH = ([0-9]*)', mtl_, re.I)
        if h11_p:
            self.WRS_path = h11_p.group(1)
        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)
            self.WRS_row = int(h11_r.group(1))
        else:
            self.WRS_row = int(h11_r1.group(1)) if h11_r1.group(1) == h11_r2.group(1) else self.WRS_row

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

            if len(result) == 1:  # e.g. [('LE71950282003121EDC00',)]
                self.EntityID = result[0][0]
            elif len(result) == 0:
                self.logger.warning("Scene ID could not be assigned because no dataset matching to the query "
                                    "parameters could be found in metadata database.")
            else:  # e.g. [('LE71950282003121EDC00',), ('LE71950282003105ASN00',)]
                self.logger.warning("Scene ID could not be assigned because %s datasets matching to the query "
                                    "parameters were found in metadata database." % len(result))
        # if self.EntityID=='':
        #     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 ''
        #     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']

        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

        # mtl.close()

    def Read_RE_metaxml(self):
        """----METHOD_4------------------------------------------------------------
        read metadata from RapidEye metafile: <dataname>metadata.xml
        """

        # self.default_attr()
        if os.path.isdir(self.FolderOrArchive):
            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
            self.Metafile = glob_res[0]
            mxml_ = open(self.Metafile, "r").read()
        else:  # archive
            mxml_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*_metadata.xml')

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

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

        # Satellite
        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"

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

        # Acquisition Parameters: Incidence Angle, SunAzimuth, SunElevation, ViewingAngle, FieldOfView, AcqDate, AcqTime
        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'
                           'Time>([0-9-]*)T([0-9:]*)[\S]*</re:acquisitionDateTime>', mxml_, re.I)
            self.IncidenceAngle = float(h5.group(1))
            self.SunAzimuth = float(h5.group(2))
            self.SunElevation = float(h5.group(3))
            # 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))])
            self.ViewingAngle = float(h5.group(5))
            self.FOV = get_FieldOfView(self.GMS_identifier)
            AcqDate = h5.group(6)
            AcqTime = h5.group(7)

        except AttributeError:
            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)  #
            self.IncidenceAngle = float(h5.group(1))
            self.SunAzimuth = float(h5.group(2))
            self.SunElevation = float(h5.group(3))
            self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))])
            self.ViewingAngle = 9999.99
            self.FOV = get_FieldOfView(self.GMS_identifier)
            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')

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

        # Additional: Projection
        h6 = re.search("<re:geodeticDatum>([\S]*)</re:geodeticDatum>[\s]*<re:projection>([\S^<\s]*)</re:projection>",
                       mxml_, re.I)
        try:
            self.additional.append(["SpatialReference", ["geodeticDatum", h6.group(1)], ["projection", h6.group(2)]])
        except AttributeError:  # Level1-B data has no geodaticDatum
            pass

        # Corrections applied + Quality
        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)
        # fuer IRIS ihre Daten
        if h7 is None:
            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)],
                                    ["atmosphericCorrectionApplied", h7.group(4)]])
        else:
            self.additional.append(["Corrections", ["radiometricCorrectionApplied", h7.group(1)],
                                    ["radiometricCalibrationVersion", h7.group(2)], ["geoCorrectionLevel", h7.group(3)],
                                    ["elevationCorrectionApplied", h7.group(4)],
                                    ["atmosphericCorrectionApplied", h7.group(5)]])
            self.Quality.append(["geomProductAccuracy[m]:", str(
                round(float(h7.group(6)), 1))])  # Estimated product horizontal CE90 uncertainty [m]

        # additional. missing lines, suspectlines, binning, shifting, masking
        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)

        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]

        # Units
        self.ScaleFactor = 1
        self.PhysUnit = "DN"

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

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

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

            # Coordinate Reference System
        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
        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)
        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])]))

        self.orbitParams = get_orbit_params(self.GMS_identifier)
        self.filter_layerdependent_metadata()
        self.spec_vals = get_special_values(self.GMS_identifier)

    def Read_ASTER_hdffile(self, subsystem):
        """#----METHOD_5------------------------------------------------------------
        read metadata from ASTER hdf
        input:
            hdffile:
            subsystem:
        output:
        :param subsystem:
        """

        #    self.default_attr()
        dat_ = open(self.FolderOrArchive, "r").read() if sys.version_info[0] < 3 else \
            open(self.FolderOrArchive, "rb").read().decode('latin-1')

        # Split meta from raster data
        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_])

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

        # EntityID
        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)]

        # Viewing Angle
        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:
            self.ViewingAngle = float(h21[0])
        else:
            self.ViewingAngle = -99.0

        # additional GainMode
        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:"]])
        for x, i in enumerate(h31[:15]):
            self.additional[-1][-2].append(i)
            self.additional[-1][-1].append(gains[i][x])

        # Units
        self.ScaleFactor = 1
        self.PhysUnit = "DN"

        # Gains/Offsets
        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)
        for x, i in