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
from typing import List, TYPE_CHECKING  # noqa F401  # flake8 issue

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 sicor.options import get_options as get_ac_options

from ..options.config import GMS_config as CFG
from ..io.input_reader import open_specific_file_within_archive, Solar_Irradiance_reader, SRF_reader
from ..io.output_writer import enviHdr_keyOrder
from ..algorithms import geoprocessing as GEOP
from ..misc import helper_functions as HLP_F
from ..misc import database_tools as DB_T
from ..misc.path_generator import path_generator, get_path_ac_options
from ..misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene, datasetid_to_sat_sen
from ..misc.exceptions import ACNotSupportedError

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

__author__ = 'Daniel Scheffler', 'Robert Behling'


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

        # unpack GMS_identifier
        self.GMS_identifier = GMS_id
        self.image_type = GMS_id.image_type
        self.Satellite = GMS_id.satellite
        self.Sensor = GMS_id.sensor
        self.Subsystem = GMS_id.subsystem
        self.proc_level = GMS_id.proc_level
        self.logger = GMS_id.logger

        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 = {}  # Dict with offset for each band (radiance)
        self.Gains = {}  # Dict with gain for each band (radiance)
        self.OffsetsRef = {}  # Dict with offset for each band for conversion to Reflectance (Landsat8)
        self.GainsRef = {}  # Dict with gain for each band for conversion to Reflectance (Landsat8)
        self.CWL = {}
        self.FWHM = {}
        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}

        self.version_gms_preprocessing = CFG.version
        self.versionalias_gms_preprocessing = CFG.versionalias

    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)."""
        if not self._AcqDateTime and self.AcqDate and self.AcqTime:
            self._AcqDateTime = datetime.datetime.strptime('%s %s%s' % (self.AcqDate, self.AcqTime, '.000000+0000'),
                                                           '%Y-%m-%d %H:%M:%S.%f%z')

        return self._AcqDateTime

    @AcqDateTime.setter
    def AcqDateTime(self, DateTime):
        # type: (datetime.datetime) -> None
        if isinstance(DateTime, str):
            self._AcqDateTime = datetime.datetime.strptime(DateTime, '%Y-%m-%d %H:%M:%S.%f%z')
        elif isinstance(DateTime, datetime.datetime):
            self._AcqDateTime = DateTime

        self.AcqDate = DateTime.strftime('%Y-%m-%d')
        self.AcqTime = DateTime.strftime('%H:%M:%S')

    @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)

    @property
    def LayerBandsAssignment_full(self):
        # type: () -> list
        """Return complete LayerBandsAssignment without excluding thermal or panchromatic bands.

        NOTE: CFG.sort_bands_by_cwl is respected, so returned list may be sorted by central wavelength
        """
        return get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False, return_fullLBA=True)

    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"

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

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

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

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

        ##########################
        # band specific metadata #
        ##########################

        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

        # 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):  # FIXME does not respect sort_bands_by_cwl
            self.additional[-1].append(ind)
        # Offsets
        for ii, (ind, band) in enumerate(zip(h92, LBA_full_sorted)):
            self.Offsets[band] = float(ind)
        # Gains
        for ii, (ind, band) in enumerate(zip(h93, LBA_full_sorted)):
            # gains in dimap file represent reciprocal 1/gain
            self.Gains[band] = 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 = dict(zip(LBA_full_sorted, 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': dict(zip(LBA_full_sorted, [1855., 1615., 1090., 1680.])),
                    'SPOT1b': dict(zip(LBA_full_sorted, [1845., 1575., 1040., 1690.])),
                    'SPOT2a': dict(zip(LBA_full_sorted, [1865., 1620., 1085., 1705.])),
                    'SPOT2b': dict(zip(LBA_full_sorted, [1865., 1615., 1090., 1670.])),
                    'SPOT3a': dict(zip(LBA_full_sorted, [1854., 1580., 1065., 1668.])),
                    'SPOT3b': dict(zip(LBA_full_sorted, [1855., 1597., 1067., 1667.])),
                    'SPOT4a': dict(zip(LBA_full_sorted, [1843., 1568., 1052., 233., 1568.])),
                    'SPOT4b': dict(zip(LBA_full_sorted, [1851., 1586., 1054., 240., 1586.])),
                    'SPOT5a': dict(zip(LBA_full_sorted, [1858., 1573., 1043., 236., 1762.])),
                    'SPOT5b': dict(zip(LBA_full_sorted, [1858., 1575., 1047., 234., 1773.]))}
                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
            sensorcode = get_GMS_sensorcode(self.GMS_identifier)[:2]
            if sensorcode in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b', 'SPOT3a', 'SPOT3b']:
                self.CWL = dict(zip(LBA_full_sorted, [545., 638., 819., 615.]))
            elif sensorcode in ['SPOT4a', 'SPOT4b']:
                self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 645.]))
            elif sensorcode == ['SPOT5a', 'SPOT5b']:
                self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 595.]))

        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"

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

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

        ##########################
        # band specific metadata #
        ##########################
        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

        # Gains and Offsets
        h4 = re.search("GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I)
        h4_ = h4.group(0)
        h4max_rad = re.findall(" LMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4min_rad = re.findall(" LMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)  # space in front IS needed
        h4max_pix = re.findall("QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        h4min_pix = re.findall("QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I)
        if not h4max_rad:
            h4max_rad = re.findall(" RADIANCE_MAXIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4min_rad = re.findall(" RADIANCE_MINIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_,
                                   re.I)  # space in front IS needed
            h4max_pix = re.findall("QUANTIZE_CAL_MAX_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
            h4min_pix = re.findall("QUANTIZE_CAL_MIN_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I)
        # additional: LMAX, LMIN, QCALMAX, QCALMIN
        self.additional.append([["LMAX"], ["LMIN"], ["QCALMAX"], ["QCALMIN"]])
        # Offsets + Gains
        Gains, Offsets = [], []
        for ii, ind in enumerate(h4min_rad):
            Gains.append(
                (float(h4max_rad[ii]) - float(h4min_rad[ii])) / (float(h4max_pix[ii]) - float(h4min_pix[ii])))
            Offsets.append(float(ind) - float(h4min_pix[ii]) * Gains[-1])
            self.additional[-1][-4].append(h4max_rad[ii])
            self.additional[-1][-3].append(h4min_rad[ii])
            self.additional[-1][-2].append(h4max_pix[ii])
            self.additional[-1][-1].append(h4min_pix[ii])
        self.Gains = {bN: Gains[i] for i, bN in enumerate(LBA_full_sorted)}
        self.Offsets = {bN: Offsets[i] for i, bN in enumerate(LBA_full_sorted)}

        # Solar irradiance, central wavelengths , full width half maximum per band
        self.wvlUnit = 'Nanometers'
        # Exact values calculated based in SRFs
        self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band()  # 3x dict
        # Provider values
        if not self.SolIrradiance:  # Preconfigured Irradiation and CWL values
            if re.search("[\D]*5", self.Satellite):
                # Landsat 5; resource for center wavelength (6 = thermal)
                self.SolIrradiance = {'1': 1957., '2': 1826., '3': 1554., '4': 1036., '5': 215., '6': 0.0, '7': 80.67}
                self.CWL = {'1': 485., '2': 560., '3': 660., '4': 830., '5': 1650., '6': 11450., '7': 2215.}
            if re.search("[\D]*7", self.Satellite):
                # Landsat 7; resource for center wavelength:
                # http://opticks.org/confluence/display/opticksDev/Sensor+Wavelength+Definitions
                # 6L(thermal),6H(thermal),B8(pan)
                self.SolIrradiance = {'1': 1969., '2': 1840., '3': 1551., '4': 1044., '5': 225.7, '6L': 0.0, '6H': 0.0,
                                      '7': 82.07, '8': 1368.}
                self.CWL = {'1': 482.5, '2': 565., '3': 660., '4': 825., '5': 1650., '6L': 11450., '6H': 11450.,
                            '7': 2215., '8': 710.}
            if re.search("[\D]*8", self.Satellite):  # Landsat 8
                # no irradiation values available
                self.CWL = {'1': 443., '2': 482.6, '3': 561.3, '4': 654.6, '5': 864.6, '6': 1609.1, '7': 2201.2,
                            '8': 591.7, '9': 1373.5, '10': 10903.6, '11': 12003.}
                #    if None in SolIrradiance:

        # Thermal Constants (K1 = [W*m-2*um-1]; K1 = [K])
        if re.search("[\D]*5", self.Satellite):
            # resource: http://www.yale.edu/ceo/Documentation/Landsat_DN_to_Kelvin.pdf
            ThermalConstK1 = [0., 0., 0., 0., 0., 607.76, 0.]
            ThermalConstK2 = [0., 0., 0., 0., 0., 1260.56, 0.]
        if re.search("[\D]*7", self.Satellite):
            # resource: http://www.yale.edu/ceo/Documentation/Landsat_DN_to_Kelvin.pdf
            ThermalConstK1 = [0., 0., 0., 0., 0., 666.09, 666.09, 0., 0.]
            ThermalConstK2 = [0., 0., 0., 0., 0., 1282.71, 1282.71, 0., 0.]
        if re.search("[\D]*8", self.Satellite):  # Landsat 8
            K1_res = re.findall("K1_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
            K2_res = re.findall("K2_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I)
            if len(K1_res) == 2:
                ThermalConstK1 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K1_res[0]), float(K1_res[1])]
                ThermalConstK2 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K2_res[0]), float(K2_res[1])]
            else:
                self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. '
                                  'Found %s' % len(K1_res))
        self.ThermalConstK1 = {bN: ThermalConstK1[i] for i, bN in enumerate(LBA_full_sorted)}
        self.ThermalConstK2 = {bN: ThermalConstK2[i] for i, bN in enumerate(LBA_full_sorted)}

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

                for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref):
                    GainsRef[LBA_full_sorted.index(ii)] = ind
                for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref):
                    OffsetsRef[LBA_full_sorted.index(ii)] = ind

                self.GainsRef = {bN: GainsRef[i] for i, bN in enumerate(LBA_full_sorted)}
                self.OffsetsRef = {bN: OffsetsRef[i] for i, bN in enumerate(LBA_full_sorted)}

        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"

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

        ##########################
        # band specific metadata #
        ##########################

        LBA_full_sorted = HLP_F.sorted_nicely(self.LayerBandsAssignment_full)

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

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

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

        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"

        # Satellite
        self.Satellite = 'Terra'

        # Sensor
        h10 = re.search("OBJECT[\s]*=[\s]*INSTRUMENTSHORTNAME[\s\S]*END_OBJECT[\s]*=[\s]*INSTRUMENTSHORTNAME", h_, re.I)
        self.Sensor = re.search("VALUE[\s]*=[\s]*[\"]([A-Za-z]*)\"", h10.group(0), re.I).group(1)

        # Subsystem
        h5 = re.search("GROUP[\s]*=[\s]*OBSERVATIONMODE[\s\S]*END_GROUP[\s]*=[\s]*OBSERVATIONMODE", h_, re.I)
        avail_subsystems = re.findall('VALUE[\s]*=[\s]*[(]\"([a-zA-Z0-9]*)\", \"([ONF]*)\"', h5.group(0), re.I)
        if subsystem not in ['VNIR1', 'VNIR2', 'SWIR', 'TIR']:
            raise ValueError('Unexpected subsystem >%s<. Unable to specify the correct ASTER subsystem to be '
                             'processed.' % subsystem)
        else:
            if subsystem == 'VNIR1' and avail_subsystems[0][1] == 'ON':
                self.nBands = 3
            if subsystem == 'VNIR2' and avail_subsystems[1][1] == 'ON':
                self.nBands = 1
            if subsystem == 'SWIR' and avail_subsystems[2][1] == 'ON':
                self.nBands = 6
            if subsystem == 'TIR' and avail_subsystems[3][1] == 'ON':
                self.nBands = 5
        self.Subsystem = subsystem

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

        # Ground resolution
        re_res_GSD = re.findall('OBJECT[\s]*=[\s]*SPATIALRESOLUTION[\s\S]*VALUE[\s]*=[\s]*[\S]{1}([1-9][0-9]), '
                                '([1-9][0-9]), ([1-9][0-9])[\s\S]*END_OBJECT[\s]*=[\s]*SPATIALRESOLUTION', h_, re.I)
        idx_subsystem = \
            0 if subsystem[:4] == 'VNIR' else \
            1 if subsystem[:4] == 'SWIR' else \
            2 if subsystem[:4] == 'TIR' else None
        self.gResolution = float(re_res_GSD[0][idx_subsystem]) \
            if re_res_GSD and idx_subsystem is not None else self.gResolution

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

        # SunAzimuth SunElevation
        h7 = re.search("OBJECT[\s]*=[\s]*SOLARDIRECTION[\s\S]*END_OBJECT[\s]*=[\s]*SOLARDIRECTION", h_, re.I)
        h71 = re.search("VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)", h7.group(0), re.I)
        self.SunAzimuth = float(h71.group(1))
        self.SunElevation = float(h71.group(2))

        # data Quality
        h8 = re.findall(r"GROUP[\s]*=[\s]*DATAQUALITY[1-9][0-4BN]?[^.]*END_GROUP[\s]*=[\s]*DATAQUALITY[1-9][0-4BN]?",
                        h_, re.I)
        h81 = re.findall('VALUE[\s]*=[\s]*([^\n]*)', "\n".join(h8), re.I)
        self.Quality.append(["NoOfBadPixel:"])
        for i in h81:
            self.Quality[-1].append(i)