Commit 31ac77d7 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed issue where SensorMapGeometryTransformer raised an exception when trying...


Fixed issue where SensorMapGeometryTransformer raised an exception when trying to resample a 3D input array.
Fixed SensorMapGeometryTransformer.to_sensor_geometry() not working for resamp_alg='bilinear'.
The test_reproject module now tests all resampling algorithms instead of only one per test.

Updated version info.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 4c725e9c
Pipeline #5069 passed with stage
in 1 minute and 16 seconds
...@@ -585,6 +585,11 @@ class SensorMapGeometryTransformer(object): ...@@ -585,6 +585,11 @@ class SensorMapGeometryTransformer(object):
self.lats = lats self.lats = lats
self.lons = lons self.lons = lons
self.swath_definition = SwathDefinition(lons=lons, lats=lats) self.swath_definition = SwathDefinition(lons=lons, lats=lats)
# use a projection string for local coordinates (https://gis.stackexchange.com/a/300407)
# -> this is needed for bilinear resampling
self.swath_definition.proj_str = '+proj=omerc +lat_0=51.6959777875 +lonc=7.0923165808 +alpha=-20.145 ' \
'+gamma=0 +k=1 +x_0=50692.579 +y_0=81723.458 +ellps=WGS84 ' \
'+towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
self.area_extent_ll = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)] self.area_extent_ll = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)]
self.area_definition = None # type: AreaDefinition self.area_definition = None # type: AreaDefinition
...@@ -733,6 +738,15 @@ class SensorMapGeometryTransformer(object): ...@@ -733,6 +738,15 @@ class SensorMapGeometryTransformer(object):
elif self.resamp_alg == 'gauss': elif self.resamp_alg == 'gauss':
opts = {k: v for k, v in self.opts.items()} opts = {k: v for k, v in self.opts.items()}
# ensure that sigmas are provided as list if data is 3-dimensional
if data.ndim != 2:
if not isinstance(opts['sigmas'], list):
opts['sigmas'] = [opts['sigmas']] * data.ndim
if not len(opts['sigmas']) == data.ndim:
raise ValueError("The 'sigmas' parameter must have the same number of values like data.ndim."
"n_sigmas: %d; data.ndim: %d" % (len(opts['sigmas']), data.ndim))
result = resample_gauss(source_geo_def, data, target_geo_def, **opts) result = resample_gauss(source_geo_def, data, target_geo_def, **opts)
elif self.resamp_alg == 'custom': elif self.resamp_alg == 'custom':
......
...@@ -19,5 +19,5 @@ ...@@ -19,5 +19,5 @@
# You should have received a copy of the GNU Lesser General Public License along # You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>. # with this program. If not, see <http://www.gnu.org/licenses/>.
__version__ = '0.14.22' __version__ = '0.14.23'
__versionalias__ = '20191016_01' __versionalias__ = '20191017_01'
...@@ -41,6 +41,7 @@ from py_tools_ds.geo.raster.reproject import SensorMapGeometryTransformer, Senso ...@@ -41,6 +41,7 @@ from py_tools_ds.geo.raster.reproject import SensorMapGeometryTransformer, Senso
tests_path = os.path.abspath(os.path.join(__path__[0], "..", "tests")) tests_path = os.path.abspath(os.path.join(__path__[0], "..", "tests"))
rsp_algs = ['nearest', 'bilinear', 'gauss']
class Test_SensorMapGeometryTransformer(TestCase): class Test_SensorMapGeometryTransformer(TestCase):
...@@ -65,61 +66,66 @@ class Test_SensorMapGeometryTransformer(TestCase): ...@@ -65,61 +66,66 @@ class Test_SensorMapGeometryTransformer(TestCase):
5267203.56579] # UR_y 5267203.56579] # UR_y
def test_to_sensor_geometry(self): def test_to_sensor_geometry(self):
SMGT = SensorMapGeometryTransformer(lons=self.lons, for rsp_alg in rsp_algs:
lats=self.lats, SMGT = SensorMapGeometryTransformer(lons=self.lons,
resamp_alg='nearest') lats=self.lats,
dem_sensors_geo = SMGT.to_sensor_geometry(self.dem_map_geo, resamp_alg=rsp_alg,
src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm) radius_of_influence=30 if rsp_alg != 'bilinear' else 45)
self.assertIsInstance(dem_sensors_geo, np.ndarray) dem_sensors_geo = SMGT.to_sensor_geometry(self.dem_map_geo,
self.assertFalse(np.array_equal(np.unique(dem_sensors_geo), np.array([0]))) src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
self.assertEquals(dem_sensors_geo.shape, (150, 1000)) self.assertIsInstance(dem_sensors_geo, np.ndarray)
self.assertFalse(np.array_equal(np.unique(dem_sensors_geo), np.array([0])))
self.assertEquals(dem_sensors_geo.shape, (150, 1000))
def test_to_sensor_geometry_3DInput(self): def test_to_sensor_geometry_3DInput(self):
SMGT = SensorMapGeometryTransformer(lons=self.lons, for rsp_alg in rsp_algs:
lats=self.lats, SMGT = SensorMapGeometryTransformer(lons=self.lons,
resamp_alg='nearest') lats=self.lats,
dem_sensors_geo = SMGT.to_sensor_geometry(np.dstack([self.dem_map_geo] * 2), resamp_alg=rsp_alg)
src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm) dem_sensors_geo = SMGT.to_sensor_geometry(np.dstack([self.dem_map_geo] * 2),
self.assertIsInstance(dem_sensors_geo, np.ndarray) src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
self.assertFalse(np.array_equal(np.unique(dem_sensors_geo), np.array([0]))) self.assertIsInstance(dem_sensors_geo, np.ndarray)
self.assertEquals(dem_sensors_geo.shape, (150, 1000, 2)) self.assertFalse(np.array_equal(np.unique(dem_sensors_geo), np.array([0])))
self.assertTrue(np.array_equal(dem_sensors_geo[:, :, 0], dem_sensors_geo[:, :, 1])) self.assertEquals(dem_sensors_geo.shape, (150, 1000, 2))
self.assertTrue(np.array_equal(dem_sensors_geo[:, :, 0], dem_sensors_geo[:, :, 1]))
def test_to_map_geometry_lonlat(self): def test_to_map_geometry_lonlat(self):
SMGT = SensorMapGeometryTransformer(lons=self.lons, for rsp_alg in rsp_algs:
lats=self.lats, SMGT = SensorMapGeometryTransformer(lons=self.lons,
resamp_alg='nearest') lats=self.lats,
resamp_alg=rsp_alg)
# to Lon/Lat
dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=4326) # to Lon/Lat
self.assertIsInstance(dem_map_geo, np.ndarray) dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=4326)
self.assertEquals(dem_map_geo.shape, (SMGT.area_definition.height, self.assertIsInstance(dem_map_geo, np.ndarray)
SMGT.area_definition.width)) self.assertEquals(dem_map_geo.shape, (SMGT.area_definition.height,
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt, SMGT.area_definition.width))
cols=dem_map_geo.shape[1], xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
rows=dem_map_geo.shape[0])) cols=dem_map_geo.shape[1],
self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]), rows=dem_map_geo.shape[0]))
np.array(self.expected_dem_area_extent_lonlat))) self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0]))) np.array(self.expected_dem_area_extent_lonlat)))
self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
with self.assertRaises(ValueError):
SMGT.to_map_geometry(self.dem_sensor_geo[:10, :10], tgt_prj=4326) # must have the shape of lons/lats with self.assertRaises(ValueError):
SMGT.to_map_geometry(self.dem_sensor_geo[:10, :10], tgt_prj=4326) # must have the shape of lons/lats
def test_to_map_geometry_utm(self): def test_to_map_geometry_utm(self):
SMGT = SensorMapGeometryTransformer(lons=self.lons, for rsp_alg in rsp_algs:
lats=self.lats, SMGT = SensorMapGeometryTransformer(lons=self.lons,
resamp_alg='nearest') lats=self.lats,
resamp_alg=rsp_alg)
# to UTM32
dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30)) # to UTM32
self.assertIsInstance(dem_map_geo, np.ndarray) dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30))
self.assertEquals(dem_map_geo.shape, (365, 975)) self.assertIsInstance(dem_map_geo, np.ndarray)
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt, self.assertEquals(dem_map_geo.shape, (365, 975))
cols=dem_map_geo.shape[1], xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
rows=dem_map_geo.shape[0])) cols=dem_map_geo.shape[1],
self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]), rows=dem_map_geo.shape[0]))
np.array(self.expected_dem_area_extent_utm))) self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0]))) np.array(self.expected_dem_area_extent_utm)))
self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
class Test_SensorMapGeometryTransformer3D(TestCase): class Test_SensorMapGeometryTransformer3D(TestCase):
...@@ -150,30 +156,32 @@ class Test_SensorMapGeometryTransformer3D(TestCase): ...@@ -150,30 +156,32 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
5267203.56579] # UR_y 5267203.56579] # UR_y
def test_to_map_geometry_lonlat_3D_geolayer(self): def test_to_map_geometry_lonlat_3D_geolayer(self):
SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D, for rsp_alg in rsp_algs:
lats=self.lats_3D, SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
# resamp_alg='nearest', lats=self.lats_3D,
resamp_alg='bilinear', # resamp_alg='nearest',
) resamp_alg=rsp_alg,
)
# to Lon/Lat
data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326) # to Lon/Lat
self.assertIsInstance(data_mapgeo_3D, np.ndarray) data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)
# only validate number of bands (height and width are validated in 2D version self.assertIsInstance(data_mapgeo_3D, np.ndarray)
# fixed numbers may fail here due to float uncertainty errors # only validate number of bands (height and width are validated in 2D version
self.assertEquals(data_mapgeo_3D.shape[2], 2) # fixed numbers may fail here due to float uncertainty errors
xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt, self.assertEquals(data_mapgeo_3D.shape[2], 2)
cols=data_mapgeo_3D.shape[1], xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
rows=data_mapgeo_3D.shape[0])) cols=data_mapgeo_3D.shape[1],
self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]), rows=data_mapgeo_3D.shape[0]))
np.array(self.expected_dem_area_extent_lonlat))) self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
np.array(self.expected_dem_area_extent_lonlat)))
def test_to_sensor_geometry(self): def test_to_sensor_geometry(self):
SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D, for rsp_alg in rsp_algs:
lats=self.lats_3D, SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
resamp_alg='nearest', lats=self.lats_3D,
) resamp_alg=rsp_alg,
dem_sensors_geo = SMGT.to_sensor_geometry(self.data_map_geo_3D, )
src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm) dem_sensors_geo = SMGT.to_sensor_geometry(self.data_map_geo_3D,
self.assertIsInstance(dem_sensors_geo, np.ndarray) src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
self.assertEquals(dem_sensors_geo.shape, (150, 1000, 2)) self.assertIsInstance(dem_sensors_geo, np.ndarray)
self.assertEquals(dem_sensors_geo.shape, (150, 1000, 2))
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