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

Merge branch 'enhancement/improve_VSWIR_spectral_fusion' into 'master'

Enhancement/improve vswir spectral fusion

See merge request !36
parents 3fc4caf7 60c9dade
Pipeline #7473 passed with stages
in 57 minutes and 34 seconds
......@@ -197,5 +197,8 @@ fabric.properties
.gitignore
.idea/**
.idea/
tests/data/test_outputs/
tests/linting/flake8.log
tests/linting/pycodestyle.log.log
tests/linting/pydocstyle.log.log
......@@ -40,8 +40,10 @@ test_enpt:
# - python setup.py install
# - cd ../../
# install sensormapgeo
- pip install sensormapgeo # TODO: remove when included in docker container
# install sensormapgeo, mvgavg
# TODO: remove when included in docker container
- pip install sensormapgeo>=0.2.1
- pip install mvgavg
# run nosetests
- make nosetests # test are called here
......@@ -108,9 +110,6 @@ test_enpt_install:
# - make install
- cd ../../
# install sensormapgeo
- pip install sensormapgeo # TODO: remove when included in docker container
# install enpt
- pip install -e .
- cd ..
......@@ -123,7 +122,6 @@ test_enpt_install:
only:
- master
- dev
- enhancement/implement_sensormapgeo
pages:
......
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$">
<excludeFolder url="file://$MODULE_DIR$/venv" />
</content>
<orderEntry type="jdk" jdkName="Python 3.6.0 (~/python_GFZ/python/bin/python3.6)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="PyDocumentationSettings">
<option name="renderExternalDocumentation" value="true" />
</component>
<component name="TestRunnerService">
<option name="projectConfiguration" value="Nosetests" />
<option name="PROJECT_TEST_RUNNER" value="Nosetests" />
</component>
</module>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="Encoding" addBOMForNewFiles="with NO BOM" />
</project>
\ No newline at end of file
<component name="InspectionProjectProfileManager">
<profile version="1.0">
<option name="myName" value="Project Default" />
<inspection_tool class="PyCompatibilityInspection" enabled="true" level="WARNING" enabled_by_default="true">
<option name="ourVersions">
<value>
<list size="3">
<item index="0" class="java.lang.String" itemvalue="3.5" />
<item index="1" class="java.lang.String" itemvalue="3.6" />
<item index="2" class="java.lang.String" itemvalue="3.7" />
</list>
</value>
</option>
</inspection_tool>
<inspection_tool class="PyPackageRequirementsInspection" enabled="true" level="WARNING" enabled_by_default="true">
<option name="ignoredPackages">
<value>
<list size="11">
<item index="0" class="java.lang.String" itemvalue="pyhdf" />
<item index="1" class="java.lang.String" itemvalue="osr" />
<item index="2" class="java.lang.String" itemvalue="ephem" />
<item index="3" class="java.lang.String" itemvalue="sklearn" />
<item index="4" class="java.lang.String" itemvalue="sicor" />
<item index="5" class="java.lang.String" itemvalue="py-tools-ds" />
<item index="6" class="java.lang.String" itemvalue="geoarray" />
<item index="7" class="java.lang.String" itemvalue="S2SCAPEM" />
<item index="8" class="java.lang.String" itemvalue="S2MSI" />
<item index="9" class="java.lang.String" itemvalue="arosics" />
<item index="10" class="java.lang.String" itemvalue="py_tools_ds" />
</list>
</value>
</option>
</inspection_tool>
<inspection_tool class="PyPep8NamingInspection" enabled="false" level="WEAK WARNING" enabled_by_default="false" />
<inspection_tool class="SpellCheckingInspection" enabled="false" level="TYPO" enabled_by_default="false">
<option name="processCode" value="true" />
<option name="processLiterals" value="true" />
<option name="processComments" value="true" />
</inspection_tool>
</profile>
</component>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="JsonSchemaMappingsProjectConfiguration">
<state>
<map>
<entry key="StyleLint">
<value>
<SchemaInfo>
<option name="name" value="StyleLint" />
<option name="relativePathToSchema" value="http://json.schemastore.org/stylelintrc" />
<option name="applicationDefined" value="true" />
<option name="patterns">
<list>
<Item>
<option name="path" value="tests/gitlab_CI_docker/context/environment_enpt.yml" />
</Item>
</list>
</option>
</SchemaInfo>
</value>
</entry>
</map>
</state>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="JavaScriptSettings">
<option name="languageLevel" value="ES6" />
</component>
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.6.0 (~/python_GFZ/python/bin/python3.6)" project-jdk-type="Python SDK" />
<component name="PythonCompatibilityInspectionAdvertiser">
<option name="version" value="3" />
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/EnPT.iml" filepath="$PROJECT_DIR$/.idea/EnPT.iml" />
</modules>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="PySciProjectComponent">
<option name="PY_SCI_VIEW" value="true" />
<option name="PY_SCI_VIEW_SUGGESTED" value="true" />
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="$PROJECT_DIR$" vcs="Git" />
</component>
</project>
\ No newline at end of file
......@@ -123,6 +123,9 @@ def get_enpt_argparser():
"('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')")
add('--ortho_resampAlg', type=str, default='bilinear',
help="Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')")
add('--vswir_overlap_algorithm', type=str, default='swir_only',
help="Algorithm specifying how to deal with the spectral bands in the VNIR/SWIR spectral overlap region "
"('order_by_wvl', 'average', 'vnir_only', 'swir_only')")
# link parser to run function
parser.set_defaults(func=run_job)
......
......@@ -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_sorted - cwl)) for cwl in wvls_l2a])
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]
......
......@@ -117,6 +117,7 @@ config_for_testing_dlr = dict(
enable_ice_retrieval=False,
CPUs=1,
ortho_resampAlg='gauss',
vswir_overlap_algorithm='swir_only'
)
......@@ -293,8 +294,9 @@ class EnPTConfig(object):
self.deadpix_P_interp_spectral = gp('deadpix_P_interp_spectral')
self.deadpix_P_interp_spatial = gp('deadpix_P_interp_spatial')
# orthorectification
# orthorectification / VSWIR fusion
self.ortho_resampAlg = gp('ortho_resampAlg')
self.vswir_overlap_algorithm = gp('vswir_overlap_algorithm')
#########################
# validate final config #
......
......@@ -62,7 +62,13 @@
},
"orthorectification": {
"resamp_alg": "bilinear" /*Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')*/
"resamp_alg": "bilinear", /*Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')*/
"vswir_overlap_algorithm": "swir_only" /*Algorithm how to output the spectral bands in the VNIR/SWIR spectral overlap region
'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)
*/
}
}
}
\ No newline at end of file
......@@ -106,6 +106,8 @@ enpt_schema_input = dict(
type='dict', required=False,
schema=dict(
resamp_alg=dict(type='string', required=False, allowed=['nearest', 'bilinear', 'gauss']),
vswir_overlap_algorithm=dict(type='string', required=False,
allowed=['order_by_wvl', 'average', 'vnir_only', 'swir_only']),
))
))
)
......@@ -157,6 +159,7 @@ parameter_mapping = dict(
# processors > orthorectification
ortho_resampAlg=('processors', 'orthorectification', 'resamp_alg'),
vswir_overlap_algorithm=('processors', 'orthorectification', 'vswir_overlap_algorithm'),
)
......
......@@ -32,10 +32,11 @@ EnPT module 'orthorectification' for transforming an EnMAP image from sensor to
based on a pixel- and band-wise coordinate-layer (geolayer).
"""
from .orthorectification import Orthorectifier
from .orthorectification import Orthorectifier, VNIR_SWIR_Stacker
__author__ = 'Daniel Scheffler'
__all__ = [
'__author__',
'Orthorectifier'
'Orthorectifier',
'VNIR_SWIR_Stacker'
]
......@@ -32,9 +32,11 @@ EnPT module 'orthorectification' for transforming an EnMAP image from sensor to
based on a pixel- and band-wise coordinate-layer (geolayer).
"""
from typing import Tuple # noqa: F401
from typing import Tuple, Union # noqa: F401
from types import SimpleNamespace
import numpy as np
from mvgavg import mvgavg
from geoarray import GeoArray
from py_tools_ds.geo.coord_trafo import transform_any_prj
from py_tools_ds.geo.projection import EPSG2WKT, prj_equal
......@@ -52,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
......@@ -96,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],
......@@ -126,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'
......@@ -150,7 +158,7 @@ class Orthorectifier(object):
enmap_grid: bool = True) -> Tuple[float, float, float, float]:
"""Get common target extent for VNIR and SWIR.
:param enmap_ImageL1:
:para enmap_ImageL1:
:param tgt_epsg:
:param enmap_grid:
:return:
......@@ -173,22 +181,114 @@ class Orthorectifier(object):
return common_extent
@staticmethod
def _get_VNIR_SWIR_stack(vnir_data, swir_data, vnir_gt, swir_gt, vnir_prj, swir_prj, wvls_vnir, wvls_swir):
"""Stack VNIR and SWIR bands with respect to their spectral overlap."""
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.')
# get band index order
wvls_vnir_plus_swir = np.hstack([wvls_vnir, wvls_swir])
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])
# stack bands ordered by wavelengths
data_stacked = np.dstack([vnir_data, swir_data])[:, :, bandidx_order]
# TODO implement correction for VNIR/SWIR spectral jump
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.
return GeoArray(data_stacked, geotransform=vnir_gt, projection=vnir_prj)
: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.')
if self.vnir.bands != len(self.wvls.vnir):
raise ValueError("The number of VNIR bands must be equal to the number of elements in 'vnir_wvls': "
"%d != %d" % (self.vnir.bands, len(self.wvls.vnir)))
if self.swir.bands != len(self.wvls.swir):
raise ValueError("The number of SWIR bands must be equal to the number of elements in 'swir_wvls': "
"%d != %d" % (self.swir.bands, len(self.wvls.swir)))
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
"""
# FIXME this has to respect nodata values - especially for pixels where one detector has no 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 + 1], 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
"""
# TODO: This should also set an output nodata value.
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(algorithm)
gA_stacked = GeoArray(data_stacked, geotransform=self.vnir.gt, projection=self.vnir.prj)
gA_stacked.meta.band_meta['wavelength'] = list(wvls)
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:
......
......@@ -13,3 +13,4 @@ utm
git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
lxml
numpy-indexed
mvgavg
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