diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a6b932a51e8432b7cb636800ca082d04e7b7a877..d33908a484fe8db81c5905c071d875ee77003e40 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -17,7 +17,6 @@ test_enpt:
# update some dependencies
# - pip install 'geoarray>=0.15.8'
# - pip install 'py_tools_ds>=0.14.23'
- - mamba install "sicor>=0.19.1"
# install sicor and perhaps switch sicor branch
# - rm -rf context/sicor
diff --git a/HISTORY.rst b/HISTORY.rst
index 46d1eb7f2a86661e753400d9ad8eb96729669524..5fc40733c81dbe98d89b2fe36a78b8773fc2a9fa 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,14 @@
History
=======
+0.20.0 (coming soon)
+--------------------
+
+* !103: Revised the EnPT orthorectification to support the completely re-implemented sensormapgeo package (v1.0.0),
+ which leads to 30-50x faster orthorectification and adds support for all GDAL-supported resampling techniques.
+* !103: Fixed #58 (Holes in the orthorectified image in case of geographic target projection and bilinear resampling).
+
+
0.19.7 (2024-01-04)
-------------------
diff --git a/enpt/cli.py b/enpt/cli.py
index c044bda66c26655bde40ccd9ef2561fb0b45df4a..cc4f30cf17361d9c0c2c39938f53cf3c381f94f9 100644
--- a/enpt/cli.py
+++ b/enpt/cli.py
@@ -140,7 +140,8 @@ def get_enpt_argparser():
help="Spatial interpolation algorithm to be used during dead pixel correction "
"('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')")
add('--ortho_resampAlg', type=str, default='bilinear',
- help="Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')")
+ help="Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss', 'cubic', 'cubic_spline', "
+ "'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')")
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')")
diff --git a/enpt/model/images/image_baseclasses.py b/enpt/model/images/image_baseclasses.py
index 982dd575f80e480216e1621794677a3e99b3b65e..d0be043ebcc1cbb1229f23f9f86647fd954d277f 100644
--- a/enpt/model/images/image_baseclasses.py
+++ b/enpt/model/images/image_baseclasses.py
@@ -124,7 +124,7 @@ class _EnMAP_Image(object):
self.data = '/path/to/image.tif' # sets self.data to GeoArray('/path/to/image.tif')
- Link a numpy.ndarray instance with self.data (remaining attributes like geocoding, projection, etc.
- are copied from the previous self.data attribute.
+ are copied from the previous self.data attribute).
self.data = numpy.array([[1,2,3],[4,5,6]])
- Set self.data to an existing instance of GeoArray
@@ -197,7 +197,7 @@ class _EnMAP_Image(object):
@property
def mask_cloudshadow(self) -> GeoArray:
- """Return the cloud shadow mask (0=no cloud shadow, 1=cloud shadow)..
+ """Return the cloud shadow mask (0=no cloud shadow, 1=cloud shadow).
:return: geoarray.GeoArray
"""
@@ -230,7 +230,7 @@ class _EnMAP_Image(object):
@property
def mask_snow(self) -> GeoArray:
- """Return the snow mask (0=no snow, 1=snow)..
+ """Return the snow mask (0=no snow, 1=snow).
:return: geoarray.GeoArray
"""
@@ -246,7 +246,7 @@ class _EnMAP_Image(object):
@property
def mask_cirrus(self) -> GeoArray:
- """Return the cirrus mask (0=none, 1=thin, 2=medium, 3=thick)..
+ """Return the cirrus mask (0=none, 1=thin, 2=medium, 3=thick).
:return: geoarray.GeoArray
"""
diff --git a/enpt/model/images/images_sensorgeo.py b/enpt/model/images/images_sensorgeo.py
index f4388d2c12c7f0c410468169c704496adccec745..4cc7120638c9d9bd50bf49746f72141b25946da9 100644
--- a/enpt/model/images/images_sensorgeo.py
+++ b/enpt/model/images/images_sensorgeo.py
@@ -372,9 +372,11 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image):
src_lats: np.ndarray,
src_epsg: int,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> np.ndarray:
- """Transform the given input raster from SWIR to VNIR or from SWIR to VNIR sensor geometry.
+ """Transform the given input raster from VNIR to SWIR sensor geometry or vice-versa.
NOTE:
@@ -391,9 +393,12 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image):
:param src_lats: geolayer latitudes corresponding to the input array
(same nbands like the source detector)
:param src_epsg: projection EPSG code of the source array
- :param resamp_alg: resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorithm ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case a 3D array
in the dimension of the source detector is passed (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
:return:
"""
detN = self.detector_name
@@ -431,7 +436,8 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image):
res_vnir=(30, 30),
res_swir=(30, 30),
resamp_alg=resamp_alg,
- # radius_of_influence=45,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata,
nprocs=self.cfg.CPUs
)
if detN == 'VNIR':
@@ -478,7 +484,9 @@ class EnMAP_VNIR_SensorGeo(EnMAP_Detector_SensorGeo):
swir_lats: np.ndarray,
swir_epsg: int,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> np.ndarray:
"""Transform the given SWIR sensor-geometry raster array into VNIR sensor geometry.
@@ -486,12 +494,15 @@ class EnMAP_VNIR_SensorGeo(EnMAP_Detector_SensorGeo):
:param swir_lons: longitude geolayer array of the SWIR
:param swir_lats: latitude geolayer array of the SWIR
:param swir_epsg: EPSG code of the SWIR when transformed to map geometry
- :param resamp_alg: resampling algorith ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorith ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case a 3D array
in the dimension of the SWIR detector is passed (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
"""
return self._transform_raster_geometry_from_other_detector(
- array_swirsensorgeo, swir_lons, swir_lats, swir_epsg, resamp_alg, respect_keystone)
+ array_swirsensorgeo, swir_lons, swir_lats, swir_epsg, resamp_alg, respect_keystone, src_nodata, tgt_nodata)
class EnMAP_SWIR_SensorGeo(EnMAP_Detector_SensorGeo):
@@ -513,7 +524,9 @@ class EnMAP_SWIR_SensorGeo(EnMAP_Detector_SensorGeo):
vnir_lats: np.ndarray,
vnir_epsg: int,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> np.ndarray:
"""Transform the given VNIR sensor-geometry raster array into SWIR sensor geometry.
@@ -521,12 +534,15 @@ class EnMAP_SWIR_SensorGeo(EnMAP_Detector_SensorGeo):
:param vnir_lons: longitude geolayer array of the VNIR
:param vnir_lats: latitude geolayer array of the VNIR
:param vnir_epsg: EPSG code of the VNIR when transformed to map geometry
- :param resamp_alg: resampling algorith ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorith ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case a 3D array
in the dimension of the VNIR detector is passed (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
"""
return self._transform_raster_geometry_from_other_detector(
- array_vnirsensorgeo, vnir_lons, vnir_lats, vnir_epsg, resamp_alg, respect_keystone)
+ array_vnirsensorgeo, vnir_lons, vnir_lats, vnir_epsg, resamp_alg, respect_keystone, src_nodata, tgt_nodata)
class EnMAPL1Product_SensorGeo(object):
@@ -798,14 +814,19 @@ class EnMAPL1Product_SensorGeo(object):
def transform_vnir_to_swir_raster(self,
array_vnirsensorgeo: np.ndarray,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> np.ndarray:
"""Transform the given array from VNIR into SWIR sensor geometry.
:param array_vnirsensorgeo: raster array in VNIR sensor geometry to be transformed into SWIR sensor geometry
- :param resamp_alg: resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorithm ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case a 3D array
in the dimension of the VNIR detector is passed (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
"""
if self.meta.vnir.lons is None or self.meta.vnir.lats is None or \
self.meta.swir.lons is None or self.meta.swir.lats is None:
@@ -817,18 +838,26 @@ class EnMAPL1Product_SensorGeo(object):
vnir_lats=self.meta.vnir.lats,
vnir_epsg=self.meta.vnir.epsg_ortho,
resamp_alg=resamp_alg,
- respect_keystone=respect_keystone)
+ respect_keystone=respect_keystone,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata
+ )
def transform_swir_to_vnir_raster(self, array_swirsensorgeo: np.ndarray,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> np.ndarray:
"""Transform the given array from SWIR into VNIR sensor geometry.
:param array_swirsensorgeo: raster array in SWIR sensor geometry to be transformed into VNIR sensor geometry
- :param resamp_alg: resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorithm ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case a 3D array
in the dimension of the VNIR detector is passed (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
"""
if self.meta.vnir.lons is None or self.meta.vnir.lats is None or \
self.meta.swir.lons is None or self.meta.swir.lats is None:
@@ -840,18 +869,26 @@ class EnMAPL1Product_SensorGeo(object):
swir_lats=self.meta.swir.lats,
swir_epsg=self.meta.swir.epsg_ortho,
resamp_alg=resamp_alg,
- respect_keystone=respect_keystone)
+ respect_keystone=respect_keystone,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata
+ )
def set_SWIRattr_with_transformedVNIRattr(self, attrName: str,
resamp_alg: str = 'nearest',
- respect_keystone: bool = False
+ respect_keystone: bool = False,
+ src_nodata: int = None,
+ tgt_nodata: int = None
) -> None:
"""Set the specified SWIR raster attribute with a VNIR attribute transformed to SWIR sensor geometry.
:param attrName: name of the attribute to be set
- :param resamp_alg: resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
+ :param resamp_alg: resampling algorithm ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:param respect_keystone: whether to use the full geoarray (all bands) in case the attribute
to be transformed is 'data' (default: False)
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
"""
self.logger.info("Transforming the '%s' attribute from VNIR to SWIR sensor geometry." % attrName)
@@ -862,7 +899,10 @@ class EnMAPL1Product_SensorGeo(object):
attr_transformed = self.transform_vnir_to_swir_raster(array_vnirsensorgeo=np.array(vnir_rasterAttr),
resamp_alg=resamp_alg,
- respect_keystone=respect_keystone)
+ respect_keystone=respect_keystone,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata,
+ )
setattr(self.swir, attrName, attr_transformed)
def run_AC(self):
diff --git a/enpt/model/metadata/metadata_sensorgeo.py b/enpt/model/metadata/metadata_sensorgeo.py
index 3d7b4c7376999afd1964f5dbfd494107a1ca9062..a95da215ed2ac49fc320095a490ba4317cd56c79 100644
--- a/enpt/model/metadata/metadata_sensorgeo.py
+++ b/enpt/model/metadata/metadata_sensorgeo.py
@@ -105,8 +105,9 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
self.epsg_ortho: Optional[int] = None # EPSG code of the orthorectified image
self.rpc_coeffs: OrderedDict = OrderedDict() # RPC coefficients for geolayer computation
self.ll_mapPoly: Optional[Polygon] = None # footprint polygon in longitude/latitude map coordinates
- self.lats: Optional[np.ndarray] = None # 2D array of latitude coordinates according to given lon/lat sampling
- self.lons: Optional[np.ndarray] = None # 2D array of longitude coordinates according to given lon/lat sampling
+ self.lats: Optional[np.ndarray] = None # 2D/3D array of latitude coordinates
+ self.lons: Optional[np.ndarray] = None # 2D/3D array of longitude coordinates
+ self.geolayer_has_keystone: Optional[bool] = None # indicates if lon/lat geolayer considers keystone (3D array)
self.unit: str = '' # radiometric unit of pixel values
self.unitcode: str = '' # code of radiometric unit
self.preview_wvls: Optional[List[float]] = None # wavelengths to be used for quicklook images
@@ -267,6 +268,7 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
self.lats = self.interpolate_corners(*self.lat_UL_UR_LL_LR, self.ncols, self.nrows)
self.lons = self.interpolate_corners(*self.lon_UL_UR_LL_LR, self.ncols, self.nrows)
+ self.geolayer_has_keystone = False
self.preview_wvls = np.array([self.wvl_center[i] for i in self.preview_bands])
# drop bad bands from metadata if desired
@@ -427,13 +429,17 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
def compute_geolayer_for_cube(self):
self.logger.info('Computing %s geolayer...' % self.detector_name)
- lons, lats = \
- RPC_3D_Geolayer_Generator(rpc_coeffs_per_band=self.rpc_coeffs,
- elevation=self.cfg.path_dem if self.cfg.path_dem else self.cfg.average_elevation,
- enmapIm_cornerCoords=tuple(zip(self.lon_UL_UR_LL_LR, self.lat_UL_UR_LL_LR)),
- enmapIm_dims_sensorgeo=(self.nrows, self.ncols),
- CPUs=self.cfg.CPUs)\
- .compute_geolayer()
+ GeolayerGen = \
+ RPC_3D_Geolayer_Generator(
+ rpc_coeffs_per_band=self.rpc_coeffs,
+ elevation=self.cfg.path_dem if self.cfg.path_dem else self.cfg.average_elevation,
+ enmapIm_cornerCoords=tuple(zip(self.lon_UL_UR_LL_LR, self.lat_UL_UR_LL_LR)),
+ enmapIm_dims_sensorgeo=(self.nrows, self.ncols),
+ CPUs=self.cfg.CPUs
+ )
+ lons, lats = GeolayerGen.compute_geolayer()
+
+ self.geolayer_has_keystone = len(GeolayerGen.bandgroups_with_unique_rpc_coeffs) > 1
return lons, lats
diff --git a/enpt/options/config.py b/enpt/options/config.py
index a5c11c7d381166afaad7145ea2c6f4940339fd5f..b9671cce5b0997939acc1dee16dffbb26aa0bbb2 100644
--- a/enpt/options/config.py
+++ b/enpt/options/config.py
@@ -97,7 +97,7 @@ config_for_testing_water = dict(
enable_keystone_correction=False,
enable_vnir_swir_coreg=False,
n_lines_to_append=None,
- ortho_resampAlg='gauss',
+ ortho_resampAlg='bilinear',
run_deadpix_P=True,
run_smile_P=False,
scale_factor_boa_ref=10000,
@@ -190,7 +190,7 @@ config_for_testing_dlr = dict(
enable_ac=True,
mode_ac='land',
CPUs=32,
- ortho_resampAlg='gauss',
+ ortho_resampAlg='bilinear',
vswir_overlap_algorithm='swir_only'
)
@@ -326,7 +326,8 @@ class EnPTConfig(object):
('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')
:key ortho_resampAlg:
- Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')
+ Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
:key target_projection_type:
Projection type of the raster output files ('UTM', 'Geographic') (default: 'UTM')
@@ -466,12 +467,6 @@ class EnPTConfig(object):
self.target_coord_grid if self.target_coord_grid else \
enmap_coordinate_grid_utm if self.target_projection_type == 'UTM' else None
- # bug warning regarding holes in bilinear resampling output
- if self.target_projection_type == 'Geographic' and self.ortho_resampAlg == 'bilinear':
- warnings.warn("There is currently a bug that causes holes in the bilinear resampling results if the "
- "target projection is 'Geographic'. It is recommended to use 'nearest' or 'gauss' instead.",
- UserWarning)
-
@staticmethod
def absPath(path):
return path if not path or os.path.isabs(path) else os.path.abspath(os.path.join(path_enptlib, path))
diff --git a/enpt/options/options_default.json b/enpt/options/options_default.json
index 2cd19cd91b53988e9403bba93a8f57d6dd4697c5..1096e7c821841662cb6ed7b22c8702b79307de6b 100644
--- a/enpt/options/options_default.json
+++ b/enpt/options/options_default.json
@@ -90,7 +90,9 @@
},
"orthorectification": {
- "resamp_alg": "gauss", /*Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')*/
+ "resamp_alg": "bilinear", /*Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss',
+ 'cubic', 'cubic_spline', 'lanczos', 'average', 'mode', 'max', 'min', 'med',
+ 'q1', 'q3')*/
"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
diff --git a/enpt/options/options_schema.py b/enpt/options/options_schema.py
index 45b0bfbf409d49ee2f962abccd126cf737110eeb..3b832b35696d37a03a9f9dca7cd9a060a242a782 100644
--- a/enpt/options/options_schema.py
+++ b/enpt/options/options_schema.py
@@ -114,7 +114,9 @@ enpt_schema_input = dict(
orthorectification=dict(
type='dict', required=False,
schema=dict(
- resamp_alg=dict(type='string', required=False, allowed=['nearest', 'bilinear', 'gauss']),
+ resamp_alg=dict(type='string', required=False, allowed=['nearest', 'bilinear', 'gauss', 'cubic',
+ 'cubic_spline', 'lanczos', 'average',
+ 'mode', 'max', 'min', 'med', 'q1', 'q3']),
vswir_overlap_algorithm=dict(type='string', required=False,
allowed=['order_by_wvl', 'average', 'vnir_only', 'swir_only']),
target_projection_type=dict(type='string', required=False, allowed=['UTM', 'Geographic']),
diff --git a/enpt/processors/atmospheric_correction/atmospheric_correction.py b/enpt/processors/atmospheric_correction/atmospheric_correction.py
index 128e49c7af9bb6ef80df8bf18c35393dc7cb748d..a75688c995e5d6c3a4ca2d674b2e009ac70ee67d 100644
--- a/enpt/processors/atmospheric_correction/atmospheric_correction.py
+++ b/enpt/processors/atmospheric_correction/atmospheric_correction.py
@@ -292,7 +292,7 @@ class AtmosphericCorrector(object):
:return: atmospherically corrected output EnMAP image containing BOA reflectance / water leaving reflectance
(an instance EnMAPL1Product_SensorGeo)
"""
- enmap_ImageL1.set_SWIRattr_with_transformedVNIRattr('mask_landwater')
+ enmap_ImageL1.set_SWIRattr_with_transformedVNIRattr('mask_landwater', src_nodata=0, tgt_nodata=0)
enmap_ImageL1.logger.info(
f"Starting atmospheric correction for VNIR and SWIR detector in '{self.cfg.mode_ac}' mode. "
diff --git a/enpt/processors/dem_preprocessing/dem_preprocessing.py b/enpt/processors/dem_preprocessing/dem_preprocessing.py
index e798dbedd88d25f1ec90ea81909684430f261592..1bccf7c3074c172bf90b499d2f3c0cc7dd601547 100644
--- a/enpt/processors/dem_preprocessing/dem_preprocessing.py
+++ b/enpt/processors/dem_preprocessing/dem_preprocessing.py
@@ -127,7 +127,7 @@ class DEM_Processor(object):
def to_sensor_geometry(self,
lons: np.ndarray,
lats: np.ndarray):
- GT = Geometry_Transformer(lons=lons, lats=lats, nprocs=self.CPUs)
+ GT = Geometry_Transformer(lons=lons, lats=lats, backend='gdal', resamp_alg='bilinear', nprocs=self.CPUs)
data_sensorgeo = GT.to_sensor_geometry(self.dem)
return GeoArray(data_sensorgeo)
diff --git a/enpt/processors/orthorectification/orthorectification.py b/enpt/processors/orthorectification/orthorectification.py
index aad2a8c0c28d4b1065077e2b7c081bb327cc919d..a43a042fdd76585e61cf7ab71a0259c03f60dd3b 100644
--- a/enpt/processors/orthorectification/orthorectification.py
+++ b/enpt/processors/orthorectification/orthorectification.py
@@ -36,6 +36,7 @@ from typing import Tuple, Union # noqa: F401
from types import SimpleNamespace
import numpy as np
+from pyproj import Geod
from mvgavg import mvgavg
from geoarray import GeoArray
from py_tools_ds.geo.coord_trafo import transform_any_prj
@@ -46,7 +47,6 @@ from ...model.images import EnMAPL1Product_SensorGeo, EnMAPL2Product_MapGeo
from ...model.metadata import EnMAP_Metadata_L2A_MapGeo
from ..spatial_transform import \
Geometry_Transformer, \
- Geometry_Transformer_3D, \
move_extent_to_coord_grid
__author__ = 'Daniel Scheffler'
@@ -72,6 +72,16 @@ class Orthorectifier(object):
raise RuntimeError('Expected a %s geolayer shape of %s or %s. Received %s.'
% (detector.detector_name, str(datashape), str(datashape[:2]), str(XY.shape)))
+ @staticmethod
+ def get_enmap_coordinate_grid_ll(lon: float, lat: float
+ ) -> (Tuple[float, float], Tuple[float, float]):
+ """Return EnMAP-like (30x30m) longitude/latitude pixel grid specs at the given position."""
+ geod = Geod(ellps="WGS84")
+ delta_lon = abs(lon - geod.fwd(lon, lat, az=90, dist=30)[0])
+ delta_lat = abs(lat - geod.fwd(lon, lat, az=0, dist=30)[1])
+
+ return (lon, lon + delta_lon), (lat, lat + delta_lat)
+
def run_transformation(self, enmap_ImageL1: EnMAPL1Product_SensorGeo) -> EnMAPL2Product_MapGeo:
self.validate_input(enmap_ImageL1)
@@ -85,36 +95,48 @@ class Orthorectifier(object):
lons_vnir, lats_vnir = enmap_ImageL1.vnir.detector_meta.lons, enmap_ImageL1.vnir.detector_meta.lats
lons_swir, lats_swir = enmap_ImageL1.swir.detector_meta.lons, enmap_ImageL1.swir.detector_meta.lats
+ if not enmap_ImageL1.vnir.detector_meta.geolayer_has_keystone and lons_vnir.ndim == 3:
+ lons_vnir, lats_vnir = lons_vnir[:, :, 0], lats_vnir[:, :, 0]
+ if not enmap_ImageL1.swir.detector_meta.geolayer_has_keystone and lons_swir.ndim == 3:
+ lons_swir, lats_swir = lons_swir[:, :, 0], lats_swir[:, :, 0]
# get target EPSG code and common extent
+ # (VNIR/SWIR overlap, i.e., INNER extent - non-overlapping parts are cleared later)
tgt_epsg = enmap_ImageL1.meta.vnir.epsg_ortho
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
- )
- if self.cfg.ortho_resampAlg == 'bilinear':
- # increase that if the resampling result contains gaps (default is 32 but this is quite slow)
- kw_init['neighbours'] = 8
-
- kw_trafo = dict(tgt_prj=tgt_epsg, tgt_extent=tgt_extent,
- tgt_coordgrid=((self.cfg.target_coord_grid['x'], self.cfg.target_coord_grid['y'])
- if self.cfg.target_coord_grid else None)
- )
- # transform VNIR and SWIR to map geometry
- GeoTransformer = Geometry_Transformer if lons_vnir.ndim == 2 else Geometry_Transformer_3D
+ # set up parameters for Geometry_Transformer initialization and execution of the transformation
+ kw_init = dict(
+ backend='gdal' if self.cfg.ortho_resampAlg != 'gauss' else 'pyresample',
+ resamp_alg=self.cfg.ortho_resampAlg,
+ nprocs=self.cfg.CPUs
+ )
+ kw_trafo = dict(
+ tgt_prj=tgt_epsg,
+ tgt_extent=tgt_extent,
+ tgt_coordgrid=((self.cfg.target_coord_grid['x'],
+ self.cfg.target_coord_grid['y'])
+ if self.cfg.target_coord_grid else
+ None),
+ src_nodata=enmap_ImageL1.vnir.data.nodata,
+ tgt_nodata=self.cfg.output_nodata_value
+ )
+ # make sure VNIR and SWIR are also transformed to the same lon/lat pixel grid
+ if self.cfg.target_projection_type == 'Geographic' and kw_trafo['tgt_coordgrid'] is None:
+ center_row, center_col = lons_vnir.shape[0] // 2, lons_vnir.shape[1] // 2
+ center_lon, center_lat = lons_vnir[center_row, center_col], lats_vnir[center_row, center_col]
+ kw_trafo['tgt_coordgrid'] = self.get_enmap_coordinate_grid_ll(center_lon, center_lat)
- # FIXME So far, the fill value is set to 0. Is this meaningful?
+ # transform VNIR and SWIR to map geometry
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=self.cfg.output_nodata_value, **kw_init)
+ GT_vnir = Geometry_Transformer(lons=lons_vnir, lats=lats_vnir, **kw_init)
vnir_mapgeo_gA = GeoArray(*GT_vnir.to_map_geometry(enmap_ImageL1.vnir.data, **kw_trafo),
nodata=self.cfg.output_nodata_value)
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=self.cfg.output_nodata_value, **kw_init)
+ GT_swir = Geometry_Transformer(lons=lons_swir, lats=lats_swir, **kw_init)
swir_mapgeo_gA = GeoArray(*GT_swir.to_map_geometry(enmap_ImageL1.swir.data, **kw_trafo),
nodata=self.cfg.output_nodata_value)
@@ -132,9 +154,9 @@ class Orthorectifier(object):
# TODO allow to set geolayer band to be used for warping of 2D arrays
# always use nearest neighbour resampling for masks and bitmasks with discrete values
- rsp_nearest_list = ['mask_landwater', 'mask_clouds', 'mask_cloudshadow', 'mask_haze', 'mask_snow',
- 'mask_cirrus', 'polymer_bitmask']
- kw_init_nearest = dict(resamp_alg='nearest', nprocs=self.cfg.CPUs)
+ rsp_nearest_list = ['mask_landwater', 'mask_clouds', 'mask_cloudshadow',
+ 'mask_haze', 'mask_snow', 'mask_cirrus', 'polymer_bitmask']
+ kw_init_nearest = dict(backend='gdal', resamp_alg='nearest', nprocs=self.cfg.CPUs)
# run the orthorectification
for attrName in ['mask_landwater', 'mask_clouds', 'mask_cloudshadow', 'mask_haze', 'mask_snow', 'mask_cirrus',
@@ -142,14 +164,18 @@ class Orthorectifier(object):
attr = getattr(enmap_ImageL1.vnir, attrName)
if attr is not None:
+ kw_init_attr = kw_init.copy() if attrName not in rsp_nearest_list else kw_init_nearest
+ kw_trafo_attr = kw_trafo.copy()
+ kw_trafo_attr['src_nodata'] = attr.nodata
+ kw_trafo_attr['tgt_nodata'] = attr.nodata
+
GT = Geometry_Transformer(
lons=lons_vnir if lons_vnir.ndim == 2 else lons_vnir[:, :, 0],
lats=lats_vnir if lats_vnir.ndim == 2 else lats_vnir[:, :, 0],
- fill_value=attr.nodata,
- **(kw_init if attrName not in rsp_nearest_list else kw_init_nearest))
+ **kw_init_attr)
enmap_ImageL1.logger.info("Orthorectifying '%s' attribute..." % attrName)
- attr_ortho = GeoArray(*GT.to_map_geometry(attr, **kw_trafo), nodata=attr.nodata)
+ attr_ortho = GeoArray(*GT.to_map_geometry(attr, **kw_trafo_attr), nodata=attr.nodata)
setattr(L2_obj, attrName, attr_ortho)
# TODO transform dead pixel map, quality test flags?
@@ -206,9 +232,13 @@ class Orthorectifier(object):
V_UL_UR_LL_LR_ll = [(V_lons[y, x], V_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, 0), (-1, -1)]]
S_UL_UR_LL_LR_ll = [(S_lons[y, x], S_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, 0), (-1, -1)]]
- # transform them to UTM
- V_UL_UR_LL_LR_prj = [transform_any_prj(4326, tgt_epsg, x, y) for x, y in V_UL_UR_LL_LR_ll]
- S_UL_UR_LL_LR_prj = [transform_any_prj(4326, tgt_epsg, x, y) for x, y in S_UL_UR_LL_LR_ll]
+ # transform them to tgt_epsg
+ if tgt_epsg != 4326:
+ V_UL_UR_LL_LR_prj = [transform_any_prj(4326, tgt_epsg, x, y) for x, y in V_UL_UR_LL_LR_ll]
+ S_UL_UR_LL_LR_prj = [transform_any_prj(4326, tgt_epsg, x, y) for x, y in S_UL_UR_LL_LR_ll]
+ else:
+ V_UL_UR_LL_LR_prj = V_UL_UR_LL_LR_ll
+ S_UL_UR_LL_LR_prj = S_UL_UR_LL_LR_ll
# separate X and Y
V_X_prj, V_Y_prj = zip(*V_UL_UR_LL_LR_prj)
@@ -216,7 +246,7 @@ class Orthorectifier(object):
# 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
- # (these pixels would be set to nodata later anyways, so we don't need to increase the extent for them)
+ # (these pixels would be set to nodata later anyway, so we don't need to increase the extent for them)
if V_lons.ndim == 3:
V_X_prj = (V_X_prj[0].max(), V_X_prj[1].min(), V_X_prj[2].max(), V_X_prj[3].min())
V_Y_prj = (V_Y_prj[0].min(), V_Y_prj[1].min(), V_Y_prj[2].max(), V_Y_prj[3].max())
diff --git a/enpt/processors/spatial_optimization/spatial_optimization.py b/enpt/processors/spatial_optimization/spatial_optimization.py
index b4832c59cc777a7ed257e4e0ce147adaa21b4872..fc9fc5d0029ebdd156d1cf391ab556da069c18e9 100644
--- a/enpt/processors/spatial_optimization/spatial_optimization.py
+++ b/enpt/processors/spatial_optimization/spatial_optimization.py
@@ -39,7 +39,7 @@ import os
from typing import Optional
import numpy as np
-from osgeo import gdal
+from osgeo import gdal # noqa
from arosics import COREG_LOCAL
from geoarray import GeoArray
@@ -49,7 +49,7 @@ from py_tools_ds.processing.progress_mon import ProgressBar
from ...options.config import EnPTConfig
from ...model.images.images_sensorgeo import EnMAPL1Product_SensorGeo
-from ..spatial_transform import Geometry_Transformer, Geometry_Transformer_3D
+from ..spatial_transform import Geometry_Transformer
class Spatial_Optimizer(object):
@@ -78,15 +78,17 @@ class Spatial_Optimizer(object):
% (self._EnMAP_bandIdx + 1))
GT = Geometry_Transformer(lons=self._EnMAP_Im.meta.vnir.lons[:, :, self._EnMAP_bandIdx],
lats=self._EnMAP_Im.meta.vnir.lats[:, :, self._EnMAP_bandIdx],
- fill_value=0,
- resamp_alg='gauss',
- radius_of_influence=30,
+ backend='gdal',
+ resamp_alg='bilinear',
nprocs=self.cfg.CPUs)
self._EnMAP_band = \
GeoArray(*GT.to_map_geometry(enmap_band_sensorgeo,
tgt_prj=self._get_optimal_utm_epsg(),
- tgt_coordgrid=((0, 15), (0, -15))),
+ tgt_coordgrid=((0, 15), (0, -15)),
+ src_nodata=self._EnMAP_Im.vnir.data.nodata,
+ tgt_nodata=0
+ ),
nodata=0)
return self._EnMAP_band
@@ -101,14 +103,17 @@ class Spatial_Optimizer(object):
self._EnMAP_Im.logger.info('Temporarily transforming EnMAP water mask to map geometry for co-registration.')
GT = Geometry_Transformer(lons=self._EnMAP_Im.meta.vnir.lons[:, :, self._EnMAP_bandIdx],
lats=self._EnMAP_Im.meta.vnir.lats[:, :, self._EnMAP_bandIdx],
- fill_value=0,
+ backend='gdal',
resamp_alg='nearest',
nprocs=self.cfg.CPUs)
self._EnMAP_mask = \
GeoArray(*GT.to_map_geometry(enmap_mask_sensorgeo,
tgt_prj=self._get_optimal_utm_epsg(),
- tgt_coordgrid=((0, 15), (0, -15))),
+ tgt_coordgrid=((0, 15), (0, -15)),
+ src_nodata=0, # 0=background
+ tgt_nodata=0
+ ),
nodata=0)
return self._EnMAP_mask
@@ -305,18 +310,18 @@ class Spatial_Optimizer(object):
enmap_ImageL1.logger.info('Transforming spatial optimization results back to sensor geometry.')
lons_band = self._EnMAP_Im.meta.vnir.lons[:, :, self._EnMAP_bandIdx]
lats_band = self._EnMAP_Im.meta.vnir.lats[:, :, self._EnMAP_bandIdx]
- GT = Geometry_Transformer_3D(lons=np.repeat(lons_band[:, :, np.newaxis], 2, axis=2),
- lats=np.repeat(lats_band[:, :, np.newaxis], 2, axis=2),
- fill_value=0,
- # resamp_alg='bilinear',
- resamp_alg='gauss',
- nprocs=self.cfg.CPUs)
+ GT = Geometry_Transformer(lons=lons_band,
+ lats=lats_band,
+ backend='gdal',
+ resamp_alg='bilinear',
+ nprocs=self.cfg.CPUs)
geolayer_sensorgeo = \
GT.to_sensor_geometry(GeoArray(np.dstack([xgrid_map_coreg,
ygrid_map_coreg]),
geotransform=self._EnMAP_band.gt,
- projection=self._EnMAP_band.prj))
+ projection=self._EnMAP_band.prj),
+ tgt_nodata=0)
enmap_ImageL1.logger.info('Applying results of spatial optimization to geolayer.')
lons_coreg, lats_coreg = transform_coordArray(prj_src=self._ref_band_prep.prj,
diff --git a/enpt/processors/spatial_transform/spatial_transform.py b/enpt/processors/spatial_transform/spatial_transform.py
index 0109c15334494186c441ad4ee4321448cbfde895..17dcbf0581b2e446d37d92032c912075c663bd99 100644
--- a/enpt/processors/spatial_transform/spatial_transform.py
+++ b/enpt/processors/spatial_transform/spatial_transform.py
@@ -28,7 +28,6 @@
# with this program. If not, see <https://www.gnu.org/licenses/>.
"""EnPT module 'spatial transform', containing everything related to spatial transformations."""
-import warnings
from typing import Union, Tuple, List, Optional, Sequence # noqa: F401
from multiprocessing import Pool, cpu_count
from collections import OrderedDict
@@ -40,8 +39,8 @@ from natsort import natsorted
import numpy_indexed as npi
from pyproj import CRS
-from sensormapgeo import SensorMapGeometryTransformer, SensorMapGeometryTransformer3D
-from sensormapgeo.transformer_2d import AreaDefinition
+from sensormapgeo import Transformer
+from sensormapgeo.pyresample_backend.transformer_2d import AreaDefinition
from py_tools_ds.geo.projection import prj_equal
from py_tools_ds.geo.coord_grid import find_nearest
from py_tools_ds.geo.coord_trafo import transform_any_prj, transform_coordArray
@@ -51,86 +50,38 @@ from ...options.config import enmap_coordinate_grid_utm
__author__ = 'Daniel Scheffler'
-class Geometry_Transformer(SensorMapGeometryTransformer):
+class Geometry_Transformer(Transformer):
# use Sentinel-2 grid (30m grid with origin at 0/0)
# EnMAP geolayer contains pixel center coordinate
def to_sensor_geometry(self,
path_or_geoarray_mapgeo: Union[str, GeoArray],
+ src_gt: Tuple[float, float, float, float, float, float] = None,
src_prj: Union[str, int] = None,
- src_extent: Tuple[float, float, float, float] = None):
+ src_nodata: int = None,
+ tgt_nodata: int = None
+ ) -> np.ndarray:
data_mapgeo = GeoArray(path_or_geoarray_mapgeo)
if not data_mapgeo.is_map_geo:
raise RuntimeError('The dataset to be transformed into sensor geometry already represents sensor geometry.')
- with warnings.catch_warnings():
- # FIXME remove that after removing pyresample pinning
- warnings.filterwarnings('ignore', category=DeprecationWarning,
- message='.*is a deprecated alias for the builtin.*')
-
- return super().to_sensor_geometry(
- data_mapgeo[:],
- src_prj=src_prj or data_mapgeo.prj,
- src_extent=src_extent or list(np.array(data_mapgeo.box.boundsMap)[[0, 2, 1, 3]]))
-
- def to_map_geometry(self,
- 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_coordgrid: Tuple[Tuple, Tuple] = None,
- area_definition: AreaDefinition = None):
- data_sensorgeo = GeoArray(path_or_geoarray_sensorgeo)
-
- if data_sensorgeo.is_map_geo and not data_sensorgeo.is_rotated:
- raise RuntimeError('The dataset to be transformed into map geometry already represents map geometry.')
-
- # run transformation (output extent/area definition etc. is internally computed from the geolayers if not given)
- with warnings.catch_warnings():
- # FIXME remove that after removing pyresample pinning
- warnings.filterwarnings('ignore', category=DeprecationWarning,
- message='.*is a deprecated alias for the builtin.*')
-
- out_data, out_gt, out_prj = \
- super(Geometry_Transformer, self).to_map_geometry(data_sensorgeo[:],
- tgt_prj=tgt_prj,
- tgt_extent=tgt_extent,
- tgt_res=tgt_res,
- tgt_coordgrid=tgt_coordgrid,
- area_definition=self.area_definition)
-
- return out_data, out_gt, out_prj
-
-
-class Geometry_Transformer_3D(SensorMapGeometryTransformer3D):
- # use Sentinel-2 grid (30m grid with origin at 0/0)
- # EnMAP geolayer contains pixel center coordinate
- def to_sensor_geometry(self,
- path_or_geoarray_mapgeo: Union[str, GeoArray],
- src_prj: Union[str, int] = None,
- src_extent: Tuple[float, float, float, float] = None):
- data_mapgeo = GeoArray(path_or_geoarray_mapgeo)
-
- if not data_mapgeo.is_map_geo:
- raise RuntimeError('The dataset to be transformed into sensor geometry already represents sensor geometry.')
-
- with warnings.catch_warnings():
- # FIXME remove that after removing pyresample pinning
- warnings.filterwarnings('ignore', category=DeprecationWarning,
- message='.*is a deprecated alias for the builtin.*')
-
- return super().to_sensor_geometry(
- data_mapgeo[:],
- src_prj=src_prj or data_mapgeo.prj,
- src_extent=src_extent or list(np.array(data_mapgeo.box.boundsMap)[[0, 2, 1, 3]]))
+ return super().to_sensor_geometry(
+ data_mapgeo[:],
+ src_gt=src_gt or data_mapgeo.gt,
+ src_prj=src_prj or data_mapgeo.prj,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata
+ )
def to_map_geometry(self,
path_or_geoarray_sensorgeo: Union[str, GeoArray, np.ndarray],
- tgt_prj: Union[str, int] = None,
+ tgt_prj: Union[str, int],
tgt_extent: Tuple[float, float, float, float] = None,
- tgt_res: Tuple[float, float] = None,
+ tgt_res: Tuple = None,
tgt_coordgrid: Tuple[Tuple, Tuple] = None,
+ src_nodata: int = None,
+ tgt_nodata: int = None,
area_definition: AreaDefinition = None
) -> Tuple[np.ndarray, tuple, str]:
data_sensorgeo = GeoArray(path_or_geoarray_sensorgeo)
@@ -139,18 +90,17 @@ class Geometry_Transformer_3D(SensorMapGeometryTransformer3D):
raise RuntimeError('The dataset to be transformed into map geometry already represents map geometry.')
# run transformation (output extent/area definition etc. is internally computed from the geolayers if not given)
- with warnings.catch_warnings():
- # FIXME remove that after removing pyresample pinning
- warnings.filterwarnings('ignore', category=DeprecationWarning,
- message='.*is a deprecated alias for the builtin.*')
-
- out_data, out_gt, out_prj = \
- super(Geometry_Transformer_3D, self).to_map_geometry(data_sensorgeo[:],
- tgt_prj=tgt_prj,
- tgt_extent=tgt_extent,
- tgt_res=tgt_res,
- tgt_coordgrid=tgt_coordgrid,
- area_definition=area_definition)
+ out_data, out_gt, out_prj = \
+ super().to_map_geometry(
+ data_sensorgeo[:],
+ tgt_prj=tgt_prj,
+ tgt_extent=tgt_extent,
+ tgt_res=tgt_res,
+ tgt_coordgrid=tgt_coordgrid,
+ src_nodata=src_nodata,
+ tgt_nodata=tgt_nodata,
+ area_definition=area_definition
+ )
return out_data, out_gt, out_prj
@@ -168,7 +118,9 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
res_vnir: Tuple[float, float],
res_swir: Tuple[float, float],
resamp_alg: str = 'nearest',
- **gt_opts) -> None:
+ src_nodata: int = None,
+ tgt_nodata: int = None,
+ nprocs: int = cpu_count()) -> None:
"""Get an instance of VNIR_SWIR_SensorGeometryTransformer.
:param lons_vnir: VNIR longitude array corresponding to the sensor geometry arrays passed later (2D or 3D)
@@ -179,9 +131,11 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
:param prj_swir: projection of the SWIR if it would be transformed to map geometry (WKT string or EPSG code)
:param res_vnir: pixel dimensions of the VNIR if it would be transformed to map geometry (X, Y)
:param res_swir: pixel dimensions of the SWIR if it would be transformed to map geometry (X, Y)
- :param resamp_alg: resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
- :param gt_opts: additional options to be passed to the Geometric_Transformer class,
- e.g., 'fill_value', 'radius_of_influence', 'nprocs'...
+ :param resamp_alg: resampling algorithm ('nearest', 'near', 'bilinear', 'cubic', 'cubic_spline',
+ 'lanczos', 'average', 'mode', 'max', 'min', 'med', 'q1', 'q3')
+ :param src_nodata: nodata value of the input array
+ :param tgt_nodata: pixel value to be used for undetermined pixels in the output
+ :param nprocs: number of CPU cores to use
"""
self.vnir_meta = dict(
lons=lons_vnir,
@@ -196,7 +150,9 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
res=res_swir
)
self.resamp_alg = resamp_alg
- self.gt_opts = gt_opts
+ self.src_nodata = src_nodata
+ self.tgt_nodata = tgt_nodata
+ self.cpus = nprocs
def transform_sensorgeo_VNIR_to_SWIR(self, data_vnirsensorgeo: np.ndarray) -> np.ndarray:
"""Transform any array in VNIR sensor geometry to SWIR sensor geometry to remove geometric shifts.
@@ -216,12 +172,15 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
def _transform_sensorgeo(self,
data2transform: np.ndarray,
- inputgeo: str) -> np.ndarray:
+ inputgeo: str,
+ band_outputgeo: int = 0
+ ) -> np.ndarray:
"""Transform the input array between VNIR and SWIR sensor geometry.
:param data2transform: input array to be transformed
:param inputgeo: 'vnir' if data2transform is given in VNIR sensor geometry
'swir' if data2transform is given in SWIR sensor geometry
+ :param band_outputgeo: band index of the band from the sensor geometry should be used to generate the output
"""
# TODO: Avoid the resampling here, maybe by replacing the lon/lat arrays by image coordinates for the source
# geometry and by image coordinate differences for the target geometry. Maybe also the proj string for
@@ -232,35 +191,40 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
src, tgt = (self.vnir_meta, self.swir_meta) if inputgeo == 'vnir' else (self.swir_meta, self.vnir_meta)
- # get 2D or 3D GeoTransformer
- if data2transform.ndim == 3 and src['lons'].ndim == 3 and src['lons'].shape[2] == data2transform.shape[2]:
- GT_src = Geometry_Transformer_3D(lons=src['lons'],
- lats=src['lats'],
- resamp_alg=self.resamp_alg,
- **self.gt_opts)
- GT_tgt = Geometry_Transformer_3D(lons=tgt['lons'],
- lats=tgt['lats'],
- resamp_alg=self.resamp_alg,
- **self.gt_opts)
+ # get GeoTransformer
+ # NOTE: We can only use a 3D geolayer to temporarily transform to map geometry but not back to sensor geometry
+ # of the other detector. For this, we need to select the band of which the sensor geometry should be used.
+ def ensure_2d(coord_array):
+ return coord_array if coord_array.ndim == 2 else coord_array[:, :, band_outputgeo]
+ if data2transform.ndim == 3 and src['lons'].ndim == 3 and src['lons'].shape[2] == data2transform.shape[2]:
+ src_lons, src_lats = src['lons'], src['lats']
+ tgt_lons, tgt_lats = ensure_2d(tgt['lons']), ensure_2d(tgt['lats'])
else:
- def ensure_2D(coord_array):
- return coord_array if coord_array.ndim == 2 else coord_array[:, :, 0]
-
- GT_src = Geometry_Transformer(lons=ensure_2D(src['lons']),
- lats=ensure_2D(src['lats']),
- resamp_alg=self.resamp_alg,
- **self.gt_opts)
- GT_tgt = Geometry_Transformer(lons=ensure_2D(tgt['lons']),
- lats=ensure_2D(tgt['lats']),
- resamp_alg=self.resamp_alg,
- **self.gt_opts)
+ src_lons, src_lats = ensure_2d(src['lons']), ensure_2d(src['lats'])
+ tgt_lons, tgt_lats = ensure_2d(tgt['lons']), ensure_2d(tgt['lats'])
+
+ GT_src = Geometry_Transformer(src_lons, src_lats, backend='gdal', resamp_alg=self.resamp_alg, nprocs=self.cpus)
+ GT_tgt = Geometry_Transformer(tgt_lons, tgt_lats, backend='gdal', resamp_alg=self.resamp_alg, nprocs=self.cpus)
# temporarily transform the input sensor geometry array to map geometry
- gA_mapgeo = GeoArray(*GT_src.to_map_geometry(data2transform, tgt_prj=src['prj']))
+ gA_mapgeo = GeoArray(
+ *GT_src.to_map_geometry(
+ data2transform,
+ tgt_prj=src['prj'],
+ src_nodata=self.src_nodata,
+ tgt_nodata=self.tgt_nodata
+ ),
+ nodata=self.tgt_nodata
+ )
- # generate the target sensor geometry array (target lons/lats define the target swath definition)
- tgt_data_sensorgeo = GT_tgt.to_sensor_geometry(gA_mapgeo)
+ # generate the target sensor geometry array (target lons/lats define the target swath definition
+ tgt_data_sensorgeo =\
+ GT_tgt.to_sensor_geometry(
+ gA_mapgeo,
+ src_nodata=gA_mapgeo._nodata, # noqa
+ tgt_nodata=self.tgt_nodata
+ )
return tgt_data_sensorgeo
@@ -686,7 +650,7 @@ def compute_mapCoords_within_sensorGeoDims(sensorgeoCoords_YX: List[Tuple[float,
enmapIm_cornerCoords: Tuple[Tuple[float, float]],
enmapIm_dims_sensorgeo: Tuple[int, int],
) -> List[Tuple[float, float]]:
- """Compute map coordinates for a given image coordinate pair of an EnMAP image in sensor geometry.
+ """Compute map coordinates for a given image coordinate-pair of an EnMAP image in sensor geometry.
:param sensorgeoCoords_YX: list of requested sensor geometry positions [(row, column), (row, column), ...]
:param rpc_coeffs: RPC coefficients describing the relation between sensor and map geometry
diff --git a/requirements.txt b/requirements.txt
index de0b5bf965b00e7335d50b61cd064b9f66447410..6174902902616447860837294c197743c7b19a52 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -16,7 +16,7 @@ pyproj>=3.4.0
py_tools_ds>=0.14.25
scikit-image
scipy
-sensormapgeo>=0.4.0
+sensormapgeo>=1.0.0
sicor>=0.19.1
tqdm
utm
diff --git a/setup.py b/setup.py
index 00904bff3e96111b2611dea89fe727855e6ef34b..ce02015556f49eb4c09fe865645d88881cace97e 100644
--- a/setup.py
+++ b/setup.py
@@ -59,7 +59,7 @@ req = [
'py_tools_ds>=0.14.25',
'scikit-image',
'scipy',
- 'sensormapgeo>=0.4.0',
+ 'sensormapgeo>=1.0.0',
'sicor>=0.19.1',
'tqdm',
'utm',
diff --git a/tests/gitlab_CI_docker/context/environment_enpt.yml b/tests/gitlab_CI_docker/context/environment_enpt.yml
index 186183d438cdd18fc644b33bf036611369438347..05c88c7b7ae078fbf5dfe375d3e91d9224e121a5 100644
--- a/tests/gitlab_CI_docker/context/environment_enpt.yml
+++ b/tests/gitlab_CI_docker/context/environment_enpt.yml
@@ -24,7 +24,7 @@ dependencies:
- py-tools-ds>=0.14.25
- scikit-image
- scipy
- - sensormapgeo>=0.4.0
+ - sensormapgeo>=1.0.0
- sicor>=0.19.1
- tqdm
- utm