Commit be82838c authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed issue #66 (Number of wavelengths does not match number bands in L2C header file).

-> required to revise metadata storage for layer dependent metadata (they are now stored as dicts instead of lists)
parent 3e0b05f0
......@@ -314,14 +314,9 @@ class L1A_object(GMS_object):
proc_opt_therm = sorted(list(set(self.dict_LayerOptTherm.values())))
assert proc_opt_therm in [['optical', 'thermal'], ['optical'], ['thermal']]
# get input parameters
off, gai, irr, zen, esd, inFill, inZero, inSaturated, k1, k2 = \
[self.MetaObj.Offsets, self.MetaObj.Gains, self.MetaObj.SolIrradiance,
90 - float(self.MetaObj.SunElevation),
self.MetaObj.EarthSunDist, self.MetaObj.spec_vals['fill'], self.MetaObj.spec_vals['zero'],
self.MetaObj.spec_vals['saturated'], self.MetaObj.ThermalConstK1, self.MetaObj.ThermalConstK2]
(rS, rE), (cS, cE) = self.arr_pos if self.arr_pos else ((0, self.shape_fullArr[0]), (0, self.shape_fullArr[1]))
inFill, inZero, inSaturated = \
self.MetaObj.spec_vals['fill'], self.MetaObj.spec_vals['zero'], self.MetaObj.spec_vals['saturated']
(rS, rE), (cS, cE) = self.arr_pos or ((0, self.shape_fullArr[0]), (0, self.shape_fullArr[1]))
# in_mem = hasattr(self, 'arr') and isinstance(self.arr, np.ndarray)
# if in_mem:
# (rS, rE), (cS, cE) =
......@@ -336,7 +331,10 @@ class L1A_object(GMS_object):
"DN2RadRef-Input data have %s bands although %s bands are specified in self.LayerBandsAssignment." \
% (self.arr.bands, len(self.LayerBandsAssignment))
# perform conversion if needed
################################
# perform conversion if needed #
################################
data_optical, data_thermal, optical_bandsList, thermal_bandsList = None, None, [], []
for optical_thermal in ['optical', 'thermal']:
if optical_thermal not in self.dict_LayerOptTherm.values():
......@@ -359,7 +357,21 @@ class L1A_object(GMS_object):
if arr_desc == 'DN':
self.log_for_fullArr_or_firstTile('Calculating %s...' % conv, subset)
# get input parameters #
########################
off = [self.MetaObj.Offsets[b] for b in self.LayerBandsAssignment]
gai = [self.MetaObj.Gains[b] for b in self.LayerBandsAssignment]
irr = [self.MetaObj.SolIrradiance[b] for b in self.LayerBandsAssignment]
zen, esd = 90 - float(self.MetaObj.SunElevation), self.MetaObj.EarthSunDist
k1 = [self.MetaObj.ThermalConstK1[b] for b in self.LayerBandsAssignment]
k2 = [self.MetaObj.ThermalConstK2[b] for b in self.LayerBandsAssignment]
OFF, GAI, IRR, K1, K2 = [list(np.array(par)[bList]) for par in [off, gai, irr, k1, k2]]
# perform conversion #
######################
res = \
GEOP.DN2Rad(inArray, OFF, GAI, inFill, inZero, inSaturated) if conv == 'Rad' else \
GEOP.DN2TOARef(inArray, OFF, GAI, IRR, zen, esd, inFill, inZero,
......@@ -410,7 +422,10 @@ class L1A_object(GMS_object):
else:
data_thermal, thermal_bandsList = res, bList
# combine results from optical and thermal data
#################################################
# combine results from optical and thermal data #
#################################################
if data_optical is not None and data_thermal is not None:
bands_opt, bands_therm = [1 if len(d.shape) == 2 else d.shape[2] for d in [data_optical, data_thermal]]
r, c, b = data_optical.shape[0], data_optical.shape[1], bands_opt + bands_therm
......
......@@ -568,9 +568,11 @@ class AtmCorr(object):
solar_irradiance = {}
for inObj in self.inObjs:
for bandN, bandIdx in inObj.arr.bandnames.items():
for bandN in inObj.arr.bandnames:
lba_key = bandN[2:] if bandN.startswith('B0') else bandN[1:]
if bandN not in solar_irradiance:
solar_irradiance[bandN] = inObj.meta_odict['SolIrradiance'][bandIdx]
solar_irradiance[bandN] = inObj.MetaObj.SolIrradiance[lba_key]
return solar_irradiance
def _meta_get_viewing_zenith(self):
......@@ -938,21 +940,28 @@ class AtmCorr(object):
ac_bandNs = [bandN for bandN in inObj.arr.bandnames if bandN in self.results.data_ac.keys()]
out_LBA = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs]
# update metadata
# update metadata #
###################
inObj.arr_desc = 'BOA_Ref'
inObj.MetaObj.bands = len(self.results.data_ac)
inObj.MetaObj.PhysUnit = 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef
inObj.MetaObj.LayerBandsAssignment = out_LBA
inObj.LayerBandsAssignment = out_LBA
inObj.MetaObj.filter_layerdependent_metadata()
inObj.meta_odict = inObj.MetaObj.to_odict() # actually auto-updated by getter
# join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config
##################################################################################
# join SURFACE REFLECTANCE as 3D int16 array, scaled to scale factor from config #
##################################################################################
# FIXME AC output nodata values = 0 -> new nodata areas but mask not updated
oF_refl, oZ_refl, oS_refl = get_outFillZeroSaturated(inObj.arr.dtype)
surf_refl = np.dstack((self.results.data_ac[bandN] for bandN in ac_bandNs))
surf_refl *= CFG.scale_factor_BOARef # scale using scale factor (output is float16)
# FIXME really set AC nodata values to GMS outZero?
surf_refl[nodata] = oZ_refl # overwrite AC nodata values with GMS outZero
# apply the original nodata mask (indicating background values)
surf_refl[np.array(inObj.mask_nodata).astype(np.int8) == 0] = oF_refl
......@@ -962,8 +971,7 @@ class AtmCorr(object):
surf_refl[
surf_refl == self.results.bad_data_value] = oF_refl # FIXME meaningful to set AC nans to -9999?
# overwrite LayerBandsAssignment and use inObj.arr setter to generate a GeoArray
inObj.LayerBandsAssignment = out_LBA
# use inObj.arr setter to generate a GeoArray
inObj.arr = surf_refl.astype(inObj.arr.dtype) # -> int16 (also converts NaNs to 0 if needed
else:
......
......@@ -62,7 +62,6 @@ class L2B_object(L2A_object):
src_dsID = sat_sen_to_datasetid(self.satellite, self.sensor)
src_cwls = self.meta_odict['wavelength']
# FIXME exclude or include thermal bands; respect sorted CWLs in context of LayerBandsAssignment
tgt_cwls = CFG.target_CWL
tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref)
# NOTE: get target LBA at L2A, because spectral characteristics of target sensor do not change after AC
tgt_LBA = get_LayerBandsAssignment(
......@@ -78,7 +77,7 @@ class L2B_object(L2A_object):
"the spectral refernce sensor.")
return
if src_cwls == tgt_cwls or (self.satellite == tgt_sat and self.sensor == tgt_sen):
if src_cwls == CFG.target_CWL or (self.satellite == tgt_sat and self.sensor == tgt_sen):
# FIXME catch the case if LayerBandsAssignments are unequal with np.take
self.logger.info("Spectral homogenization has been skipped because the current spectral characteristics "
"are already equal to the target sensor's.")
......@@ -93,7 +92,7 @@ class L2B_object(L2A_object):
if method == 'LI' or CFG.datasetid_spectral_ref is None:
# linear interpolation (if intended by user or in case of custom spectral characteristics of target sensor)
# -> no classifier for that case available -> linear interpolation
im = SpH.interpolate_cube(self.arr, src_cwls, tgt_cwls, kind='linear')
im = SpH.interpolate_cube(self.arr, src_cwls, CFG.target_CWL, kind='linear')
if CFG.spechomo_estimate_accuracy:
self.logger.warning("Unable to compute any error information in case spectral homogenization algorithm "
......@@ -115,7 +114,7 @@ class L2B_object(L2A_object):
compute_errors=CFG.spechomo_estimate_accuracy,
bandwise_errors=CFG.spechomo_bandwise_accuracy,
fallback_argskwargs=dict(
args=dict(source_CWLs=src_cwls, target_CWLs=tgt_cwls,),
args=dict(source_CWLs=src_cwls, target_CWLs=CFG.target_CWL,),
kwargs=dict(kind='linear')
))
......@@ -124,8 +123,9 @@ class L2B_object(L2A_object):
###################
self.LayerBandsAssignment = tgt_LBA
self.meta_odict['wavelength'] = list(tgt_cwls)
self.meta_odict['bands'] = len(tgt_cwls)
self.meta_odict['wavelength'] = list(CFG.target_CWL)
self.meta_odict['bandwidths'] = list(CFG.target_FWHM)
self.meta_odict['bands'] = len(CFG.target_CWL)
if 'band names' in self.meta_odict: # FIXME bug workaround
del self.meta_odict['band names'] # TODO
......@@ -139,8 +139,8 @@ class L2B_object(L2A_object):
if self.ac_errors and self.ac_errors.ndim == 3:
self.logger.info("Performing linear interpolation for 'AC errors' array to match target sensor bands "
"number..")
outarr = \
interp1d(np.array(src_cwls), self.ac_errors, axis=2, kind='linear', fill_value='extrapolate')(tgt_cwls)
outarr = interp1d(np.array(src_cwls), self.ac_errors,
axis=2, kind='linear', fill_value='extrapolate')(CFG.target_CWL)
self.ac_errors = outarr.astype(np.int16)
......
......@@ -193,6 +193,7 @@ def SRF_reader(GMS_identifier, no_thermal=None, no_pan=None, v=False):
LayerBandsAssignment = META.get_LayerBandsAssignment(GMS_identifier, no_thermal=no_thermal, no_pan=no_pan)
SRF_dict = collections.OrderedDict()
SRF_dir = PG.get_path_srf_file(GMS_identifier)
if os.path.isdir(SRF_dir):
for bandname in LayerBandsAssignment:
SRF_path = PG.get_path_srf_file(GMS_identifier, bandname=bandname)
......
......@@ -38,7 +38,7 @@ from sicor.options import get_options as get_ac_options
from ..misc.logging import GMS_logger as DatasetLogger
from ..model.mgrs_tile import MGRS_tile
from ..model.metadata import METADATA, get_dict_LayerOptTherm, metaDict_to_metaODict
from ..model.metadata import METADATA, get_dict_LayerOptTherm, metaDict_to_metaODict, layerdependent_metadata
from ..model.dataset import Dataset
from ..misc import path_generator as PG
from ..misc import database_tools as DB_T
......@@ -739,9 +739,7 @@ class GMS_object(Dataset):
del self.pathGen # must be refreshed because subsystem is now ''
self.close_GMS_loggers() # must also be refreshed because it depends on pathGen
for attrN in ['SolIrradiance', 'CWL', 'FWHM', 'Offsets', 'OffsetsRef', 'Gains', 'GainsRef',
'ThermalConstK1', 'ThermalConstK2', 'ViewingAngle_arrProv', 'IncidenceAngle_arrProv']:
for attrN in layerdependent_metadata:
# combine values from separate subsystems to a single value
attrDic_fullLBA = {}
for GMS_obj in list_GMS_objs:
......
......@@ -63,20 +63,20 @@ class METADATA(object):
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.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.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
......@@ -170,6 +170,15 @@ class METADATA(object):
'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
......@@ -238,6 +247,42 @@ class METADATA(object):
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)
......@@ -245,15 +290,15 @@ class METADATA(object):
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):
for ii, ind in enumerate(h91): # FIXME does not respect sort_bands_by_cwl
self.additional[-1].append(ind)
# Offsets
for ii, ind in enumerate(h92):
self.Offsets.append(float(ind))
for ii, (ind, band) in enumerate(zip(h92, LBA_full_sorted)):
self.Offsets[band] = float(ind)
# Gains
for ii, ind in enumerate(h93):
for ii, (ind, band) in enumerate(zip(h93, LBA_full_sorted)):
# gains in dimap file represent reciprocal 1/gain
self.Gains.append(1 / float(ind))
self.Gains[band] = 1 / float(ind)
# Solar irradiance, central wavelengths , full width half maximum per band
self.wvlUnit = 'Nanometers'
......@@ -276,58 +321,33 @@ class METADATA(object):
h10_ = h10.group(0)
h101 = re.findall("<SOLAR_IRRADIANCE_VALUE>([^<]*)</SOLAR_IRRADIANCE_VALUE>", h10_, re.I)
if h101:
self.SolIrradiance = 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': [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.]}
'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
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
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()
......@@ -350,7 +370,7 @@ class METADATA(object):
else: # archive
mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt')
# Processing Level
# 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)
......@@ -395,72 +415,7 @@ class METADATA(object):
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
# 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))
......@@ -482,7 +437,7 @@ class METADATA(object):
self.Quality.append(i_.split(" = "))
x += 1
# Additional: coordinate system, geom. Resolution
# 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"))
......@@ -508,24 +463,6 @@ class METADATA(object):
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)
......@@ -562,7 +499,8 @@ class METADATA(object):
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})
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]
......@@ -599,6 +537,105 @@ class METADATA(object):
# 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,