First draft to avoid deadlock caused by an incompatible build of OpenMP.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

sensormapgeo/transformer_2d.py

@@ -100,9 +100,8 @@ class SensorMapGeometryTransformer(object):
 
         # NOTE: If pykdtree is built with OpenMP support (default) the number of threads is controlled with the
         #       standard OpenMP environment variable OMP_NUM_THREADS. The nprocs argument has no effect on pykdtree.
+        os.environ['OMP_NUM_THREADS'] = '%d' % (self.opts['nprocs'] if 'nprocs' in self.opts else 1)
         if 'nprocs' in self.opts:
-            if self.opts['nprocs'] > 1:
-                os.environ['OMP_NUM_THREADS'] = '%d' % opts['nprocs']
             del self.opts['nprocs']
 
         self.lats = lats

sensormapgeo/transformer_3d.py

@@ -26,7 +26,12 @@
 
 from typing import Union, Tuple
 import multiprocessing
+from timeout_decorator import timeout
+from concurrent.futures import TimeoutError as Timeout
+from warnings import warn
+import platform
 
+from pebble import ProcessPool, ProcessExpired
 import numpy as np
 from pyproj import CRS
 from py_tools_ds.geo.coord_trafo import transform_any_prj
@@ -224,8 +229,43 @@ class SensorMapGeometryTransformer3D(object):
             # NOTE: We use imap here as it directly returns the results when available (works like a generator).
             #       This saves a lot of memory compared with map. We also don't use an initializer to share the input
             #       arrays because this would allocate the memory for the input arrays of all bands at once.
-            with multiprocessing.Pool(self.opts['nprocs']) as pool:
-                result = [res for res in pool.imap_unordered(self._to_map_geometry_2D, args)]
+            try:
+                with ProcessPool() as pool:
+                    future = pool.map(self._to_map_geometry_2D, args, timeout=10)
+                    result = [i for i in future.result()]
+
+            except (Timeout, ProcessExpired):
+                msg = "Switched multiprocessing algorithm from 'bands' to 'tiles' due to a timeout in 'bands' mode. "
+                if platform.system() == 'Linux':
+                    msg += "Consider using the LLVM instead of the GNU build of OpenMP to fix this issue (install, " \
+                           "e.g., by 'conda install -c conda-forge _openmp_mutex=*=1_llvm'."
+                warn(msg, RuntimeWarning)
+
+                init_opts = self.opts.copy()
+                for argset in args:
+                    argset.update(dict(init_opts=init_opts))
+
+                result = [self._to_map_geometry_2D(argsdict) for argsdict in args]
+
+            # from time import time
+            # t0 = time()
+            # self._to_map_geometry_2D(args[0])
+            # print(args[0]['data'].shape, time() - t0)
+            #
+            # ncols, nbands = data.shape[0], data.shape[2]
+            # nprocs = self.opts['nprocs'] if nbands > self.opts['nprocs'] else nbands
+            # max_seconds = None if self.mp_alg == 'tiles' else \
+            #     (0.002 * ncols * nbands / nprocs * 10)
+            # print('timeout=%s' % max_seconds)
+            # @timeout(seconds=max_seconds, timeout_exception=TimeoutError)
+            # def run():
+            #     with multiprocessing.Pool(self.opts['nprocs']) as pool:
+            #         return [res for res in pool.imap_unordered(self._to_map_geometry_2D, args)]
+            #
+            # try:
+            #     result = run()
+            # except TimeoutError:
+            #     print('CATCHED!')
 
         else:
             result = [self._to_map_geometry_2D(argsdict) for argsdict in args]

tests/test_sensormapgeo.py

@@ -39,6 +39,7 @@ from sensormapgeo import SensorMapGeometryTransformer, SensorMapGeometryTransfor
 
 tests_path = os.path.abspath(os.path.join(__path__[0], "..", "tests"))
 rsp_algs = ['nearest', 'bilinear', 'gauss']
+mp_algs = ['bands', 'tiles']
 
 
 class Test_SensorMapGeometryTransformer(TestCase):
@@ -194,10 +195,12 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
 
     def test_to_map_geometry_lonlat_3D_geolayer(self):
         for rsp_alg in rsp_algs:
+            for mp_alg in mp_algs:
                 SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                       lats=self.lats_3D,
                                                       # resamp_alg='nearest',
                                                       resamp_alg=rsp_alg,
+                                                      mp_alg=mp_alg
                                                       )
 
                 # to Lon/Lat
@@ -223,12 +226,45 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
                                            rtol=0.01))
                 self.assertEqual(self.data_sensor_geo_3D.dtype, data_mapgeo_3D.dtype)
 
+    # def test_to_map_geometry_lonlat_3D_geolayer_bands(self):
+    #     SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
+    #                                           lats=self.lats_3D,
+    #                                           # resamp_alg='nearest',
+    #                                           resamp_alg='gauss',
+    #                                           mp_alg='bands'
+    #                                           )
+    #
+    #     # 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(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_lonlat)))
+    #
+    #     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_map_geometry_utm_3D_geolayer(self):
         for rsp_alg in rsp_algs:
+            for mp_alg in mp_algs:
                 SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                       lats=self.lats_3D,
                                                       # resamp_alg='nearest',
                                                       resamp_alg=rsp_alg,
+                                                      mp_alg=mp_alg
                                                       )
 
                 # to Lon/Lat
@@ -260,9 +296,11 @@ class Test_SensorMapGeometryTransformer3D(TestCase):
 
     def test_to_sensor_geometry(self):
         for rsp_alg in rsp_algs:
+            for mp_alg in mp_algs:
                 SMGT = SensorMapGeometryTransformer3D(lons=self.lons_3D,
                                                       lats=self.lats_3D,
                                                       resamp_alg=rsp_alg,
+                                                      mp_alg=mp_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)