From cf1d44f0c41177ed390497a71032a1749baad7a9 Mon Sep 17 00:00:00 2001
From: Daniel Scheffler <danschef@gfz-potsdam.de>
Date: Thu, 23 Apr 2020 00:27:16 +0200
Subject: [PATCH 1/3] Revised computation of the common VNIR/SWIR extent within
 orthorectification (fixes issue #34). All pixels that have values in VNIR or
 SWIR only are not set to nodata in the L2A output. Nodata values of masks are
 now set.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>
---
 HISTORY.rst                                   |  4 ++
 enpt/model/images/image_baseclasses.py        | 13 ++--
 enpt/options/config.py                        |  2 +-
 .../orthorectification/orthorectification.py  | 71 +++++++++++++------
 4 files changed, 63 insertions(+), 27 deletions(-)

diff --git a/HISTORY.rst b/HISTORY.rst
index 17c93c30..85dae5e3 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,10 @@
 History
 =======
 
+* Revised computation of the common VNIR/SWIR extent within orthorectification (fixes issue #34).
+* All pixels that have values in VNIR or SWIR only are not set to nodata in the L2A output (fixes issue #34).
+* Nodata values of masks are now set.
+
 0.12.3 (2020-04-21)
 -------------------
 
diff --git a/enpt/model/images/image_baseclasses.py b/enpt/model/images/image_baseclasses.py
index 4e0e44bf..e1021b09 100644
--- a/enpt/model/images/image_baseclasses.py
+++ b/enpt/model/images/image_baseclasses.py
@@ -161,7 +161,7 @@ class _EnMAP_Image(object):
 
     @mask_landwater.setter
     def mask_landwater(self, *geoArr_initArgs):
-        self._mask_landwater = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_landwater')
+        self._mask_landwater = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_landwater', nodataVal=0)
 
     @mask_landwater.deleter
     def mask_landwater(self):
@@ -177,7 +177,7 @@ class _EnMAP_Image(object):
 
     @mask_clouds.setter
     def mask_clouds(self, *geoArr_initArgs):
-        self._mask_clouds = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_clouds')
+        self._mask_clouds = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_clouds', nodataVal=0)
 
     @mask_clouds.deleter
     def mask_clouds(self):
@@ -193,7 +193,8 @@ class _EnMAP_Image(object):
 
     @mask_cloudshadow.setter
     def mask_cloudshadow(self, *geoArr_initArgs):
-        self._mask_cloudshadow = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cloudshadow')
+        self._mask_cloudshadow = \
+            self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cloudshadow', nodataVal=0)
 
     @mask_cloudshadow.deleter
     def mask_cloudshadow(self):
@@ -209,7 +210,7 @@ class _EnMAP_Image(object):
 
     @mask_haze.setter
     def mask_haze(self, *geoArr_initArgs):
-        self._mask_haze = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_haze')
+        self._mask_haze = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_haze', nodataVal=0)
 
     @mask_haze.deleter
     def mask_haze(self):
@@ -225,7 +226,7 @@ class _EnMAP_Image(object):
 
     @mask_snow.setter
     def mask_snow(self, *geoArr_initArgs):
-        self._mask_snow = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_snow')
+        self._mask_snow = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_snow', nodataVal=0)
 
     @mask_snow.deleter
     def mask_snow(self):
@@ -241,7 +242,7 @@ class _EnMAP_Image(object):
 
     @mask_cirrus.setter
     def mask_cirrus(self, *geoArr_initArgs):
-        self._mask_cirrus = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cirrus')
+        self._mask_cirrus = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cirrus', nodataVal=0)
 
     @mask_cirrus.deleter
     def mask_cirrus(self):
diff --git a/enpt/options/config.py b/enpt/options/config.py
index 4793ef94..f35c6cb0 100644
--- a/enpt/options/config.py
+++ b/enpt/options/config.py
@@ -135,7 +135,7 @@ class EnPTConfig(object):
              path to JSON file containing configuration parameters or a string in JSON format
 
         :key CPUs:
-             number of CPU cores to be used for processing (default: "None" -> use all available
+             number of CPU cores to be used for processing (default: "None" -> use all available)
 
         :key path_l1b_enmap_image:
             input path of the EnMAP L1B image to be processed
diff --git a/enpt/processors/orthorectification/orthorectification.py b/enpt/processors/orthorectification/orthorectification.py
index 55c15d87..e8342893 100644
--- a/enpt/processors/orthorectification/orthorectification.py
+++ b/enpt/processors/orthorectification/orthorectification.py
@@ -92,19 +92,23 @@ class Orthorectifier(object):
         tgt_extent = self._get_common_extent(enmap_ImageL1, tgt_epsg, enmap_grid=True)
         kw_init = dict(resamp_alg=self.cfg.ortho_resampAlg,
                        nprocs=self.cfg.CPUs,
-                       radius_of_influence=30 if not self.cfg.ortho_resampAlg == 'bilinear' else 45)
+                       # neighbours=8,
+                       radius_of_influence=30 if not self.cfg.ortho_resampAlg == 'bilinear' else 45
+                       )
         kw_trafo = dict(tgt_prj=tgt_epsg, tgt_extent=tgt_extent)
 
         # 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...')
+        enmap_ImageL1.logger.info("Orthorectifying VNIR data using '%s' resampling algorithm..."
+                                  % self.cfg.ortho_resampAlg)
         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...')
+        enmap_ImageL1.logger.info("Orthorectifying SWIR data using '%s' resampling algorithm..."
+                                  % self.cfg.ortho_resampAlg)
         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)
@@ -130,11 +134,23 @@ class Orthorectifier(object):
 
             if attr is not None:
                 enmap_ImageL1.logger.info("Orthorectifying '%s' attribute..." % attrName)
-                attr_ortho = GeoArray(*GT_2D.to_map_geometry(attr, **kw_trafo))
+                attr_ortho = GeoArray(*GT_2D.to_map_geometry(attr, **kw_trafo), nodata=attr.nodata)
                 setattr(L2_obj, attrName, attr_ortho)
 
         # TODO transform dead pixel map, quality test flags?
 
+        # set all pixels to nodata that don't have values in all bands #
+        ################################################################
+
+        enmap_ImageL1.logger.info("Setting all pixels to nodata that have values in the VNIR or the SWIR only...")
+        mask_nodata_common = np.all(np.dstack([vnir_mapgeo_gA.mask_nodata[:],
+                                               swir_mapgeo_gA.mask_nodata[:]]), axis=2)
+        L2_obj.data[~mask_nodata_common] = L2_obj.data.nodata
+
+        for attr_gA in [L2_obj.mask_landwater, L2_obj.mask_clouds, L2_obj.mask_cloudshadow, L2_obj.mask_haze,
+                        L2_obj.mask_snow, L2_obj.mask_cirrus]:
+            attr_gA[~mask_nodata_common] = attr_gA.nodata
+
         # metadata adjustments #
         ########################
 
@@ -170,23 +186,37 @@ class Orthorectifier(object):
         :param enmap_grid:
         :return:
         """
-        vnir_coords_utm = np.array([transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y)
-                                    for x, y in zip(enmap_ImageL1.vnir.detector_meta.lon_UL_UR_LL_LR,
-                                                    enmap_ImageL1.vnir.detector_meta.lat_UL_UR_LL_LR,)])
-        swir_coords_utm = np.array([transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y)
-                                    for x, y in zip(enmap_ImageL1.swir.detector_meta.lon_UL_UR_LL_LR,
-                                                    enmap_ImageL1.swir.detector_meta.lat_UL_UR_LL_LR, )])
-        all_coords_utm = np.dstack([vnir_coords_utm, swir_coords_utm])
-
-        common_extent = (all_coords_utm[:, 0, :].min(),  # xmin
-                         all_coords_utm[:, 1, :].min(),  # ymin
-                         all_coords_utm[:, 0, :].max(),  # xmax
-                         all_coords_utm[:, 1, :].max())  # ymax
-
+        # use only the first geolayer bands  # FIXME respect keystone by using the outermost coords of ALL bands?
+        V_lons, V_lats = enmap_ImageL1.meta.vnir.lons[:, :, 0], enmap_ImageL1.meta.vnir.lats[:, :, 0]
+        S_lons, S_lats = enmap_ImageL1.meta.swir.lons[:, :, 0], enmap_ImageL1.meta.swir.lats[:, :, 0]
+
+        # get Lon/Lat corner coordinates of geolayers
+        V_UL_UR_LL_LR_ll = [(V_lons[y, x], V_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, -1), (-1, 0)]]
+        S_UL_UR_LL_LR_ll = [(S_lons[y, x], S_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, -1), (-1, 0)]]
+
+        # transform them to UTM
+        V_UL_UR_LL_LR_utm = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in V_UL_UR_LL_LR_ll]
+        S_UL_UR_LL_LR_utm = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in S_UL_UR_LL_LR_ll]
+
+        # separate X and Y
+        V_X_utm, V_Y_utm = zip(*V_UL_UR_LL_LR_utm)
+        S_X_utm, S_Y_utm = zip(*S_UL_UR_LL_LR_utm)
+
+        # use inner coordinates as common extent
+        xmin_utm = max([min(V_X_utm), min(S_X_utm)])
+        ymin_utm = max([min(V_Y_utm), min(S_Y_utm)])
+        xmax_utm = min([max(V_X_utm), max(S_X_utm)])
+        ymax_utm = min([max(V_Y_utm), max(S_Y_utm)])
+        common_extent_utm = (xmin_utm, ymin_utm, xmax_utm, ymax_utm)
+
+        # move the extent to the EnMAP coordinate grid
         if enmap_grid:
-            common_extent = move_extent_to_EnMAP_grid(common_extent)
+            common_extent_utm = move_extent_to_EnMAP_grid(common_extent_utm)
+
+        enmap_ImageL1.logger.info('Computed common target extent of orthorectified image (xmin, ymin, xmax, ymax in '
+                                  'EPSG %s): %s' % (tgt_epsg, str(common_extent_utm)))
 
-        return common_extent
+        return common_extent_utm
 
 
 class VNIR_SWIR_Stacker(object):
@@ -295,7 +325,8 @@ class VNIR_SWIR_Stacker(object):
         else:
             raise ValueError(algorithm)
 
-        gA_stacked = GeoArray(data_stacked, geotransform=self.vnir.gt, projection=self.vnir.prj)
+        gA_stacked = GeoArray(data_stacked,
+                              geotransform=self.vnir.gt, projection=self.vnir.prj, nodata=self.vnir.nodata)
         gA_stacked.meta.band_meta['wavelength'] = list(wvls)
 
         return gA_stacked
-- 
GitLab


From fb749b992ec73272c65482d44faf6a646c744e69 Mon Sep 17 00:00:00 2001
From: Daniel Scheffler <danschef@gfz-potsdam.de>
Date: Thu, 23 Apr 2020 00:53:14 +0200
Subject: [PATCH 2/3] Updated HISTORY.rst.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>
---
 HISTORY.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/HISTORY.rst b/HISTORY.rst
index 85dae5e3..082ed2ab 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,10 +2,14 @@
 History
 =======
 
+0.12.4 (2020-04-??)
+-------------------
+
 * Revised computation of the common VNIR/SWIR extent within orthorectification (fixes issue #34).
 * All pixels that have values in VNIR or SWIR only are not set to nodata in the L2A output (fixes issue #34).
 * Nodata values of masks are now set.
 
+
 0.12.3 (2020-04-21)
 -------------------
 
-- 
GitLab


From 5ce9b745a4ff93afad33efac87fe87fbf284b2bb Mon Sep 17 00:00:00 2001
From: Daniel Scheffler <danschef@gfz-potsdam.de>
Date: Thu, 23 Apr 2020 12:37:06 +0200
Subject: [PATCH 3/3] Fixed computation of VNIR/SWIR overlap extent for 2D
 geolayers. This computation now also respects deviations in per-band
 geolayers due to keystone or misregistrations.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>
---
 HISTORY.rst                                   |  3 ++-
 .../orthorectification/orthorectification.py  | 19 ++++++++++++++-----
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/HISTORY.rst b/HISTORY.rst
index 082ed2ab..16ead845 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -5,7 +5,8 @@ History
 0.12.4 (2020-04-??)
 -------------------
 
-* Revised computation of the common VNIR/SWIR extent within orthorectification (fixes issue #34).
+* Revised computation of the common VNIR/SWIR extent within orthorectification (fixes issue #34). This computation now
+  also respects deviations in per-band geolayers due to keystone or misregistrations.
 * All pixels that have values in VNIR or SWIR only are not set to nodata in the L2A output (fixes issue #34).
 * Nodata values of masks are now set.
 
diff --git a/enpt/processors/orthorectification/orthorectification.py b/enpt/processors/orthorectification/orthorectification.py
index e8342893..12c46a9b 100644
--- a/enpt/processors/orthorectification/orthorectification.py
+++ b/enpt/processors/orthorectification/orthorectification.py
@@ -149,7 +149,8 @@ class Orthorectifier(object):
 
         for attr_gA in [L2_obj.mask_landwater, L2_obj.mask_clouds, L2_obj.mask_cloudshadow, L2_obj.mask_haze,
                         L2_obj.mask_snow, L2_obj.mask_cirrus]:
-            attr_gA[~mask_nodata_common] = attr_gA.nodata
+            if attr_gA is not None:
+                attr_gA[~mask_nodata_common] = attr_gA.nodata
 
         # metadata adjustments #
         ########################
@@ -186,9 +187,9 @@ class Orthorectifier(object):
         :param enmap_grid:
         :return:
         """
-        # use only the first geolayer bands  # FIXME respect keystone by using the outermost coords of ALL bands?
-        V_lons, V_lats = enmap_ImageL1.meta.vnir.lons[:, :, 0], enmap_ImageL1.meta.vnir.lats[:, :, 0]
-        S_lons, S_lats = enmap_ImageL1.meta.swir.lons[:, :, 0], enmap_ImageL1.meta.swir.lats[:, :, 0]
+        # get geolayers - 2D for dummy data format else 3D
+        V_lons, V_lats = enmap_ImageL1.meta.vnir.lons, enmap_ImageL1.meta.vnir.lats
+        S_lons, S_lats = enmap_ImageL1.meta.swir.lons, enmap_ImageL1.meta.swir.lats
 
         # get Lon/Lat corner coordinates of geolayers
         V_UL_UR_LL_LR_ll = [(V_lons[y, x], V_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, -1), (-1, 0)]]
@@ -202,7 +203,15 @@ class Orthorectifier(object):
         V_X_utm, V_Y_utm = zip(*V_UL_UR_LL_LR_utm)
         S_X_utm, S_Y_utm = zip(*S_UL_UR_LL_LR_utm)
 
-        # use inner coordinates as common extent
+        # in case of 3D geolayers, the corner coordinates have multiple values for multiple bands
+        # -> use the innermost coordinates to avoid pixels with VNIR-only/SWIR-only values due to keystone
+        if V_lons.ndim == 3:
+            V_X_utm = (V_X_utm[0].max(), V_X_utm[1].min(), V_X_utm[2].max(), V_X_utm[3].min())
+            V_Y_utm = (V_Y_utm[0].min(), V_Y_utm[1].min(), V_Y_utm[2].max(), V_Y_utm[3].max())
+            S_X_utm = (S_X_utm[0].max(), S_X_utm[1].min(), S_X_utm[2].max(), S_X_utm[3].min())
+            S_Y_utm = (S_Y_utm[0].min(), S_Y_utm[1].min(), S_Y_utm[2].max(), S_Y_utm[3].max())
+
+        # use inner coordinates of VNIR and SWIR as common extent
         xmin_utm = max([min(V_X_utm), min(S_X_utm)])
         ymin_utm = max([min(V_Y_utm), min(S_Y_utm)])
         xmax_utm = min([max(V_X_utm), max(S_X_utm)])
-- 
GitLab