Merge pull request #21 from ec-jrc/development

Added full correlation distance decay (CDD) interpolation method, using a map of CDD values.
Speeds up angular distance weighted (ADW) interpolation by using numba.
Added option to split computation in sub-regions to save memory

.github/workflows/ci_env.yml

@@ -33,7 +33,7 @@ jobs:
       run: |
         conda list
         pip install poppler-utils
-        conda install 'poppler=>22.10.0'
+        conda install -c conda-forge 'poppler=>22.10.0'
     - name: Check installation
       shell: bash -el {0}
       run: |

.gitignore

@@ -5,7 +5,8 @@
 /build
 pyg2p.egg-info
 /dist
-
+*.nbc
+*.nbi
 /.python-version
 .coverage
 .coverage.*

README.md

@@ -670,16 +670,44 @@ Attributes p, leafsize and eps for the kd tree algorithm are default in scipy li
 | leafsize  | 10                   |
 
 #### ADW
-It's the Angular Distance Weighted (ADW) algorithm with scipy.kd_tree, using 4 neighbours.
-If @adw_broadcasting is set to true, computations will run in full broadcasting mode but requires more memory
+It's the Angular Distance Weighted (ADW) algorithm by Shepard et al. 1968, with scipy.kd_tree using 11 neighbours.
+If @use_broadcasting is set to true, computations will run in full broadcasting mode but requires more memory
+If @num_of_splits is set to any number, computations will be split on subset and then recollected into the final map, to save mamory (do not set it if you have enought memory to run interpolation)
 
 ```json
 {
 "Interpolation": {
   "@latMap": "/dataset/maps/europe5km/lat.map",
   "@lonMap": "/dataset/maps/europe5km/long.map",
   "@mode": "adw",
-  "@adw_broadcasting": false}
+  "@use_broadcasting": false,
+  "@num_of_splits": 10}
+}
+```
+
+#### CDD
+It's the Correlation Distance Decay (CDD) modified implementation of the Angular Distance Weighted algorithm, with scipy.kd_tree using 11 neighbours. It needs a map of CDD values for each point, to be specified in the field @cdd_map
+@cdd_mode can be one of the following values: "Hofstra", "NewEtAl" or "MixHofstraShepard"
+In case of mode "MixHofstraShepard", @cdd_options allows to customize the parameters of Hofstra and Shepard algorithm ("weights_mode": can be "All" or "OnlyTOP10" to take 10 higher values only in the interpolation of each point).
+If @use_broadcasting is set to true, computations will run in full broadcasting mode but requires more memory
+If @num_of_splits is set to any number, computations will be split on subset and then recollected into the final map, to save mamory (do not set it if you have enought memory to run interpolation)
+
+```json
+{
+"Interpolation": {
+  "@latMap": "/dataset/maps/europe5km/lat.map",
+  "@lonMap": "/dataset/maps/europe5km/long.map",
+  "@mode": "cdd",
+  "@cdd_map": "/dataset/maps/europe5km/cdd_map.nc",
+  "@cdd_mode": "MixHofstraShepard",
+  "@cdd_options": {
+    "m_const": 4,
+    "min_num_of_station": 4,
+    "radius_ratio": 0.3333333333333333,
+    "weights_mode": "All"
+  },
+  "@use_broadcasting": false,
+  "@num_of_splits": 10}
 }
 ```
 

requirements.txt

@@ -1,7 +1,8 @@
-netCDF4>=1.5.3
+netCDF4>=1.5.3,<=1.6.4
 ujson>=5.4.0
 numpy>=1.18.2
-scipy>=1.4.1
+numba>=0.54.0
+scipy>=1.6.0
 # GDAL
 numexpr>=2.7.1
 importlib-metadata<5.0.0

src/pyg2p/VERSION

@@ -1 +1 @@
-3.2.6
+3.2.7

src/pyg2p/main/api.py

@@ -168,7 +168,11 @@ def _define_exec_params(self):
         self._vars['outMaps.clone'] = self.api_conf['OutMaps']['cloneMap']
         interpolation_conf = self.api_conf['OutMaps']['Interpolation']
         self._vars['interpolation.mode'] = interpolation_conf.get('mode', self.default_values['interpolation.mode'])
+        self._vars['interpolation.cdd_map'] = interpolation_conf.get('cdd_map', '')
+        self._vars['interpolation.cdd_mode'] = interpolation_conf.get('cdd_mode', '')
+        self._vars['interpolation.cdd_options'] = interpolation_conf.get('cdd_options', None)
         self._vars['interpolation.use_broadcasting'] = interpolation_conf.get('use_broadcasting', False)
+        self._vars['interpolation.num_of_splits'] = interpolation_conf.get('num_of_splits', None)        
         self._vars['interpolation.rotated_target'] = interpolation_conf.get('rotated_target', False)
         if not self._vars['interpolation.dir'] and self.api_conf.get('intertableDir'):
             self._vars['interpolation.dirs']['user'] = self.api_conf['intertableDir']

src/pyg2p/main/context.py

@@ -360,7 +360,11 @@ def _define_exec_params(self):
         self._vars['outMaps.clone'] = exec_conf['OutMaps']['@cloneMap']
         interpolation_conf = exec_conf['OutMaps']['Interpolation']
         self._vars['interpolation.mode'] = interpolation_conf.get('@mode', self.default_values['interpolation.mode'])
-        self._vars['interpolation.adw_broadcasting'] = interpolation_conf.get('@adw_broadcasting', False)
+        self._vars['interpolation.cdd_map'] = interpolation_conf.get('@cdd_map', '')
+        self._vars['interpolation.cdd_mode'] = interpolation_conf.get('@cdd_mode', '')
+        self._vars['interpolation.cdd_options'] = interpolation_conf.get('@cdd_options', None)
+        self._vars['interpolation.use_broadcasting'] = interpolation_conf.get('@use_broadcasting', False)
+        self._vars['interpolation.num_of_splits'] = interpolation_conf.get('@num_of_splits', None)
         self._vars['interpolation.rotated_target'] = interpolation_conf.get('@rotated_target', False)
         if not self._vars['interpolation.dir'] and interpolation_conf.get('@intertableDir'):
             # get from JSON

src/pyg2p/main/interpolation/__init__.py

@@ -20,7 +20,7 @@
 class Interpolator(Loggable):
     _LOADED_INTERTABLES = {}
     _prefix = 'I'
-    scipy_modes_nnear = {'nearest': 1, 'invdist': 4, 'adw': 4, 'cdd': 4, 'bilinear': 4, 'triangulation': 3, 'bilinear_delaunay': 4}
+    scipy_modes_nnear = {'nearest': 1, 'invdist': 4, 'adw': 11, 'cdd': 11, 'bilinear': 4, 'triangulation': 3, 'bilinear_delaunay': 4}
     suffixes = {'grib_nearest': 'grib_nearest', 'grib_invdist': 'grib_invdist',
                 'nearest': 'scipy_nearest', 'invdist': 'scipy_invdist', 'adw': 'scipy_adw', 'cdd': 'scipy_cdd',
                 'bilinear': 'scipy_bilinear', 'triangulation': 'scipy_triangulation', 'bilinear_delaunay': 'scipy_bilinear_delaunay'}
@@ -31,7 +31,11 @@ def __init__(self, exec_ctx, mv_input):
         self._mv_grib = mv_input
         self.interpolate_with_grib = exec_ctx.is_with_grib_interpolation
         self._mode = exec_ctx.get('interpolation.mode')
-        self._adw_broadcasting = exec_ctx.get('interpolation.adw_broadcasting')
+        self._cdd_map = exec_ctx.get('interpolation.cdd_map')
+        self._cdd_mode = exec_ctx.get('interpolation.cdd_mode')
+        self._cdd_options = exec_ctx.get('interpolation.cdd_options')
+        self._use_broadcasting = exec_ctx.get('interpolation.use_broadcasting')
+        self._num_of_splits = exec_ctx.get('interpolation.num_of_splits')
         self._source_filename = pyg2p.util.files.filename(exec_ctx.get('input.file'))
         self._suffix = self.suffixes[self._mode]
         self._intertable_dirs = exec_ctx.get('interpolation.dirs')
@@ -167,6 +171,10 @@ def interpolate_scipy(self, latgrib, longrib, v, grid_id, grid_details=None):
                 # Global_3arcmin DEBUG
                 latefas=self._target_coords.lats[1800-int(DEBUG_MAX_LAT*20):1800-int(DEBUG_MIN_LAT*20), 3600+int(DEBUG_MIN_LON*20):3600+int(DEBUG_MAX_LON*20)]
                 lonefas=self._target_coords.lons[1800-int(DEBUG_MAX_LAT*20):1800-int(DEBUG_MIN_LAT*20), 3600+int(DEBUG_MIN_LON*20):3600+int(DEBUG_MAX_LON*20)]
+                #latefas-=0.008369999999992217
+                # lonefas-=0.00851999999999431
+                #lonefas-=0.024519999999977227   
+
             else:
                 # European_1arcmin DEBUG
                 selection_lats = np.logical_and(self._target_coords.lats[:,0]>=DEBUG_MIN_LAT,self._target_coords.lats[:,0]<=DEBUG_MAX_LAT)
@@ -194,9 +202,18 @@ def interpolate_scipy(self, latgrib, longrib, v, grid_id, grid_details=None):
             # latgrib = np.array([ 7.39050803,  8.72151493,  7.82210701,  7.35906546])
             # longrib = np.array([49.16690015, 48.11557968, 48.27217824, 49.70238655])
             # v = np.array([100, 133, 166, 200  ])
-            latgrib = np.array([ 8,  8,  8,  8])
-            longrib = np.array([45, 48.5, 49, 50])
-            v = np.array([200, 100, 100, 100  ])
+            # latgrib = np.array([ 8,  8,  8,  8])
+            # longrib = np.array([45, 48.5, 49, 50])
+            # v = np.array([200, 100, 100, 100  ])
+
+            # OR load data points for the TEST from file
+            #data = np.genfromtxt('/media/sf_VMSharedFolder/pyg2p_adw_cdd_test/pr199106180600_idw.txt', delimiter='\t', skip_header=1)
+            #data = np.genfromtxt('/media/sf_VMSharedFolder/pyg2p_adw_cdd_test/pr199106170600_20230714101901.txt', delimiter='\t', skip_header=1)
+            data = np.genfromtxt('/media/sf_VMSharedFolder/test_split/tn202401010600_20240213140643.txt', delimiter='\t', skip_header=1)
+            longrib = data[:,0]
+            latgrib = data[:,1]
+            v = data[:,2]
+
             intertable_id, intertable_name = 'DEBUG_ADW','DEBUG_ADW.npy'
 
         nnear = self.scipy_modes_nnear[self._mode]
@@ -212,11 +229,22 @@ def interpolate_scipy(self, latgrib, longrib, v, grid_id, grid_details=None):
                 self._log('Trying to interpolate without grib lat/lons. Probably a malformed grib!', 'ERROR')
                 raise ApplicationException.get_exc(5000)
 
+            # CR: to use float32 computations uncomment here:
+            # longrib=np.float32(longrib)
+            # latgrib=np.float32(latgrib)
+            # lonefas=np.float32(lonefas)
+            # latefas=np.float32(latefas)
+            # v=np.float32(v)
+            # self.mv_out=np.float32(self.mv_out)
+
             self._log('\nInterpolating table not found\n Id: {}\nWill create file: {}'.format(intertable_id, intertable_name), 'WARN')
             scipy_interpolation = ScipyInterpolation(longrib, latgrib, grid_details, v.ravel(), nnear, self.mv_out,
                                           self._mv_grib, target_is_rotated=self._rotated_target_grid,
-                                          parallel=self.parallel, mode=self._mode, use_broadcasting=self._adw_broadcasting)
-            _, weights, indexes = scipy_interpolation.interpolate(lonefas, latefas)
+                                          parallel=self.parallel, mode=self._mode, 
+                                          cdd_map=self._cdd_map, cdd_mode=self._cdd_mode, cdd_options = self._cdd_options,
+                                          use_broadcasting=self._use_broadcasting,
+                                          num_of_splits=self._num_of_splits)
+            _, weights, indexes = scipy_interpolation.interpolate(lonefas, latefas)            
             result = self._interpolate_scipy_append_mv(v, weights, indexes, nnear)
 
             # saving interpolation lookup table