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