Commit 2259c234 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Moved Orthorecifier._get_VNIR_SWIR_stack() to new class "VNIR_SWIR_Stacker"....


Moved Orthorecifier._get_VNIR_SWIR_stack() to new class "VNIR_SWIR_Stacker". Fixed issue that caused the Orthorectifier to always return float64 data (fixed by an update of sensormapgeo). EnMAP_Metadata_L2A_MapGeo now respects the real list of wavelengths included in the L2A product.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 04a1f122
......@@ -41,7 +41,7 @@ test_enpt:
# - cd ../../
# install sensormapgeo
- pip install sensormapgeo # TODO: remove when included in docker container
- pip install sensormapgeo>=0.2.1 # TODO: remove when included in docker container
# run nosetests
- make nosetests # test are called here
......@@ -109,7 +109,7 @@ test_enpt_install:
- cd ../../
# install sensormapgeo
- pip install sensormapgeo # TODO: remove when included in docker container
- pip install sensormapgeo>=0.2.1 # TODO: remove when included in docker container
# install enpt
- pip install -e .
......@@ -123,7 +123,6 @@ test_enpt_install:
only:
- master
- dev
- enhancement/implement_sensormapgeo
pages:
......
......@@ -627,12 +627,14 @@ class EnMAP_Metadata_L2A_MapGeo(object):
def __init__(self,
config: EnPTConfig,
meta_l1b: EnMAP_Metadata_L1B_SensorGeo,
wvls_l2a: Union[List, np.ndarray],
dims_mapgeo: Tuple[int, int, int],
logger=None):
"""EnMAP Metadata class for the metadata of the complete EnMAP L2A product in map geometry incl. VNIR and SWIR.
:param config: EnPT configuration object
:param meta_l1b: metadata object of the L1B dataset in sensor geometry
:param wvls_l2a: list of center wavelengths included in the L2A product
:param dims_mapgeo: dimensions of the EnMAP raster data in map geometry, e.g., (1024, 1000, 218)
:param logger: instance of logging.logger or subclassed
"""
......@@ -670,9 +672,9 @@ 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_vnir_plus_swir = np.hstack([self._meta_l1b.vnir.wvl_center, self._meta_l1b.swir.wvl_center])
wvls_sorted = np.array(sorted(wvls_vnir_plus_swir))
bandidx_order = np.array([np.argmin(np.abs(wvls_vnir_plus_swir - cwl)) for cwl in wvls_sorted])
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_l2a - cwl)) for cwl in wvls_sorted])
self.wvl_center = np.hstack([meta_l1b.vnir.wvl_center, meta_l1b.swir.wvl_center])[bandidx_order]
self.fwhm = np.hstack([meta_l1b.vnir.fwhm, meta_l1b.swir.fwhm])[bandidx_order]
......
......@@ -33,6 +33,7 @@ based on a pixel- and band-wise coordinate-layer (geolayer).
"""
from typing import Tuple, Union # noqa: F401
from types import SimpleNamespace
import numpy as np
from mvgavg import mvgavg
......@@ -53,7 +54,7 @@ __author__ = 'Daniel Scheffler'
class Orthorectifier(object):
def __init__(self, config: EnPTConfig = None):
def __init__(self, config: EnPTConfig):
"""Create an instance of Orthorectifier."""
self.cfg = config
......@@ -97,18 +98,24 @@ class Orthorectifier(object):
# transform VNIR and SWIR to map geometry
GeoTransformer = Geometry_Transformer if lons_vnir.ndim == 2 else Geometry_Transformer_3D
# FIXME So far, the fill value is set to 0. Is this meaningful?
enmap_ImageL1.logger.info('Orthorectifying VNIR data...')
GT_vnir = GeoTransformer(lons=lons_vnir, lats=lats_vnir, **kw_init)
vnir_mapgeo, vnir_gt, vnir_prj = GT_vnir.to_map_geometry(enmap_ImageL1.vnir.data[:], **kw_trafo)
GT_vnir = GeoTransformer(lons=lons_vnir, lats=lats_vnir, fill_value=0, **kw_init)
vnir_mapgeo_gA = GeoArray(*GT_vnir.to_map_geometry(enmap_ImageL1.vnir.data[:], **kw_trafo),
nodata=0)
enmap_ImageL1.logger.info('Orthorectifying SWIR data...')
GT_swir = GeoTransformer(lons=lons_swir, lats=lats_swir, **kw_init)
swir_mapgeo, swir_gt, swir_prj = GT_swir.to_map_geometry(enmap_ImageL1.swir.data[:], **kw_trafo)
GT_swir = GeoTransformer(lons=lons_swir, lats=lats_swir, fill_value=0, **kw_init)
swir_mapgeo_gA = GeoArray(*GT_swir.to_map_geometry(enmap_ImageL1.swir.data[:], **kw_trafo),
nodata=0)
# combine VNIR and SWIR
enmap_ImageL1.logger.info('Merging VNIR and SWIR data...')
L2_obj.data = self._get_VNIR_SWIR_stack(vnir_mapgeo, swir_mapgeo, vnir_gt, swir_gt, vnir_prj, swir_prj,
enmap_ImageL1.meta.vnir.wvl_center, enmap_ImageL1.meta.swir.wvl_center)
L2_obj.data = VNIR_SWIR_Stacker(vnir=vnir_mapgeo_gA,
swir=swir_mapgeo_gA,
vnir_wvls=enmap_ImageL1.meta.vnir.wvl_center,
swir_wvls=enmap_ImageL1.meta.swir.wvl_center
).compute_stack(algorithm=self.cfg.vswir_overlap_algorithm)
# TODO allow to set geolayer band to be used for warping of 2D arrays
GT_2D = Geometry_Transformer(lons=lons_vnir if lons_vnir.ndim == 2 else lons_vnir[:, :, 0],
......@@ -127,11 +134,11 @@ class Orthorectifier(object):
enmap_ImageL1.logger.info('Generating L2A metadata...')
L2_obj.meta = EnMAP_Metadata_L2A_MapGeo(config=self.cfg,
meta_l1b=enmap_ImageL1.meta,
wvls_l2a=L2_obj.data.meta.band_meta['wavelength'],
dims_mapgeo=L2_obj.data.shape,
logger=L2_obj.logger)
L2_obj.meta.add_band_statistics(L2_obj.data)
L2_obj.data.meta.band_meta['wavelength'] = list(L2_obj.meta.wvl_center)
L2_obj.data.meta.band_meta['bandwidths'] = list(L2_obj.meta.fwhm)
L2_obj.data.meta.global_meta['wavelength_units'] = 'nanometers'
......@@ -174,81 +181,106 @@ class Orthorectifier(object):
return common_extent
def _get_VNIR_SWIR_stack(self,
vnir_data: Union[np.ndarray, GeoArray],
swir_data: Union[np.ndarray, GeoArray],
vnir_gt: tuple,
swir_gt: tuple,
vnir_prj: str,
swir_prj: str,
wvls_vnir: np.ndarray,
wvls_swir: np.ndarray
) -> GeoArray:
"""Stack VNIR and SWIR bands with respect to their spectral overlap.
:param vnir_data:
:param swir_data:
:param vnir_gt:
:param swir_gt:
:param vnir_prj:
:param swir_prj:
:param wvls_vnir:
:param wvls_swir:
:return:
class VNIR_SWIR_Stacker(object):
def __init__(self,
vnir: GeoArray,
swir: GeoArray,
vnir_wvls: Union[list, np.ndarray],
swir_wvls: Union[list, np.ndarray])\
-> None:
"""Get an instance of VNIR_SWIR_Stacker.
:param vnir:
:param swir:
:param vnir_wvls:
:param swir_wvls:
"""
self.vnir = vnir
self.swir = swir
self.wvls = SimpleNamespace(vnir=vnir_wvls, swir=swir_wvls)
self.wvls.vswir = np.hstack([self.wvls.vnir, self.wvls.swir])
self.wvls.vswir_sorted = np.array(sorted(self.wvls.vswir))
self._validate_input()
def _validate_input(self):
if self.vnir.gt != self.swir.gt:
raise ValueError((self.vnir.gt, self.swir.gt), 'VNIR and SWIR geoinformation should be equal.')
if not prj_equal(self.vnir.prj, self.swir.prj):
raise ValueError((self.vnir.prj, self.swir.prj), 'VNIR and SWIR projection should be equal.')
def _get_stack_order_by_wvl(self) -> Tuple[np.ndarray, np.ndarray]:
"""Stack bands ordered by wavelengths."""
bandidx_order = np.array([np.argmin(np.abs(self.wvls.vswir - cwl))
for cwl in self.wvls.vswir_sorted])
return np.dstack([self.vnir[:], self.swir[:]])[:, :, bandidx_order], self.wvls.vswir_sorted
def _get_stack_average(self, filterwidth: int = 3) -> Tuple[np.ndarray, np.ndarray]:
"""Stack bands and use averaging to compute the spectral information in the VNIR/SWIR overlap.
:param filterwidth: number of bands to be included in the averaging - must be an uneven number
"""
if vnir_gt != swir_gt:
raise ValueError((vnir_gt, swir_gt), 'VNIR and SWIR geoinformation should be equal.')
if not prj_equal(vnir_prj, swir_prj):
raise ValueError((vnir_prj, swir_prj), 'VNIR and SWIR projection should be equal.')
if self.cfg.vswir_overlap_algorithm in ['order_by_wvl', 'average']:
# get band index order
wvls_vswir = np.hstack([wvls_vnir, wvls_swir])
wvls_vswir_sorted = np.array(sorted(wvls_vswir))
bandidx_order = np.array([np.argmin(np.abs(wvls_vswir - cwl))
for cwl in wvls_vswir_sorted])
# stack bands ordered by wavelengths
data_stacked = np.dstack([vnir_data, swir_data])[:, :, bandidx_order]
if self.cfg.vswir_overlap_algorithm == 'average':
# get wavelenghts and indices of overlapping bands
wvls_overlap_vnir = wvls_vnir[wvls_vnir > wvls_swir.min()]
wvls_overlap_swir = wvls_swir[wvls_swir < wvls_vnir.max()]
wvls_overlap_all = np.array(sorted(np.hstack([wvls_overlap_vnir,
wvls_overlap_swir])))
bandidxs_overlap = np.array([np.argmin(np.abs(wvls_vswir_sorted - cwl))
for cwl in wvls_overlap_all])
# apply a spectral moving average to the overlapping VNIR/SWIR bands
filterwidth = 3 # number of bands to be included in the averaging - must be an uneven number
bandidxs2average = np.array([bandidxs_overlap.min() - int((filterwidth - 1) / 2)] +
list(bandidxs_overlap) +
[bandidxs_overlap.max() + int((filterwidth - 1) / 2)])
data2average = data_stacked[:, :, bandidxs2average]
data_stacked[:, :, bandidxs_overlap] = mvgavg(data2average,
n=filterwidth,
axis=2).astype(data_stacked.dtype)
elif self.cfg.vswir_overlap_algorithm == 'vnir_only':
# get band index order
wvls_swir_cut = wvls_swir[wvls_swir > wvls_vnir.max()]
# wvls_vswir_sorted = np.hstack([wvls_vnir, wvls_swir_cut])
idx_swir_firstband = np.argmin(np.abs(wvls_swir - wvls_swir_cut.min()))
# stack bands while removing overlapping SWIR bands
data_stacked = np.dstack([vnir_data, swir_data[:, :, idx_swir_firstband:]])
elif self.cfg.vswir_overlap_algorithm == 'swir_only':
# get band index order
wvls_vnir_cut = wvls_vnir[wvls_vnir < wvls_swir.min()]
# wvls_vswir_sorted = np.hstack([wvls_vnir_cut, wvls_swir])
idx_vnir_lastband = np.argmin(np.abs(wvls_vnir - wvls_vnir_cut.max()))
# stack bands while removing overlapping VNIR bands
data_stacked = np.dstack([vnir_data[:, :, :idx_vnir_lastband], swir_data])
data_stacked = self._get_stack_order_by_wvl()[0]
# get wavelenghts and indices of overlapping bands
wvls_overlap_vnir = self.wvls.vnir[self.wvls.vnir > self.wvls.swir.min()]
wvls_overlap_swir = self.wvls.swir[self.wvls.swir < self.wvls.vnir.max()]
wvls_overlap_all = np.array(sorted(np.hstack([wvls_overlap_vnir,
wvls_overlap_swir])))
bandidxs_overlap = np.array([np.argmin(np.abs(self.wvls.vswir_sorted - cwl))
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)] +
list(bandidxs_overlap) +
[bandidxs_overlap.max() + int((filterwidth - 1) / 2)])
data2average = data_stacked[:, :, bandidxs2average]
data_stacked[:, :, bandidxs_overlap] = mvgavg(data2average,
n=filterwidth,
axis=2).astype(data_stacked.dtype)
return data_stacked, self.wvls.vswir_sorted
def _get_stack_vnir_only(self) -> Tuple[np.ndarray, np.ndarray]:
"""Stack bands while removing overlapping SWIR bands."""
wvls_swir_cut = self.wvls.swir[self.wvls.swir > self.wvls.vnir.max()]
wvls_vswir_sorted = np.hstack([self.wvls.vnir, wvls_swir_cut])
idx_swir_firstband = np.argmin(np.abs(self.wvls.swir - wvls_swir_cut.min()))
return np.dstack([self.vnir[:], self.swir[:, :, idx_swir_firstband:]]), wvls_vswir_sorted
def _get_stack_swir_only(self) -> Tuple[np.ndarray, np.ndarray]:
"""Stack bands while removing overlapping VNIR bands."""
wvls_vnir_cut = self.wvls.vnir[self.wvls.vnir < self.wvls.swir.min()]
wvls_vswir_sorted = np.hstack([wvls_vnir_cut, self.wvls.swir])
idx_vnir_lastband = np.argmin(np.abs(self.wvls.vnir - wvls_vnir_cut.max()))
return np.dstack([self.vnir[:, :, :idx_vnir_lastband], self.swir[:]]), wvls_vswir_sorted
def compute_stack(self, algorithm: str) -> GeoArray:
"""Stack VNIR and SWIR bands with respect to their spectral overlap.
:param algorithm: 'order_by_wvl': keep spectral bands unchanged, order bands by wavelength
'average': average the spectral information within the overlap
'vnir_only': only use the VNIR bands (cut overlapping SWIR bands)
'swir_only': only use the SWIR bands (cut overlapping VNIR bands)
:return: the stacked data cube as GeoArray instance
"""
if algorithm == 'order_by_wvl':
data_stacked, wvls = self._get_stack_order_by_wvl()
elif algorithm == 'average':
data_stacked, wvls = self._get_stack_average()
elif algorithm == 'vnir_only':
data_stacked, wvls = self._get_stack_vnir_only()
elif algorithm == 'swir_only':
data_stacked, wvls = self._get_stack_swir_only()
else:
raise ValueError(self.cfg.vswir_overlap_algorithm)
raise ValueError(algorithm)
gA_stacked = GeoArray(data_stacked, geotransform=self.vnir.gt, projection=self.vnir.prj)
gA_stacked.meta.band_meta['wavelength'] = list(wvls)
return GeoArray(data_stacked, geotransform=vnir_gt, projection=vnir_prj)
return gA_stacked
......@@ -55,7 +55,7 @@ class Radiometric_Transformer(object):
self.earthSunDist = config.path_earthSunDist # path of model for earth sun distance
def transform_TOARad2TOARef(self, enmap_ImageL1: EnMAPL1Product_SensorGeo):
"""Transform top-of-atmosphere radiance to top-of-atmosphere reflectance.
"""Transform top-of-atmosphere radiance to top-of-atmosphere reflectance (16-bit signed-integer).
NOTE: The following formula is used:
toaRef = (scale_factor * math.pi * toaRad * earthSunDist**2) /
......
......@@ -128,7 +128,8 @@ class Geometry_Transformer_3D(SensorMapGeometryTransformer3D):
path_or_geoarray_sensorgeo: Union[str, GeoArray, np.ndarray],
tgt_prj: Union[str, int] = None,
tgt_extent: Tuple[float, float, float, float] = None,
tgt_res: Tuple[float, float] = None):
tgt_res: Tuple[float, float] = None
) -> Tuple[np.ndarray, tuple, str]:
data_sensorgeo = GeoArray(path_or_geoarray_sensorgeo)
if data_sensorgeo.is_map_geo:
......
......@@ -44,7 +44,7 @@ with open("enpt/version.py", encoding='utf-8') as version_file:
exec(version_file.read(), version)
requirements = [ # put package requirements here
'numpy', 'pandas', 'scipy', 'geoarray>=0.8.11', 'py_tools_ds>=0.14.23', 'arosics>=0.9.2', 'sensormapgeo',
'numpy', 'pandas', 'scipy', 'geoarray>=0.8.11', 'py_tools_ds>=0.14.23', 'arosics>=0.9.2', 'sensormapgeo>=0.2.1',
'cerberus', 'jsmin', 'matplotlib', 'tqdm', 'utm', 'lxml', 'numpy-indexed', 'mvgavg'
# 'sicor', # pip install git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
]
......
......@@ -26,7 +26,7 @@ dependencies:
- bokeh
# sensormapgeo
- pyresample
- pyresample>=1.11.0
# arosics
- pyfftw
......@@ -46,7 +46,7 @@ dependencies:
- scipy
- geoarray>=0.8.11
- py_tools_ds>=0.14.25
- sensormapgeo
- sensormapgeo>=0.2.1
- cerberus
- jsmin
- tqdm
......
Supports Markdown
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