Merge branch 'enhancement/speed_up_3Dconv' into 'master'

Enhancement/speed up 3Dconv

See merge request !2

HISTORY.rst

@@ -2,14 +2,40 @@
 History
 =======
 
-0.3.0 (coming soon)
--------------------
+0.x.x (2020-05-18)
+------------------
+
+* Converted all type hints to Python 3.6 style. Dropped Python 3.5 support. Fixed code duplicate.
+* Split sensormapgeo module into transformer_2d and transformer_3d.
+* SensorMapGeometryTransformer.compute_areadefinition_sensor2map() now directly uses pyresample instead of GDAL if the
+  target resolution is given.
+* SensorMapGeometryTransformer3D.to_map_geometry() now computes a common area definition only ONCE which saves
+  computation time and increases stability.
+* The computation of the common extent in 3D geolayers now works properly if target projection is not set to LonLat.
+* Added paramter tgt_coordgrid to to_map_geometry methods to automatically move the output extent to a given coordinate
+  grid.
+* compute_areadefinition_sensor2map() now also adds 1 pixel around the output extent in the pyresample version just
+  like in the GDAL version.
+* Added some input validation.
+
+
+0.2.2 (2020-03-10)
+------------------
+
+* Fix for always returning 0.1.0 when calling sensormapgeo.__version__.
+
+
+0.2.1 (2020-03-10)
+------------------
 
 * Fix for always returning returning float64 output data type in case of bilinear resampling.
 * Added output data type verification to tests.
 * Fix for an exception if the output of get_proj4info() contains trailing white spaces
   (fixed by an update of py_tools_ds).
-* Fix for always returning 0.1.0 when calling sensormapgeo.__version__.
+* Improved tests.
+* Set channel priority to strict.
+* Force libgdal to be installed from conda-forge.
+* Fixed broken documentation link
 
 
 0.2.0 (2020-01-06)

sensormapgeo/__init__.py

@@ -24,7 +24,8 @@
 
 """Top-level package for sensormapgeo."""
 
-from .sensormapgeo import SensorMapGeometryTransformer, SensorMapGeometryTransformer3D
+from .transformer_2d import SensorMapGeometryTransformer
+from .transformer_3d import SensorMapGeometryTransformer3D
 from .version import __version__, __versionalias__   # noqa (E402 + F401)
 
 __all__ = [

sensormapgeo/transformer_3d.py

+# -*- coding: utf-8 -*-
+
+# sensormapgeo, Transform remote sensing images between sensor and map geometry.
+#
+# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, danschef@gfz-potsdam.de)
+#
+# This software was developed within the context of the EnMAP project supported
+# by the DLR Space Administration with funds of the German Federal Ministry of
+# Economic Affairs and Energy (on the basis of a decision by the German Bundestag:
+# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG.
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# 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/>.
+
+"""Module to transform 3D arrays between sensor and map geometry (using band-wise Lon/Lat arrays)."""
+
+from typing import Union, Tuple
+import multiprocessing
+
+import numpy as np
+from py_tools_ds.geo.projection import WKT2EPSG, EPSG2WKT, proj4_to_WKT
+from py_tools_ds.geo.coord_trafo import get_proj4info, transform_any_prj
+
+from .transformer_2d import \
+    SensorMapGeometryTransformer, _corner_coords_lonlat_to_extent, \
+    _move_extent_to_coordgrid, _get_validated_tgt_res
+from pyresample import AreaDefinition
+
+
+class SensorMapGeometryTransformer3D(object):
+    def __init__(self,
+                 lons: np.ndarray,
+                 lats: np.ndarray,
+                 resamp_alg: str = 'nearest',
+                 radius_of_influence: int = 30,
+                 mp_alg: str = 'auto',
+                 **opts) -> None:
+        """Get an instance of SensorMapGeometryTransformer.
+
+        :param lons:    3D longitude array corresponding to the 3D sensor geometry array
+        :param lats:    3D latitude array corresponding to the 3D sensor geometry array
+
+        :Keyword Arguments:  (further documentation here: https://pyresample.readthedocs.io/en/latest/swath.html)
+            - resamp_alg:           resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
+            - radius_of_influence:  <float> Cut off distance in meters (default: 30)
+                                    NOTE: keyword is named 'radius' in case of bilinear resampling
+            - mp_alg                multiprocessing algorithm
+                                    'bands': parallelize over bands using multiprocessing lib
+                                    'tiles': parallelize over tiles using OpenMP
+                                    'auto': automatically choose the algorithm
+            - sigmas:               <list of floats or float> [ONLY 'gauss'] List of sigmas to use for the gauss
+                                    weighting of each channel 1 to k, w_k = exp(-dist^2/sigma_k^2). If only one channel
+                                    is resampled sigmas is a single float value.
+            - neighbours:           <int> [ONLY 'bilinear', 'gauss'] Number of neighbours to consider for each grid
+                                    point when searching the closest corner points
+            - epsilon:              <float> Allowed uncertainty in meters. Increasing uncertainty reduces execution time
+            - weight_funcs:         <list of function objects or function object> [ONLY 'custom'] List of weight
+                                    functions f(dist) to use for the weighting of each channel 1 to k. If only one
+                                    channel is resampled weight_funcs is a single function object.
+            - fill_value:           <int or None> Set undetermined pixels to this value (default: 0).
+                                    If fill_value is None a masked array is returned with undetermined pixels masked
+            - reduce_data:          <bool> Perform initial coarse reduction of source dataset in order to reduce
+                                    execution time
+            - nprocs:               <int>, Number of processor cores to be used
+            - segments:             <int or None> Number of segments to use when resampling.
+                                    If set to None an estimate will be calculated
+            - with_uncert:          <bool> [ONLY 'gauss' and 'custom'] Calculate uncertainty estimates
+                                    NOTE: resampling function has 3 return values instead of 1: result, stddev, count
+        """
+        # validation
+        if lons.ndim != 3:
+            raise ValueError('Expected a 3D longitude array. Received a %dD array.' % lons.ndim)
+        if lats.ndim != 3:
+            raise ValueError('Expected a 3D latitude array. Received a %dD array.' % lats.ndim)
+        if lons.shape != lats.shape:
+            raise ValueError((lons.shape, lats.shape), "'lons' and 'lats' are expected to have the same shape.")
+
+        self.lats = lats
+        self.lons = lons
+        self.resamp_alg = resamp_alg
+        self.radius_of_influence = radius_of_influence
+        self.opts = opts
+
+        # define number of CPUs to use (but avoid sub-multiprocessing)
+        #   -> parallelize either over bands or over image tiles
+        #      bands: multiprocessing uses multiprocessing.Pool, implemented in to_map_geometry / to_sensor_geometry
+        #      tiles: multiprocessing uses OpenMP implemented in pykdtree which is used by pyresample
+        self.opts['nprocs'] = opts.get('nprocs', multiprocessing.cpu_count())
+        self.mp_alg = ('bands' if self.lons.shape[2] >= opts['nprocs'] else 'tiles') if mp_alg == 'auto' else mp_alg
+
+    @staticmethod
+    def _to_map_geometry_2D(kwargs_dict: dict
+                            ) -> Tuple[np.ndarray, tuple, str, int]:
+        assert [var is not None for var in (_global_shared_lons, _global_shared_lats, _global_shared_data)]
+
+        SMGT2D = SensorMapGeometryTransformer(lons=_global_shared_lons[:, :, kwargs_dict['band_idx']],
+                                              lats=_global_shared_lats[:, :, kwargs_dict['band_idx']],
+                                              resamp_alg=kwargs_dict['resamp_alg'],
+                                              radius_of_influence=kwargs_dict['radius_of_influence'],
+                                              **kwargs_dict['init_opts'])
+        data_mapgeo, out_gt, out_prj = SMGT2D.to_map_geometry(data=_global_shared_data[:, :, kwargs_dict['band_idx']],
+                                                              area_definition=kwargs_dict['area_definition'])
+
+        return data_mapgeo, out_gt, out_prj, kwargs_dict['band_idx']
+
+    def _get_common_target_extent(self,
+                                  tgt_epsg: int,
+                                  tgt_coordgrid: Tuple[Tuple, Tuple] = None):
+        if tgt_epsg == 4326:
+            corner_coords_ll = [[self.lons[0, 0, :].min(), self.lats[0, 0, :].max()],  # common UL_xy
+                                [self.lons[0, -1, :].max(), self.lats[0, -1, :].max()],  # common UR_xy
+                                [self.lons[-1, 0, :].min(), self.lats[-1, 0, :].min()],  # common LL_xy
+                                [self.lons[-1, -1, :].max(), self.lats[-1, -1, :].min()],  # common LR_xy
+                                ]
+            common_tgt_extent = _corner_coords_lonlat_to_extent(corner_coords_ll, tgt_epsg)
+        else:
+            # get Lon/Lat corner coordinates of geolayers
+            UL_UR_LL_LR_ll = [(self.lons[y, x], self.lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, 0), (-1, -1)]]
+
+            # transform them to target projection
+            UL_UR_LL_LR_prj = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in UL_UR_LL_LR_ll]
+
+            # separate X and Y
+            X_prj, Y_prj = zip(*UL_UR_LL_LR_prj)
+
+            # 3D geolayers, i.e., the corner coordinates have multiple values for multiple bands
+            # -> use the outermost coordinates to be sure all data is included
+            X_prj = (X_prj[0].min(), X_prj[1].max(), X_prj[2].min(), X_prj[3].max())
+            Y_prj = (Y_prj[0].max(), Y_prj[1].max(), Y_prj[2].min(), Y_prj[3].min())
+
+            # get the extent
+            common_tgt_extent = (min(X_prj), min(Y_prj), max(X_prj), max(Y_prj))
+
+        if tgt_coordgrid:
+            common_tgt_extent = _move_extent_to_coordgrid(common_tgt_extent, *tgt_coordgrid)
+
+        return common_tgt_extent
+
+    def _get_common_area_definition(self,
+                                    data: np.ndarray,
+                                    tgt_prj: Union[str, int],
+                                    tgt_extent: Tuple[float, float, float, float] = None,
+                                    tgt_res: Tuple = None,
+                                    tgt_coordgrid: Tuple[Tuple, Tuple] = None
+                                    ) -> AreaDefinition:
+        # get common target extent
+        tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
+        tgt_extent = tgt_extent or self._get_common_target_extent(tgt_epsg, tgt_coordgrid=tgt_coordgrid)
+
+        SMGT2D = SensorMapGeometryTransformer(lons=self.lons[:, :, 0],  # does NOT affect the computed area definition
+                                              lats=self.lats[:, :, 0],  # -> only needed for __init__
+                                              resamp_alg=self.resamp_alg,
+                                              radius_of_influence=self.radius_of_influence,
+                                              **self.opts)
+        common_area_definition = SMGT2D.compute_areadefinition_sensor2map(data=data[:, :, 0],
+                                                                          tgt_prj=tgt_prj,
+                                                                          tgt_extent=tgt_extent,
+                                                                          tgt_res=tgt_res)
+
+        return common_area_definition
+
+    def to_map_geometry(self,
+                        data: np.ndarray,
+                        tgt_prj: Union[str, int],
+                        tgt_extent: Tuple[float, float, float, float] = None,
+                        tgt_res: Tuple[float, float] = None,
+                        tgt_coordgrid: Tuple[Tuple, Tuple] = None,
+                        area_definition: AreaDefinition = None
+                        ) -> Tuple[np.ndarray, tuple, str]:
+        """Transform the input sensor geometry array into map geometry.
+
+        :param data:            3D numpy array (representing sensor geometry) to be warped to map geometry
+        :param tgt_prj:         target projection (WKT or 'epsg:1234' or <EPSG_int>)
+        :param tgt_extent:      extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
+        :param tgt_res:         target X/Y resolution (e.g., (30, 30))
+        :param tgt_coordgrid:   target coordinate grid ((x, x), (y, y)):
+                                if given, the output extent is moved to this coordinate grid
+        :param area_definition: an instance of pyresample.geometry.AreaDefinition;
+                                OVERRIDES tgt_prj, tgt_extent, tgt_res and tgt_coordgrid; saves computation time
+        """
+        if data.ndim != 3:
+            raise ValueError(data.ndim, "'data' must have 3 dimensions.")
+
+        if data.shape != self.lons.shape:
+            raise ValueError(data.shape, 'Expected a sensor geometry data array with %d rows, %d columns and %d bands.'
+                             % self.lons.shape)
+
+        if not tgt_prj and not area_definition:
+            raise ValueError(tgt_prj, 'Target projection must be given if area_definition is not given.')
+
+        if tgt_coordgrid:
+            tgt_res = _get_validated_tgt_res(tgt_coordgrid, tgt_res)
+
+        init_opts = self.opts.copy()
+        if self.mp_alg == 'bands':
+            del init_opts['nprocs']  # avoid sub-multiprocessing
+
+        # get common area_definition
+        if not area_definition:
+            area_definition = self._get_common_area_definition(data, tgt_prj, tgt_extent, tgt_res, tgt_coordgrid)
+
+        args = [dict(
+            resamp_alg=self.resamp_alg,
+            radius_of_influence=self.radius_of_influence,
+            init_opts=init_opts,
+            area_definition=area_definition,
+            band_idx=band
+        ) for band in range(data.shape[2])]
+
+        if self.opts['nprocs'] > 1 and self.mp_alg == 'bands':
+            with multiprocessing.Pool(self.opts['nprocs'],
+                                      initializer=_initializer,
+                                      initargs=(self.lats, self.lons, data)) as pool:
+                result = pool.map(self._to_map_geometry_2D, args)
+        else:
+            _initializer(self.lats, self.lons, data)
+            result = [self._to_map_geometry_2D(argsdict) for argsdict in args]
+
+        band_inds = list(np.array(result)[:, -1])
+        data_mapgeo = np.dstack([result[band_inds.index(i)][0] for i in range(data.shape[2])])
+        out_gt = result[0][1]
+        out_prj = result[0][2]
+
+        return data_mapgeo, out_gt, out_prj
+
+    @staticmethod
+    def _to_sensor_geometry_2D(kwargs_dict: dict
+                               ) -> (np.ndarray, int):
+        assert [var is not None for var in (_global_shared_lons, _global_shared_lats, _global_shared_data)]
+
+        SMGT2D = SensorMapGeometryTransformer(lons=_global_shared_lons[:, :, kwargs_dict['band_idx']],
+                                              lats=_global_shared_lats[:, :, kwargs_dict['band_idx']],
+                                              resamp_alg=kwargs_dict['resamp_alg'],
+                                              radius_of_influence=kwargs_dict['radius_of_influence'],
+                                              **kwargs_dict['init_opts'])
+        data_sensorgeo = SMGT2D.to_sensor_geometry(data=_global_shared_data[:, :, kwargs_dict['band_idx']],
+                                                   src_prj=kwargs_dict['src_prj'],
+                                                   src_extent=kwargs_dict['src_extent'])
+
+        return data_sensorgeo, kwargs_dict['band_idx']
+
+    def to_sensor_geometry(self,
+                           data: np.ndarray,
+                           src_prj: Union[str, int],
+                           src_extent: Tuple[float, float, float, float]
+                           ) -> np.ndarray:
+        """Transform the input map geometry array into sensor geometry
+
+        :param data:        3D numpy array (representing map geometry) to be warped to sensor geometry
+        :param src_prj:     projection of the input map geometry array (WKT or 'epsg:1234' or <EPSG_int>)
+        :param src_extent:  extent coordinates of input map geometry array (LL_x, LL_y, UR_x, UR_y) in the src_prj
+        """
+        if data.ndim != 3:
+            raise ValueError(data.ndim, "'data' must have 3 dimensions.")
+
+        init_opts = self.opts.copy()
+        if self.mp_alg == 'bands':
+            del init_opts['nprocs']  # avoid sub-multiprocessing
+
+        args = [dict(
+            resamp_alg=self.resamp_alg,
+            radius_of_influence=self.radius_of_influence,
+            init_opts=init_opts,
+            src_prj=src_prj,
+            src_extent=src_extent,
+            band_idx=band
+        ) for band in range(data.shape[2])]
+
+        if self.opts['nprocs'] > 1 and self.mp_alg == 'bands':
+            with multiprocessing.Pool(self.opts['nprocs'],
+                                      initializer=_initializer,
+                                      initargs=(self.lats, self.lons, data)) as pool:
+                result = pool.map(self._to_sensor_geometry_2D, args)
+        else:
+            _initializer(self.lats, self.lons, data)
+            result = [self._to_sensor_geometry_2D(argsdict) for argsdict in args]
+
+        band_inds = list(np.array(result)[:, -1])
+        data_sensorgeo = np.dstack([result[band_inds.index(i)][0] for i in range(data.shape[2])])
+
+        return data_sensorgeo
+
+
+_global_shared_lats = None
+_global_shared_lons = None
+_global_shared_data = None
+
+
+def _initializer(lats, lons, data):
+    """Declare global variables needed for SensorMapGeometryTransformer3D.to_map_geometry and to_sensor_geometry.
+
+    :param lats:
+    :param lons:
+    :param data:
+    """
+    global _global_shared_lats, _global_shared_lons, _global_shared_data
+    _global_shared_lats = lats
+    _global_shared_lons = lons
+    _global_shared_data = data

setup.py

@@ -53,7 +53,6 @@ setup(
         'License :: OSI Approved :: GNU General Public License v3 (GPLv3)',
         'Natural Language :: English',
         'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.5',
         'Programming Language :: Python :: 3.6',
         'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',

tests/test_sensormapgeo.py

@@ -47,20 +47,25 @@ class Test_SensorMapGeometryTransformer(TestCase):
         self.dem_sensor_geo = LoadFile(os.path.join(tests_path, 'data', 'dem_sensor_geo.bsq'))
         self.lons = LoadFile(os.path.join(tests_path, 'data', 'lons_full_vnir.bsq'))
         self.lats = LoadFile(os.path.join(tests_path, 'data', 'lats_full_vnir.bsq'))
-        self.dem_area_extent_coarse_subset_utm = [622613.864409047,  # LL_x
+        self.dem_area_extent_coarse_subset_utm = (622613.864409047,  # LL_x
                                                   5254111.40255343,  # LL_x
                                                   660473.864409047,  # LL_x
-                                                  5269351.40255343]  # UR_y
+                                                  5269351.40255343)  # UR_y
 
-        self.expected_dem_area_extent_lonlat = [10.685733901515151,  # LL_x
+        self.expected_dem_area_extent_lonlat = (10.685733901515151,  # LL_x
                                                 47.44113415492957,  # LL_y
                                                 11.073066098484848,  # UR_x
-                                                47.54576584507042]  # UR_y
+                                                47.54576584507042)  # UR_y
 
-        self.expected_dem_area_extent_utm = [626938.928052,  # LL_x
+        self.expected_dem_area_extent_utm = (626938.928052,  # LL_x
                                              5256253.56579,  # LL_y
                                              656188.928052,  # UR_x
-                                             5267203.56579]  # UR_y
+                                             5267203.56579)  # UR_y
+
+        self.expected_dem_area_extent_utm_ongrid = (626910,  # LL_x
+                                                    5256240,  # LL_y
+                                                    656190,  # UR_x
+                                                    5267220)  # UR_y
 
     def test_to_sensor_geometry(self):
         for rsp_alg in rsp_algs:
@@ -96,6 +101,11 @@ class Test_SensorMapGeometryTransformer(TestCase):
 
             # to Lon/Lat
             dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=4326)
+
+            # from geoarray import GeoArray
+            # GeoArray(dem_map_geo, dem_gt, dem_prj)\
+            #     .save(os.path.join(tests_path, 'test_output', 'resampled_pyresample_ll.bsq'))
+
             self.assertIsInstance(dem_map_geo, np.ndarray)
             self.assertEqual(dem_map_geo.shape, (SMGT.area_definition.height,
                                                  SMGT.area_definition.width))
@@ -120,14 +130,22 @@ class Test_SensorMapGeometryTransformer(TestCase):
                                                 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))
+            # dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632, tgt_res=(30, 30))
+            dem_map_geo, dem_gt, dem_prj = SMGT.to_map_geometry(self.dem_sensor_geo, tgt_prj=32632,
+                                                                tgt_res=(30, 30),
+                                                                # tgt_extent=self.expected_dem_area_extent_utm,
+                                                                tgt_coordgrid=((0, 30), (0, 30))
+                                                                )
+            # from geoarray import GeoArray
+            # GeoArray(dem_map_geo, dem_gt, dem_prj).save(os.path.join(tests_path, 'test_output', 'resampled_gdal.bsq'))
+
             self.assertIsInstance(dem_map_geo, np.ndarray)
-            self.assertEqual(dem_map_geo.shape, (365, 975))
+            self.assertEqual(dem_map_geo.shape, (366, 976))
             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)))
+                                                    np.array(self.expected_dem_area_extent_utm_ongrid)))
             self.assertFalse(np.array_equal(np.unique(dem_map_geo), np.array([0])))
             self.assertTrue(np.isclose(np.mean(dem_map_geo[dem_map_geo != 0]),
                                        np.mean(self.dem_sensor_geo),
@@ -144,23 +162,30 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
 
         self.data_map_geo_3D = np.dstack([dem_map_geo, dem_map_geo])
         self.data_sensor_geo_3D = np.dstack([dem_sensor_geo, dem_sensor_geo])
-        self.lons_3D = np.dstack([lons, lons])  # TODO use different lons per band here
-        self.lats_3D = np.dstack([lats, lats])  # TODO use different lats per band here
+        self.lons_3D = np.dstack([lons, lons + .002])  # assume slightly different coordinates in both bands
+        self.lats_3D = np.dstack([lats, lats + .002])
 
-        self.dem_area_extent_coarse_subset_utm = [622613.864409047,  # LL_x
+        self.dem_area_extent_coarse_subset_utm = (622613.864409047,  # LL_x
                                                   5254111.40255343,  # LL_x
                                                   660473.864409047,  # LL_x
-                                                  5269351.40255343]  # UR_y
+                                                  5269351.40255343)  # UR_y
 
-        self.expected_dem_area_extent_lonlat = [10.685733901515151,  # LL_x
+        # this differs from the 2D version because the geolayer in the second band has slightly different coordinates
+        self.expected_dem_area_extent_lonlat = (10.685733901515151,  # LL_x
                                                 47.44113415492957,  # LL_y
-                                                11.073066098484848,  # UR_x
-                                                47.54576584507042]  # UR_y
+                                                11.075064739115845,  # UR_x
+                                                47.54772559829233)  # UR_y
 
-        self.expected_dem_area_extent_utm = [626938.928052,  # LL_x
+        self.expected_dem_area_extent_utm = (626938.928052,  # LL_x
                                              5256253.56579,  # LL_y
                                              656188.928052,  # UR_x
-                                             5267203.56579]  # UR_y
+                                             5267203.56579)  # UR_y
+
+        # this differs from the 2D version because the geolayer in the second band has slightly different coordinates
+        self.expected_dem_area_extent_utm_ongrid = (626910,  # LL_x
+                                                    5256240,  # LL_y
+                                                    656340,  # UR_x
+                                                    5267430)  # UR_y
 
     def test_to_map_geometry_lonlat_3D_geolayer(self):
         for rsp_alg in rsp_algs:
@@ -172,11 +197,15 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
 
             # to Lon/Lat
             data_mapgeo_3D, dem_gt, dem_prj = SMGT.to_map_geometry(self.data_sensor_geo_3D, tgt_prj=4326)
+
+            # from geoarray import GeoArray
+            # GeoArray(data_mapgeo_3D, dem_gt, dem_prj)\
+            #     .save(os.path.join(tests_path, 'test_output', 'resampled_3D_02_ll.bsq'))
+
             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.assertGreater(data_mapgeo_3D.shape[0], self.data_sensor_geo_3D.shape[0])
-            self.assertGreater(data_mapgeo_3D.shape[1], self.data_sensor_geo_3D.shape[1])
+            self.assertGreater(np.dot(*data_mapgeo_3D.shape[:2]), np.dot(*self.data_sensor_geo_3D.shape[:2]))
             self.assertEqual(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],
@@ -189,6 +218,41 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
                                        rtol=0.01))
             self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)
 
+    def test_to_map_geometry_utm_3D_geolayer(self):
+        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=32632,
+                                                                   tgt_res=(30, 30),
+                                                                   # tgt_extent=self.expected_dem_area_extent_utm,
+                                                                   tgt_coordgrid=((0, 30), (0, 30))
+                                                                   )
+            # from geoarray import GeoArray
+            # GeoArray(data_mapgeo_3D, dem_gt, dem_prj)\
+            #     .save(os.path.join(tests_path, 'test_output', 'resampled_3D_02_pyresample.bsq'))
+
+            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.assertGreater(np.dot(*data_mapgeo_3D.shape[:2]), np.dot(*self.data_sensor_geo_3D.shape[:2]))
+            self.assertEqual(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_utm_ongrid)))
+
+            self.assertTrue(np.isclose(np.mean(data_mapgeo_3D[data_mapgeo_3D != 0]),
+                                       np.mean(self.data_sensor_geo_3D),
+                                       rtol=0.01))
+            self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)
+
     def test_to_sensor_geometry(self):
         for rsp_alg in rsp_algs:
             SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,