Commit 668469f6 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'enhancement/improve_L2A_XML' into 'master'

Enhancement/improve l2 a xml

Closes #53

See merge request !44
parents 34cf6481 9b69ffac
Pipeline #8746 passed with stages
in 13 minutes and 37 seconds
......@@ -2,6 +2,12 @@
History
=======
0.12.x (2020-05-??)
-------------------
* L2A METADATA.XML now contains correct band characteristics, band statistics and information about the merged
VNIR/SWIR cube (fixes issue #53).
0.12.7 (2020-05-12)
-------------------
......
......@@ -43,6 +43,7 @@ from geoarray import GeoArray
from .metadata_sensorgeo import EnMAP_Metadata_L1B_SensorGeo
from ...options.config import EnPTConfig, enmap_xres
from ..srf import SRF
from ...version import __version__
__author__ = ['Daniel Scheffler', 'Stéphane Guillaso', 'André Hollstein']
......@@ -71,6 +72,11 @@ class EnMAP_Metadata_L2A_MapGeo(object):
self.band_stds: Optional[np.ndarray] = None # band-wise standard deviations in unscaled values
self.fileinfos: list = [] # file informations for each file beloning to the EnMAP L2A product
# input validation
if not len(wvls_l2a) == dims_mapgeo[2]:
raise ValueError("The list of center wavelength must be qual to the number of bands in the map geometry "
"dimensions.")
self.proc_level = 'L2A'
self.observation_datetime: datetime = meta_l1b.observation_datetime # Date and Time of observation
# FIXME VZA may be negative in DLR data
......@@ -101,11 +107,11 @@ class EnMAP_Metadata_L2A_MapGeo(object):
# fuse band-wise metadata (sort all band-wise metadata by wavelengths but band number keeps as it is)
# get band index order
wvls_sorted = np.array(sorted(np.hstack([self._meta_l1b.vnir.wvl_center,
self._meta_l1b.swir.wvl_center])))
bandidx_order = np.array([np.argmin(np.abs(wvls_sorted - cwl)) for cwl in wvls_l2a])
wvls_vswir = np.hstack([self._meta_l1b.vnir.wvl_center,
self._meta_l1b.swir.wvl_center])
bandidx_order = np.array([np.argmin(np.abs(wvls_vswir - cwl)) for cwl in wvls_l2a])
self.wvl_center = np.hstack([meta_l1b.vnir.wvl_center, meta_l1b.swir.wvl_center])[bandidx_order]
self.wvl_center = np.array(wvls_l2a)
self.fwhm = np.hstack([meta_l1b.vnir.fwhm, meta_l1b.swir.fwhm])[bandidx_order]
self.gains = np.full((dims_mapgeo[2],), 100) # implies reflectance scaled between 0 and 10000
self.offsets = np.zeros((dims_mapgeo[2],))
......@@ -131,6 +137,8 @@ class EnMAP_Metadata_L2A_MapGeo(object):
self.nrows = dims_mapgeo[0]
self.ncols = dims_mapgeo[1]
self.nwvl = dims_mapgeo[2]
assert self.nwvl == len(self.wvl_center)
common_UL_UR_LL_LR = self.get_common_UL_UR_LL_LR()
self.lon_UL_UR_LL_LR = [lon for lon, lat in common_UL_UR_LL_LR]
self.lat_UL_UR_LL_LR = [lat for lon, lat in common_UL_UR_LL_LR]
......@@ -170,9 +178,9 @@ class EnMAP_Metadata_L2A_MapGeo(object):
def add_band_statistics(self, datastack_vnir_swir: Union[np.ndarray, GeoArray]):
R, C, B = datastack_vnir_swir.shape
# NOTE: DEVIDE by gains to reflectance in percent
self.band_means = np.mean(datastack_vnir_swir.reshape(1, R * C, B), axis=1) / self.gains
self.band_stds = np.mean(datastack_vnir_swir.reshape(1, R * C, B), axis=1) / self.gains
# NOTE: DEVIDE by gains to get reflectance in percent
self.band_means = np.mean(datastack_vnir_swir.reshape((R * C, B)), axis=0) / self.gains
self.band_stds = np.std(datastack_vnir_swir.reshape((R * C, B)), axis=0) / self.gains
def add_product_fileinformation(self, filepaths: List[str], sizes: List[int] = None, versions: List[str] = None):
self.fileinfos = []
......@@ -223,32 +231,40 @@ class EnMAP_Metadata_L2A_MapGeo(object):
self.logger.warning('Currently, the L2A metadata XML file does not contain all relevant keys and contains '
'not updated values!') # FIXME
def uk(path, value):
xml.find(path).text = str(value)
def create_SubElement(_parent, _tag, attrib={}, _text=None, nsmap=None, **_extra):
result = ElementTree.SubElement(_parent, _tag, attrib, nsmap, **_extra)
result.text = _text
return result
############
# metadata #
############
xml.find("metadata/schema/processingLevel").text = self.proc_level
xml.find("metadata/name").text = self.filename_metaxml
# xml.find("metadata/comment").text = 'EnMAP Level 0 Product of datatake 987' # FIXME hardcoded
uk("metadata/schema/processingLevel", self.proc_level)
uk("metadata/name", self.filename_metaxml)
# uk("metadata/comment").text = 'EnMAP Level 0 Product of datatake 987' # FIXME hardcoded
##############
# processing #
##############
xml.find("processing/terrainCorrection").text = 'Yes' # FIXME hardcoded {Yes, No}
xml.find("processing/ozoneValue").text = 'NA' # FIXME {[200-500], NA}
xml.find("processing/season").text = 'NA' # FIXME {summer, winter, NA}
xml.find("processing/productFormat").text = 'GeoTIFF+Metadata' # FIXME hardcoded
uk("processing/terrainCorrection", 'Yes') # FIXME hardcoded {Yes, No}
uk("processing/ozoneValue", 'NA') # FIXME {[200-500], NA})
uk("processing/season", 'NA') # FIXME {summer, winter, NA}
uk("processing/productFormat", 'GeoTIFF+Metadata') # FIXME hardcoded
# {BSQ+Metadata, BIL+Metadata, BIP+Metadata, JPEG2000+Metadata, GeoTiff+Metadata}
xml.find("processing/mapProjection").text = 'UTM_Zone_of_Scene_Center' # FIXME hardcoded
uk("processing/mapProjection", 'UTM_Zone_of_Scene_Center') # FIXME hardcoded
# {UTM_Zone_of_Scene_Center, UTM_Zone_of_Scene_Center(-1), UTM_Zone_of_Scene_Center(+1),
# UTM_Zone_of_Datatake_Center, Geographic, European_Projection_LAEA, NA}
xml.find("processing/DEMDBVersion").text = 'SRTM-C_v4' # FIXME hardcoded
uk("processing/DEMDBVersion", 'SRTM-C_v4') # FIXME hardcoded
# {SRTM-C-X_vv.rr, best-of-DEM_vv.rr, DEM-derivedfrom-Tandem-X_vv.rr, ASTER-GDEM_vv.rr, NA}
xml.find("processing/correctionType").text = 'NA' # FIXME hardcoded {Combined, Land_Mode, Water_Mode, NA}
xml.find("processing/cirrusHazeRemoval").text = 'NA' # FIXME hardcoded {Yes, No}
xml.find("processing/bandInterpolation").text = 'NA' # FIXME hardcoded {Yes, No}
xml.find("processing/waterType").text = 'NA' # FIXME hardcoded {Clear, Turbid, Highly_Turbid, NA}
uk("processing/correctionType", 'NA') # FIXME hardcoded {Combined, Land_Mode, Water_Mode, NA}
uk("processing/cirrusHazeRemoval", 'NA') # FIXME hardcoded {Yes, No}
uk("processing/bandInterpolation", 'NA') # FIXME hardcoded {Yes, No}
uk("processing/waterType", 'NA') # FIXME hardcoded {Clear, Turbid, Highly_Turbid, NA}
########
# base #
......@@ -256,20 +272,27 @@ class EnMAP_Metadata_L2A_MapGeo(object):
# TODO update corner coordinates? DLR just uses the same coordinates like in L1B
# xml.find("base/spatialCoverage" % lbl).text =
xml.find("base/format").text = 'ENMAP_%s' % self.proc_level
xml.find("base/level").text = self.proc_level
xml.find("base/size").text = 'NA' # FIXME Size of product. Attribute unit {byte, Kbyte, Mbyte, Gbyte}
uk("base/format", 'ENMAP_%s' % self.proc_level)
uk("base/level", self.proc_level)
uk("base/size", 'NA') # FIXME Size of product. Attribute unit {byte, Kbyte, Mbyte, Gbyte}
############
# specific #
############
xml.find("specific/code").text = self.proc_level
bi = "specific/bandCharacterisation/bandID/"
for ele, gain in zip(xml.findall(bi + "GainOfBand"), self.gains):
ele.text = str(gain)
for ele, offset in zip(xml.findall(bi + "OffsetOfBand"), self.offsets):
ele.text = str(offset)
uk("specific/code", self.proc_level)
# delete old band characterisation (as it contains a different set of bands)
bChar_root = xml.find('specific/bandCharacterisation')
bChar_root.clear()
# recreate sub-elements for bandCharacterisation with respect to the current set of bands
for i in range(self.nwvl):
sub = ElementTree.SubElement(bChar_root, 'bandID', number=str(i + 1))
for k, val in zip(['wavelengthCenterOfBand', 'FWHMOfBand', 'GainOfBand', 'OffsetOfBand'],
[self.wvl_center[i], self.fwhm[i], self.gains[i], self.offsets[i]]):
ele = ElementTree.SubElement(sub, k)
ele.text = str(val)
###########
# product #
......@@ -279,39 +302,87 @@ class EnMAP_Metadata_L2A_MapGeo(object):
raise ValueError('Product file informations must be added before writing metadata. '
'Call add_product_fileinformation() before!')
# image #
#########
# recreate sub-elements for image (L1B XML contains sub-elements for VNIR and SWIR here, L2A is merged)
im_root = xml.find('product/image')
im_root.clear()
merge = ElementTree.SubElement(im_root, 'merge')
for subEleArgs in [
# tag, attribute dictionary, text
('channels', {}, str(self.nwvl)),
('name', {}, self.filename_data),
('size', {'unit': "Kbyte"}, str([i for i in self.fileinfos if i['name'] == self.filename_data][0]['size'])),
('version', {}, __version__), # we use the EnPT version here
('format', {}, 'binary'),
('dimension',),
('dimensionGeographic',),
('qlChannelsSWIR',),
('qlChannelsVNIR',),
]:
create_SubElement(merge, *subEleArgs)
dim_root = xml.find('product/image/merge/dimension')
for subEleArgs in [
# tag, attribute dictionary, text
('columns', {}, str(self.ncols)),
('rows', {}, str(self.nrows)),
]:
create_SubElement(dim_root, *subEleArgs)
dimgeo_root = xml.find('product/image/merge/dimensionGeographic')
for subEleArgs in [
# tag, attribute dictionary, text
('longitude', {'unit': 'DEG'}, 'NA'), # FIXME
('latitude', {'unit': 'DEG'}, 'NA'), # FIXME 0.3314405
]:
create_SubElement(dimgeo_root, *subEleArgs)
qlswir_root = xml.find('product/image/merge/qlChannelsSWIR')
for subEleArgs in [
# tag, attribute dictionary, text
('red', {}, str(self.preview_wvls_vnir[0])),
('green', {}, str(self.preview_wvls_vnir[1])),
('blue', {}, str(self.preview_wvls_vnir[2])),
]:
create_SubElement(qlswir_root, *subEleArgs)
qlvnir_root = xml.find('product/image/merge/qlChannelsVNIR')
for subEleArgs in [
# tag, attribute dictionary, text
('red', {}, str(self.preview_wvls_vnir[0])),
('green', {}, str(self.preview_wvls_vnir[1])),
('blue', {}, str(self.preview_wvls_vnir[2])),
]:
create_SubElement(qlvnir_root, *subEleArgs)
# quicklook #
#############
from . import L2A_product_props_DLR
for detName, detMetaL1B in zip(['VNIR', 'SWIR'], [self._meta_l1b.vnir, self._meta_l1b.swir]):
lbl = L2A_product_props_DLR['xml_detector_label'][detName]
# FIXME DLR uses L0 filenames for VNIR/SWIR separately?!
xml.find("product/image/%s/name" % lbl).text = detMetaL1B.filename_data
# FIXME this is the size of the VNIR/SWIR stack
size = [F['size'] for F in self.fileinfos if os.path.splitext(F['name'])[0].endswith('-SPECTRAL_IMAGE')][0]
xml.find("product/image/%s/size" % lbl).text = str(size)
# FIXME DLR data dimensions equal neither L2A data nor L1B data
xml.find("product/image/%s/channels" % lbl).text = str(detMetaL1B.nwvl)
xml.find("product/image/%s/dimension/rows" % lbl).text = str(self.nrows)
xml.find("product/image/%s/dimension/columns" % lbl).text = str(self.ncols)
# xml.find("product/image/%s/dimension/dimensionGeographic/longitude" % lbl).text = 'NA' # TODO
# xml.find("product/image/%s/dimension/dimensionGeographic/latitude" % lbl).text = 'NA'
fN_quicklook = self.filename_quicklook_vnir if detName == 'VNIR' else self.filename_quicklook_swir
size_quicklook = [F['size'] for F in self.fileinfos
if os.path.splitext(F['name'])[0].endswith('-QL_%s' % detName)][0]
xml.find("product/quicklook/%s/name" % lbl).text = fN_quicklook
xml.find("product/quicklook/%s/size" % lbl).text = str(size_quicklook)
xml.find("product/quicklook/%s/dimension/rows" % lbl).text = str(self.nrows)
xml.find("product/quicklook/%s/dimension/columns" % lbl).text = str(self.ncols)
# xml.find("product/quicklook/%s/dimension/dimensionGeographic/longitude" % lbl).text = 'NA'
# xml.find("product/quicklook/%s/dimension/dimensionGeographic/latitude" % lbl).text = 'NA'
uk("product/quicklook/%s/name" % lbl, fN_quicklook)
uk("product/quicklook/%s/size" % lbl, str(size_quicklook))
uk("product/quicklook/%s/dimension/rows" % lbl, str(self.nrows))
uk("product/quicklook/%s/dimension/columns" % lbl, str(self.ncols))
# uk("product/quicklook/%s/dimension/dimensionGeographic/longitude" % lbl, 'NA') # FIXME
# uk("product/quicklook/%s/dimension/dimensionGeographic/latitude" % lbl, 'NA') # FIXME
# productFileInformation
########################
# productFileInformation #
##########################
# get L1B product file information
l1b_fileinfos = xmlSubtree2dict(xml, 'product/productFileInformation/')
# clear old L1B file information in XML
pFI_root = xml.findall('product/productFileInformation')[0]
pFI_root = xml.find('product/productFileInformation')
pFI_root.clear()
# recreate sub-elements for productFileInformation according to L2A file information
......@@ -333,10 +404,13 @@ class EnMAP_Metadata_L2A_MapGeo(object):
ele = ElementTree.SubElement(sub, k, **kw)
ele.text = str(fileInfo[k])
# ortho #
#########
# TODO update product/ortho/projection
# {UTM_ZoneX_North, UTM_ZoneX_South (where X in {1..60}), Geographic, LAEA-ETRS89, NA}
xml.find('product/ortho/resolution').text = str(enmap_xres)
xml.find('product/ortho/resampling').text = self.cfg.ortho_resampAlg
uk('product/ortho/resolution', enmap_xres)
uk('product/ortho/resampling', self.cfg.ortho_resampAlg)
# band statistics
#################
......@@ -345,11 +419,21 @@ class EnMAP_Metadata_L2A_MapGeo(object):
raise ValueError('Band statistics have not yet been computed. Compute them first by calling '
'add_band_statistics()!')
bs = "specific/bandStatistics/bandID/"
for ele, mean in zip(xml.findall(bs + "meanReflectance"), self.band_means):
ele.text = str(mean)
for ele, std in zip(xml.findall(bs + "stdDeviation"), self.band_stds):
ele.text = str(std)
# delete old L1B band statistics
bStats_root = xml.find('product/bandStatistics')
bStats_root.clear()
# recreate sub-elements for bandStatistics with respect to the current set of bands
for i in range(self.nwvl):
sub = ElementTree.SubElement(bStats_root, 'bandID', number=str(i + 1))
for k, val in zip(['wavelength', 'mean', 'stdDeviation'],
[self.wvl_center[i], self.band_means[i], self.band_stds[i]]):
ele = ElementTree.SubElement(sub, k)
ele.text = str(val)
#######################
# generate XML string #
#######################
xml_string = ElementTree.tostring(xml, encoding='unicode', pretty_print=True)
......
......@@ -291,9 +291,9 @@ class VNIR_SWIR_Stacker(object):
for cwl in wvls_overlap_all])
# apply a spectral moving average to the overlapping VNIR/SWIR band
bandidxs2average = np.array([bandidxs_overlap.min() - int((filterwidth - 1) / 2)] +
bandidxs2average = np.array([np.min(bandidxs_overlap) - int((filterwidth - 1) / 2)] +
list(bandidxs_overlap) +
[bandidxs_overlap.max() + int((filterwidth - 1) / 2)])
[np.max(bandidxs_overlap) + int((filterwidth - 1) / 2)])
data2average = data_stacked[:, :, bandidxs2average]
data_stacked[:, :, bandidxs_overlap] = mvgavg(data2average,
n=filterwidth,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment