Added first versions of additional algorithms to deal with the spectral overlap of VNIR and SWIR.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

bin/enpt_cli.py

@@ -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='average',
+        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)

enpt/options/config.py

@@ -293,8 +293,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 #

enpt/options/options_default.json

@@ -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": "average" /*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

enpt/options/options_schema.py

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

enpt/processors/orthorectification/orthorectification.py

@@ -32,9 +32,10 @@ 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
 
 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
@@ -150,7 +151,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 +174,81 @@ 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."""
+    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:
+        """
         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
+        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])
+
+        else:
+            raise ValueError(self.cfg.vswir_overlap_algorithm)
 
         return GeoArray(data_stacked, geotransform=vnir_gt, projection=vnir_prj)

requirements.txt

@@ -13,3 +13,4 @@ utm
 git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
 lxml
 numpy-indexed
+mvgavg

setup.py

@@ -45,7 +45,7 @@ with open("enpt/version.py", encoding='utf-8') as version_file:
 
 requirements = [  # put package requirements here
     'numpy', 'pandas', 'scipy', 'geoarray>=0.8.11', 'py_tools_ds>=0.14.23', 'arosics>=0.9.2', 'sensormapgeo',
-    'cerberus', 'jsmin', 'matplotlib', 'tqdm', 'utm', 'lxml', 'numpy-indexed'
+    'cerberus', 'jsmin', 'matplotlib', 'tqdm', 'utm', 'lxml', 'numpy-indexed', 'mvgavg'
     # 'sicor', # pip install git+https://gitext.gfz-potsdam.de/EnMAP/sicor.git
 ]
 

tests/gitlab_CI_docker/context/environment_enpt.yml

@@ -52,6 +52,7 @@ dependencies:
     - tqdm
     - utm
     - numpy-indexed
+    - mvgavg
     - sphinx-argparse
     - pycodestyle
     - flake8