Merge branch 'bugfix/fix_bilinear_sensorgeo_gauss_3D_trafo' into 'master'

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

See merge request !14

py_tools_ds/geo/raster/reproject.py

@@ -585,6 +585,11 @@ class SensorMapGeometryTransformer(object):
         self.lats = lats
         self.lons = lons
         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_definition = None  # type: AreaDefinition
 
@@ -733,6 +738,15 @@ class SensorMapGeometryTransformer(object):
 
         elif self.resamp_alg == 'gauss':
             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)
 
         elif self.resamp_alg == 'custom':

py_tools_ds/version.py

@@ -19,5 +19,5 @@
 # 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/>.
 
-__version__ = '0.14.22'
-__versionalias__ = '20191016_01'
+__version__ = '0.14.23'
+__versionalias__ = '20191017_01'

tests/test_geo/test_raster/test_reproject.py

@@ -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"))
+rsp_algs = ['nearest', 'bilinear', 'gauss']
 
 
 class Test_SensorMapGeometryTransformer(TestCase):
@@ -65,61 +66,66 @@ class Test_SensorMapGeometryTransformer(TestCase):
                                              5267203.56579]  # UR_y
 
     def test_to_sensor_geometry(self):
-        SMGT = SensorMapGeometryTransformer(lons=self.lons,
-                                            lats=self.lats,
-                                            resamp_alg='nearest')
-        dem_sensors_geo = SMGT.to_sensor_geometry(self.dem_map_geo,
-                                                  src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
-        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))
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer(lons=self.lons,
+                                                lats=self.lats,
+                                                resamp_alg=rsp_alg,
+                                                radius_of_influence=30 if rsp_alg != 'bilinear' else 45)
+            dem_sensors_geo = SMGT.to_sensor_geometry(self.dem_map_geo,
+                                                      src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
+            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):
-        SMGT = SensorMapGeometryTransformer(lons=self.lons,
-                                            lats=self.lats,
-                                            resamp_alg='nearest')
-        dem_sensors_geo = SMGT.to_sensor_geometry(np.dstack([self.dem_map_geo] * 2),
-                                                  src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
-        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, 2))
-        self.assertTrue(np.array_equal(dem_sensors_geo[:, :, 0], dem_sensors_geo[:, :, 1]))
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer(lons=self.lons,
+                                                lats=self.lats,
+                                                resamp_alg=rsp_alg)
+            dem_sensors_geo = SMGT.to_sensor_geometry(np.dstack([self.dem_map_geo] * 2),
+                                                      src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
+            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, 2))
+            self.assertTrue(np.array_equal(dem_sensors_geo[:, :, 0], dem_sensors_geo[:, :, 1]))
 
     def test_to_map_geometry_lonlat(self):
-        SMGT = SensorMapGeometryTransformer(lons=self.lons,
-                                            lats=self.lats,
-                                            resamp_alg='nearest')
-
-        # to Lon/Lat
-        dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=4326)
-        self.assertIsInstance(dem_map_geo, np.ndarray)
-        self.assertEquals(dem_map_geo.shape, (SMGT.area_definition.height,
-                                              SMGT.area_definition.width))
-        xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
-                                                                               cols=dem_map_geo.shape[1],
-                                                                               rows=dem_map_geo.shape[0]))
-        self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
-                                                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
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer(lons=self.lons,
+                                                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)
+            self.assertIsInstance(dem_map_geo, np.ndarray)
+            self.assertEquals(dem_map_geo.shape, (SMGT.area_definition.height,
+                                                  SMGT.area_definition.width))
+            xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
+                                                                                   cols=dem_map_geo.shape[1],
+                                                                                   rows=dem_map_geo.shape[0]))
+            self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
+                                                    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
 
     def test_to_map_geometry_utm(self):
-        SMGT = SensorMapGeometryTransformer(lons=self.lons,
-                                            lats=self.lats,
-                                            resamp_alg='nearest')
-
-        # to UTM32
-        dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30))
-        self.assertIsInstance(dem_map_geo, np.ndarray)
-        self.assertEquals(dem_map_geo.shape, (365, 975))
-        xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
-                                                                               cols=dem_map_geo.shape[1],
-                                                                               rows=dem_map_geo.shape[0]))
-        self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
-                                                np.array(self.expected_dem_area_extent_utm)))
-        self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer(lons=self.lons,
+                                                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))
+            self.assertIsInstance(dem_map_geo, np.ndarray)
+            self.assertEquals(dem_map_geo.shape, (365, 975))
+            xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
+                                                                                   cols=dem_map_geo.shape[1],
+                                                                                   rows=dem_map_geo.shape[0]))
+            self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
+                                                    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):
@@ -150,30 +156,32 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
                                              5267203.56579]  # UR_y
 
     def test_to_map_geometry_lonlat_3D_geolayer(self):
-        SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
-                                              lats=self.lats_3D,
-                                              # resamp_alg='nearest',
-                                              resamp_alg='bilinear',
-                                              )
-
-        # to Lon/Lat
-        data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)
-        self.assertIsInstance(data_mapgeo_3D, np.ndarray)
-        # only validate number of bands (height and width are validated in 2D version
-        #   fixed numbers may fail here due to float uncertainty errors
-        self.assertEquals(data_mapgeo_3D.shape[2], 2)
-        xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
-                                                                               cols=data_mapgeo_3D.shape[1],
-                                                                               rows=data_mapgeo_3D.shape[0]))
-        self.assertTrue(False not in np.isclose(np.array([xmin, ymin, xmax, ymax]),
-                                                np.array(self.expected_dem_area_extent_lonlat)))
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
+                                                  lats=self.lats_3D,
+                                                  # 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)
+            self.assertIsInstance(data_mapgeo_3D, np.ndarray)
+            # only validate number of bands (height and width are validated in 2D version
+            #   fixed numbers may fail here due to float uncertainty errors
+            self.assertEquals(data_mapgeo_3D.shape[2], 2)
+            xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=dem_gt,
+                                                                                   cols=data_mapgeo_3D.shape[1],
+                                                                                   rows=data_mapgeo_3D.shape[0]))
+            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):
-        SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
-                                              lats=self.lats_3D,
-                                              resamp_alg='nearest',
-                                              )
-        dem_sensors_geo = SMGT.to_sensor_geometry(self.data_map_geo_3D,
-                                                  src_prj=32632, src_extent=self.dem_area_extent_coarse_subset_utm)
-        self.assertIsInstance(dem_sensors_geo, np.ndarray)
-        self.assertEquals(dem_sensors_geo.shape, (150, 1000,  2))
+        for rsp_alg in rsp_algs:
+            SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
+                                                  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)
+            self.assertIsInstance(dem_sensors_geo, np.ndarray)
+            self.assertEquals(dem_sensors_geo.shape, (150, 1000,  2))