SensorMapGeometryTransformer3D.to_map_geometry() now computes a common area...

SensorMapGeometryTransformer3D.to_map_geometry() now computes a common area definition only ONCE which saves computation time and increases stability.
Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

sensormapgeo/transformer_3d.py

@@ -32,6 +32,7 @@ from py_tools_ds.geo.projection import WKT2EPSG, proj4_to_WKT
 from py_tools_ds.geo.coord_trafo import get_proj4info
 
 from .transformer_2d import SensorMapGeometryTransformer, _corner_coords_lonlat_to_extent
+from pyresample import AreaDefinition
 
 
 class SensorMapGeometryTransformer3D(object):
@@ -105,9 +106,7 @@ class SensorMapGeometryTransformer3D(object):
                                               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']],
-                                                              tgt_prj=kwargs_dict['tgt_prj'],
-                                                              tgt_extent=kwargs_dict['tgt_extent'],
-                                                              tgt_res=kwargs_dict['tgt_res'])
+                                                              area_definition=kwargs_dict['area_definition'])
 
         return data_mapgeo, out_gt, out_prj, kwargs_dict['band_idx']
 
@@ -121,17 +120,41 @@ class SensorMapGeometryTransformer3D(object):
 
         return 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) -> 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)
+
+        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 = None) -> Tuple[np.ndarray, tuple, str]:
+                        tgt_res: 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 area_definition: an instance of pyresample.geometry.AreaDefinition;
+                                OVERRIDES tgt_prj, tgt_extent and tgt_res; saves computation time
         """
         if data.ndim != 3:
             raise ValueError(data.ndim, "'data' must have 3 dimensions.")
@@ -140,21 +163,19 @@ class SensorMapGeometryTransformer3D(object):
             raise ValueError(data.shape, 'Expected a sensor geometry data array with %d rows, %d columns and %d bands.'
                              % self.lons.shape)
 
-        # 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)
-
         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)
+
         args = [dict(
             resamp_alg=self.resamp_alg,
             radius_of_influence=self.radius_of_influence,
             init_opts=init_opts,
-            tgt_prj=tgt_prj,
-            tgt_extent=tgt_extent,
-            tgt_res=tgt_res,
+            area_definition=area_definition,
             band_idx=band
         ) for band in range(data.shape[2])]