From cc28855eb3984929cd27f0b596ec5da2e665d78d Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 24 May 2023 14:52:21 +0200 Subject: [PATCH 01/17] Github action gdal installation fix --- .github/workflow/ci_env.yml | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/.github/workflow/ci_env.yml b/.github/workflow/ci_env.yml index 45c1397..2be6572 100644 --- a/.github/workflow/ci_env.yml +++ b/.github/workflow/ci_env.yml @@ -14,10 +14,15 @@ jobs: with: auto-update-conda: true python-version: 3.7 - - name: Install conda env + - name: Install python and gcc shell: bash -el {0} run: | - conda install -c conda-forge python=3.7 gcc=12.1.0 gdal eccodes + conda install -c conda-forge python=3.7 + conda install -c conda-forge gcc=12.1.0 + - name: Install gdal + shell: bash -el {0} + run: | + conda install -c conda-forge gdal eccodes - name: Install dependencies shell: bash -el {0} run: | From af8053383f9455bf1476369ed959f8bcc46b32a1 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 24 May 2023 14:54:25 +0200 Subject: [PATCH 02/17] Renamed workflow folder --- .github/workflows/ci_env.yml | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 .github/workflows/ci_env.yml diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml new file mode 100644 index 0000000..2be6572 --- /dev/null +++ b/.github/workflows/ci_env.yml @@ -0,0 +1,34 @@ +name: Pyg2p Unit Tests + +on: [push] + +jobs: + tests: + runs-on: ubuntu-20.04 + strategy: + max-parallel: 5 + + steps: + - uses: actions/checkout@v3 + - uses: conda-incubator/setup-miniconda@v2 + with: + auto-update-conda: true + python-version: 3.7 + - name: Install python and gcc + shell: bash -el {0} + run: | + conda install -c conda-forge python=3.7 + conda install -c conda-forge gcc=12.1.0 + - name: Install gdal + shell: bash -el {0} + run: | + conda install -c conda-forge gdal eccodes + - name: Install dependencies + shell: bash -el {0} + run: | + pip install -r requirements.txt + - name: Test with pytest + shell: bash -el {0} + run: | + pip install pytest pytest-cov + pytest \ No newline at end of file From b768242a79bf919f9ebb057d194824625d4490a3 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 24 May 2023 14:54:44 +0200 Subject: [PATCH 03/17] Renamed workflow folder --- .github/workflow/ci_env.yml | 34 ---------------------------------- 1 file changed, 34 deletions(-) delete mode 100644 .github/workflow/ci_env.yml diff --git a/.github/workflow/ci_env.yml b/.github/workflow/ci_env.yml deleted file mode 100644 index 2be6572..0000000 --- a/.github/workflow/ci_env.yml +++ /dev/null @@ -1,34 +0,0 @@ -name: Pyg2p Unit Tests - -on: [push] - -jobs: - tests: - runs-on: ubuntu-20.04 - strategy: - max-parallel: 5 - - steps: - - uses: actions/checkout@v3 - - uses: conda-incubator/setup-miniconda@v2 - with: - auto-update-conda: true - python-version: 3.7 - - name: Install python and gcc - shell: bash -el {0} - run: | - conda install -c conda-forge python=3.7 - conda install -c conda-forge gcc=12.1.0 - - name: Install gdal - shell: bash -el {0} - run: | - conda install -c conda-forge gdal eccodes - - name: Install dependencies - shell: bash -el {0} - run: | - pip install -r requirements.txt - - name: Test with pytest - shell: bash -el {0} - run: | - pip install pytest pytest-cov - pytest \ No newline at end of file From 213efef56ba42a99811aabb5d8d3653b90827f97 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 24 May 2023 16:31:43 +0200 Subject: [PATCH 04/17] Fixed importlib-metadata requirement. Added debug infos in CI github actions --- .github/workflows/ci_env.yml | 6 ++++++ requirements.txt | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 2be6572..268ec5b 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -27,6 +27,12 @@ jobs: shell: bash -el {0} run: | pip install -r requirements.txt + - name: Check installation + shell: bash -el {0} + run: | + gdal-config --version + python -c "from osgeo import gdal; print(gdal.__version__)" + conda list - name: Test with pytest shell: bash -el {0} run: | diff --git a/requirements.txt b/requirements.txt index bf2f40e..e1cbf81 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ numpy>=1.18.2 scipy>=1.4.1 # GDAL numexpr>=2.7.1 -# eccodes-python # replacing old GRIB API +importlib-metadata<5.0.0 # dask packages for parallel interpolation dask[bag] From bcbabca1ccba5b2824fcd7ab57a3659b31545228 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Wed, 24 May 2023 17:10:44 +0200 Subject: [PATCH 05/17] Update ci_env.yml --- .github/workflows/ci_env.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 268ec5b..85c1185 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -27,6 +27,7 @@ jobs: shell: bash -el {0} run: | pip install -r requirements.txt + pip install . - name: Check installation shell: bash -el {0} run: | @@ -37,4 +38,4 @@ jobs: shell: bash -el {0} run: | pip install pytest pytest-cov - pytest \ No newline at end of file + pytest From 5461e9c060339e188d48d2f5571d8c0cd2a094ca Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Wed, 24 May 2023 17:33:15 +0200 Subject: [PATCH 06/17] Update ci_env.yml --- .github/workflows/ci_env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 85c1185..484e07f 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -31,9 +31,9 @@ jobs: - name: Check installation shell: bash -el {0} run: | + conda list gdal-config --version python -c "from osgeo import gdal; print(gdal.__version__)" - conda list - name: Test with pytest shell: bash -el {0} run: | From 837379a11469a02403d6c21827590171c201cb84 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Wed, 24 May 2023 18:03:55 +0200 Subject: [PATCH 07/17] Update ci_env.yml for poppler issue --- .github/workflows/ci_env.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 484e07f..2a59559 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -28,6 +28,12 @@ jobs: run: | pip install -r requirements.txt pip install . + - name: Workaround for poppler lib issue + shell: bash -el {0} + run: | + conda list + pip install poppler-utils + conda install poppler - name: Check installation shell: bash -el {0} run: | From f0b3545932add8c167d8f635222be85eb037c356 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Thu, 25 May 2023 08:47:19 +0200 Subject: [PATCH 08/17] Update ci_env.yml --- .github/workflows/ci_env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 2a59559..a1d0b76 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -33,7 +33,7 @@ jobs: run: | conda list pip install poppler-utils - conda install poppler + conda install poppler=>22.10.0 - name: Check installation shell: bash -el {0} run: | From e00ddc81e90b6db5e6ce9981450bf75be8a147d5 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Thu, 25 May 2023 09:41:28 +0200 Subject: [PATCH 09/17] Update ci_env.yml --- .github/workflows/ci_env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index a1d0b76..1743625 100644 --- a/.github/workflows/ci_env.yml +++ b/.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 'poppler=>22.10.0' - name: Check installation shell: bash -el {0} run: | From a2b938fd3360a13efe6af0ec82d2e027b990e6bb Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 7 Jun 2023 15:23:21 +0200 Subject: [PATCH 10/17] added flexibility to use n!=4 in inverse distance interpolation --- src/pyg2p/main/interpolation/scipy_interpolation_lib.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py index 5c91710..41ec1bf 100644 --- a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py +++ b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py @@ -436,7 +436,9 @@ def _build_weights_invdist(self, distances, indexes, nnear): dist_leq_min_upper_bound = np.logical_and(~dist_leq_1e_10, distances[:, 0] <= self.min_upper_bound) # distances <= 1e-10 : take exactly the point, weight = 1 - weights[dist_leq_1e_10] = np.array([1., 0., 0., 0.]) + onlyfirst_array = np.zeros(nnear) + onlyfirst_array[0] = 1 + weights[dist_leq_1e_10] = onlyfirst_array idxs[dist_leq_1e_10] = indexes[dist_leq_1e_10] result[dist_leq_1e_10] = z[indexes[dist_leq_1e_10][:, 0]] From 2032a1b4a8c0f68cbaac6214ec444743f4906cb2 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 21 Jun 2023 14:46:51 +0200 Subject: [PATCH 11/17] Added Angular distance weighting (ADW) interpolation. Fixed issue with variable order --- src/pyg2p/VERSION | 2 +- src/pyg2p/main/context.py | 2 +- src/pyg2p/main/interpolation/__init__.py | 6 +- .../interpolation/scipy_interpolation_lib.py | 149 +++++++++++++++--- src/pyg2p/main/readers/netcdf.py | 14 +- 5 files changed, 141 insertions(+), 32 deletions(-) diff --git a/src/pyg2p/VERSION b/src/pyg2p/VERSION index 9b7a431..448ada3 100644 --- a/src/pyg2p/VERSION +++ b/src/pyg2p/VERSION @@ -1 +1 @@ -3.2.4 \ No newline at end of file +3.2.5 \ No newline at end of file diff --git a/src/pyg2p/main/context.py b/src/pyg2p/main/context.py index 02dbcf1..14681d3 100644 --- a/src/pyg2p/main/context.py +++ b/src/pyg2p/main/context.py @@ -11,7 +11,7 @@ class Context: - allowed_interp_methods = ('grib_nearest', 'grib_invdist', 'nearest', 'invdist', 'bilinear', 'triangulation', 'bilinear_delaunay') + allowed_interp_methods = ('grib_nearest', 'grib_invdist', 'nearest', 'invdist', 'adw', 'cdd', 'bilinear', 'triangulation', 'bilinear_delaunay') default_values = {'interpolation.mode': 'grib_nearest', 'outMaps.unitTime': '24'} def __getitem__(self, param): diff --git a/src/pyg2p/main/interpolation/__init__.py b/src/pyg2p/main/interpolation/__init__.py index 0f24089..4904e52 100644 --- a/src/pyg2p/main/interpolation/__init__.py +++ b/src/pyg2p/main/interpolation/__init__.py @@ -20,9 +20,9 @@ class Interpolator(Loggable): _LOADED_INTERTABLES = {} _prefix = 'I' - scipy_modes_nnear = {'nearest': 1, 'invdist': 4, 'bilinear': 4, 'triangulation': 3, 'bilinear_delaunay': 4} + scipy_modes_nnear = {'nearest': 1, 'invdist': 4, 'adw': 4, 'cdd': 4, 'bilinear': 4, 'triangulation': 3, 'bilinear_delaunay': 4} suffixes = {'grib_nearest': 'grib_nearest', 'grib_invdist': 'grib_invdist', - 'nearest': 'scipy_nearest', 'invdist': 'scipy_invdist', + 'nearest': 'scipy_nearest', 'invdist': 'scipy_invdist', 'adw': 'scipy_adw', 'cdd': 'scipy_cdd', 'bilinear': 'scipy_bilinear', 'triangulation': 'scipy_triangulation', 'bilinear_delaunay': 'scipy_bilinear_delaunay'} _format_intertable = 'tbl{prognum}_{source_file}_{target_size}_{suffix}.npy.gz'.format @@ -127,7 +127,7 @@ def _read_intertable(self, tbl_fullpath): coeffs = intertable['coeffs'] return indexes[0], indexes[1], indexes[2], indexes[3], indexes[4], indexes[5], coeffs[0], coeffs[1], coeffs[2], coeffs[3] else: - # self._mode in ('invdist', 'nearest', 'bilinear', 'triangulation', 'bilinear_delaunay'): + # self._mode in ('invdist', 'adw', 'cdd', 'nearest', 'bilinear', 'triangulation', 'bilinear_delaunay'): # return indexes and weighted distances (only used with nnear > 1) indexes = intertable['indexes'] coeffs = intertable['coeffs'] diff --git a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py index 41ec1bf..2a01e1a 100644 --- a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py +++ b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py @@ -337,22 +337,23 @@ def interpolate(self, target_lons, target_lats): else: if self.mode == 'invdist': # return results, distances, indexes - result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear) - else: - if self.mode == 'bilinear' and self.nnear == 4: # bilinear interpolation only supported with nnear = 4 - # BILINEAR INTERPOLATION - result, weights, indexes = self._build_weights_bilinear(distances, indexes, efas_locations, self.nnear) - else: - if self.mode == 'triangulation': - # linear barycentric interpolation on Delaunay triangulation - result, weights, indexes = self._build_weights_triangulation(use_bilinear=False) - else: - if self.mode == 'bilinear_delaunay': - # bilinear interpolation on Delaunay triangulation - result, weights, indexes = self._build_weights_triangulation(use_bilinear=True) - else: - raise ApplicationException.get_exc(INVALID_INTERPOL_METHOD, - f"interpolation method not supported (mode = {self.mode}, nnear = {self.nnear})") + result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear) + elif self.mode == 'adw': + result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='Shepard') + elif self.mode == 'cdd': + result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='CDD') + elif self.mode == 'bilinear' and self.nnear == 4: # bilinear interpolation only supported with nnear = 4 + # BILINEAR INTERPOLATION + result, weights, indexes = self._build_weights_bilinear(distances, indexes, efas_locations, self.nnear) + elif self.mode == 'triangulation': + # linear barycentric interpolation on Delaunay triangulation + result, weights, indexes = self._build_weights_triangulation(use_bilinear=False) + elif self.mode == 'bilinear_delaunay': + # bilinear interpolation on Delaunay triangulation + result, weights, indexes = self._build_weights_triangulation(use_bilinear=True) + else: + raise ApplicationException.get_exc(INVALID_INTERPOL_METHOD, + f"interpolation method not supported (mode = {self.mode}, nnear = {self.nnear})") stdout.write('End scipy interpolation: {}\n'.format(now_string())) @@ -419,8 +420,52 @@ def _build_nn(self, distances, indexes): stdout.flush() return result, idxs - # Invdist Optimized version (using broadcast) - def _build_weights_invdist(self, distances, indexes, nnear): + # # Invdist Optimized version (using broadcast) + # def _build_weights_invdist(self, distances, indexes, nnear): + # z = self.z + # result = mask_it(np.empty((len(distances),) + np.shape(z[0])), self._mv_target, 1) + # weights = empty((len(distances),) + (nnear,)) + # idxs = empty((len(indexes),) + (nnear,), fill_value=z.size, dtype=int) + # num_cells = result.size + # back_char, _ = progress_step_and_backchar(num_cells) + + # stdout.write('Skipping neighbors at distance > {}\n'.format(self.min_upper_bound)) + # stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 0, 5, 0, 0)) + # stdout.flush() + + # dist_leq_1e_10 = distances[:, 0] <= 1e-10 + # dist_leq_min_upper_bound = np.logical_and(~dist_leq_1e_10, distances[:, 0] <= self.min_upper_bound) + + # # distances <= 1e-10 : take exactly the point, weight = 1 + # onlyfirst_array = np.zeros(nnear) + # onlyfirst_array[0] = 1 + # weights[dist_leq_1e_10] = onlyfirst_array + # idxs[dist_leq_1e_10] = indexes[dist_leq_1e_10] + # result[dist_leq_1e_10] = z[indexes[dist_leq_1e_10][:, 0]] + + # w = np.zeros_like(distances) + # w[dist_leq_min_upper_bound] = 1. / distances[dist_leq_min_upper_bound] ** 2 + # stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 1, 5, 100/5)) + # stdout.flush() + # sums = np.sum(w[dist_leq_min_upper_bound], axis=1, keepdims=True) + # stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 2, 5, 2*100/5)) + # stdout.flush() + # weights[dist_leq_min_upper_bound] = w[dist_leq_min_upper_bound] / sums + # stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 3, 5, 3*100/5)) + # stdout.flush() + # wz = np.einsum('ij,ij->i', weights, z[indexes]) + # stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 4, 5, 4*100/5)) + # stdout.flush() + # idxs[dist_leq_min_upper_bound] = indexes[dist_leq_min_upper_bound] + # result[dist_leq_min_upper_bound] = wz[dist_leq_min_upper_bound] + # stdout.write('{}Building coeffs: {}/{} (100%)\n'.format(back_char, 5, 5)) + # stdout.flush() + # return result, weights, idxs + + # Inverse distance weights (IDW) interpolation, with optional Angular distance weighting (ADW) factor + # if adw_type == None -> simple IDW + # if adw_type == Shepard -> Shepard 1968 original formulation + def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None): z = self.z result = mask_it(np.empty((len(distances),) + np.shape(z[0])), self._mv_target, 1) weights = empty((len(distances),) + (nnear,)) @@ -446,12 +491,77 @@ def _build_weights_invdist(self, distances, indexes, nnear): w[dist_leq_min_upper_bound] = 1. / distances[dist_leq_min_upper_bound] ** 2 stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 1, 5, 100/5)) stdout.flush() + + if (adw_type=='Shepard'): + # initialize weights + weight_directional = np.zeros_like(distances) + + # get "s" vector as the inverse distance with power = 1 (in "w" we have the invdist power 2) + # actually s(d) should be + # 1/d when 0 < d < r'/3 + # (27/4r')*(d/r'-1) when r'/3 < d < r' + # 0 when r' < d + # The orginal method consider a range of 4 to 10 data points e removes distant point using the rule: + # pi*r= 7(N/A) where N id the number of points and A is the area of the largest poligon enclosed by data points + # r'(C^n) = di(n+1) where di is the Point having the minimum distance major than the first n Points + # so that + # r' = r' C^4 = di(5) when n(C) <= 4 + # r' = r when 4 < n(C) <= 10 + # r' = r' C^10 = di(11) when 10 < n(C) + # + # TODO: check if the implementation needs to be updated accordingly + # here we take only the inverse of the distance in "s" + s = np.zeros_like(distances) + s[dist_leq_min_upper_bound] = 1. / distances[dist_leq_min_upper_bound] + + self.lat_inALL = self.target_latsOR.ravel() + self.lon_inALL = self.target_lonsOR.ravel() + + for i in range(nnear): + xi = self.latgrib[indexes[:,i]] + yi = self.longrib[indexes[:,i]] + + xi_diff = self.lat_inALL - xi + yi_diff = self.lon_inALL - yi + + Di = np.sqrt(xi_diff**2 + yi_diff**2) + + xj_diff = self.lat_inALL[:, np.newaxis] - self.latgrib[indexes] + yj_diff = self.lon_inALL[:, np.newaxis] - self.longrib[indexes] + + Dj = np.sqrt(xj_diff**2 + yj_diff**2) + # TODO: check here: we could use "distances", but we are actually evaluating the cos funcion on lat and lon values, that + # approximates the real angle... to use "distances" we need to operate in 3d version of the points + # in theory it should be OK to use the solution above, otherwise we can change it to + # Di = distances[i] + # Dj = distances + # and formula cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj) + (z-zi)*(z-zj)] / di*dj + + # cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj)] / di*dj. + cos_theta = (xi_diff[:,np.newaxis] * xj_diff + yi_diff[:,np.newaxis] * yj_diff) / (Di[:,np.newaxis] * Dj) + + numerator = np.sum((1 - cos_theta) * s, axis=1) + denominator = np.sum(s, axis=1) + + weight_directional[:, i] = numerator / denominator + + # update weights with directional ones + w = np.multiply(w,1+weight_directional) + elif (adw_type=='CDD'): + raise ApplicationException.get_exc(INVALID_INTERPOL_METHOD, + f"interpolation method not implemented yet (mode = {self.mode}, nnear = {self.nnear}, adw_type = {adw_type})") + elif (adw_type!=None): + raise ApplicationException.get_exc(INVALID_INTERPOL_METHOD, + f"interpolation method not supported (mode = {self.mode}, nnear = {self.nnear}, adw_type = {adw_type})") + + #normalize weights sums = np.sum(w[dist_leq_min_upper_bound], axis=1, keepdims=True) stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 2, 5, 2*100/5)) stdout.flush() weights[dist_leq_min_upper_bound] = w[dist_leq_min_upper_bound] / sums stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 3, 5, 3*100/5)) stdout.flush() + wz = np.einsum('ij,ij->i', weights, z[indexes]) stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 4, 5, 4*100/5)) stdout.flush() @@ -459,6 +569,7 @@ def _build_weights_invdist(self, distances, indexes, nnear): result[dist_leq_min_upper_bound] = wz[dist_leq_min_upper_bound] stdout.write('{}Building coeffs: {}/{} (100%)\n'.format(back_char, 5, 5)) stdout.flush() + return result, weights, idxs # take additional points from the KDTree close to the current point and replace the wrong ones with a new ones @@ -744,7 +855,7 @@ def _build_weights_triangulation(self, use_bilinear = False): longrib_max = normalized_longrib.max() longrib_min = normalized_longrib.min() - # evaluate an approx_grib_resolution by using 10 times the first longidure values + # evaluate an approx_grib_resolution by using 10 times the first longitude values # to check if the whole globe is covered approx_grib_resolution = abs(normalized_longrib[0]-normalized_longrib[1])*1.5 is_global_map = (360-(longrib_max-longrib_min)) Date: Thu, 22 Jun 2023 13:36:09 +0200 Subject: [PATCH 12/17] Added further optimizations and ADW unit test --- .../interpolation/scipy_interpolation_lib.py | 73 ++++++++++++------- tests/test_interpolation.py | 20 +++++ 2 files changed, 65 insertions(+), 28 deletions(-) diff --git a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py index 2a01e1a..a8a8cf5 100644 --- a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py +++ b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py @@ -465,7 +465,7 @@ def _build_nn(self, distances, indexes): # Inverse distance weights (IDW) interpolation, with optional Angular distance weighting (ADW) factor # if adw_type == None -> simple IDW # if adw_type == Shepard -> Shepard 1968 original formulation - def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None): + def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use_broadcasting = False): z = self.z result = mask_it(np.empty((len(distances),) + np.shape(z[0])), self._mv_target, 1) weights = empty((len(distances),) + (nnear,)) @@ -517,34 +517,51 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None): self.lat_inALL = self.target_latsOR.ravel() self.lon_inALL = self.target_lonsOR.ravel() - for i in range(nnear): - xi = self.latgrib[indexes[:,i]] - yi = self.longrib[indexes[:,i]] - - xi_diff = self.lat_inALL - xi - yi_diff = self.lon_inALL - yi - + # start_time = time.time() + if not use_broadcasting: + for i in range(nnear): + xj_diff = self.lat_inALL[:, np.newaxis] - self.latgrib[indexes] + yj_diff = self.lon_inALL[:, np.newaxis] - self.longrib[indexes] + + Dj = np.sqrt(xj_diff**2 + yj_diff**2) + # TODO: check here: we could use "distances", but we are actually evaluating the cos funcion on lat and lon values, that + # approximates the real angle... to use "distances" we need to operate in 3d version of the points + # in theory it should be OK to use the solution above, otherwise we can change it to + # Di = distances[i] + # Dj = distances + # and formula cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj) + (z-zi)*(z-zj)] / di*dj + + # cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj)] / di*dj. + #cos_theta = (xi_diff[:,np.newaxis] * xj_diff + yi_diff[:,np.newaxis] * yj_diff) / (Di[:,np.newaxis] * Dj) + cos_theta = (xj_diff[:,i,np.newaxis] * xj_diff + yj_diff[:,i,np.newaxis] * yj_diff) / (Dj[:,i,np.newaxis] * Dj) + + numerator = np.sum((1 - cos_theta) * s, axis=1) + denominator = np.sum(s, axis=1) + + weight_directional[dist_leq_min_upper_bound, i] = numerator[dist_leq_min_upper_bound] / denominator[dist_leq_min_upper_bound] + else: + # All in broadcasting: + # this algorithm uses huge amount of memory and deas not speed up much on standard Virtual Machine + # TODO: check if it is worth to use it on HPC + # Compute xi and yi for all elements + xi = self.latgrib[indexes] + yi = self.longrib[indexes] + # Compute xi_diff and yi_diff for all elements + xi_diff = self.lat_inALL[:, np.newaxis] - xi + yi_diff = self.lon_inALL[:, np.newaxis] - yi + # Compute Di for all elements Di = np.sqrt(xi_diff**2 + yi_diff**2) - - xj_diff = self.lat_inALL[:, np.newaxis] - self.latgrib[indexes] - yj_diff = self.lon_inALL[:, np.newaxis] - self.longrib[indexes] - - Dj = np.sqrt(xj_diff**2 + yj_diff**2) - # TODO: check here: we could use "distances", but we are actually evaluating the cos funcion on lat and lon values, that - # approximates the real angle... to use "distances" we need to operate in 3d version of the points - # in theory it should be OK to use the solution above, otherwise we can change it to - # Di = distances[i] - # Dj = distances - # and formula cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj) + (z-zi)*(z-zj)] / di*dj - - # cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj)] / di*dj. - cos_theta = (xi_diff[:,np.newaxis] * xj_diff + yi_diff[:,np.newaxis] * yj_diff) / (Di[:,np.newaxis] * Dj) - - numerator = np.sum((1 - cos_theta) * s, axis=1) - denominator = np.sum(s, axis=1) - - weight_directional[:, i] = numerator / denominator - + # Compute cos_theta for all elements + cos_theta = (xi_diff[:, :, np.newaxis] * xi_diff[:, np.newaxis, :] + yi_diff[:, :, np.newaxis] * yi_diff[:, np.newaxis, :]) / (Di[:, :, np.newaxis] * Di[:, np.newaxis, :]) + numerator = np.sum((1 - cos_theta) * s[:, np.newaxis, :], axis=2) + denominator = np.sum(s[:, np.newaxis, :], axis=2) + # Compute weight_directional for all elements + weight_directional[dist_leq_min_upper_bound, :] = numerator[dist_leq_min_upper_bound, :] / denominator[dist_leq_min_upper_bound, :] + + # end_time = time.time() + # elapsed_time = end_time - start_time + # print(f"Elapsed time for computation: {elapsed_time:.6f} seconds") + # update weights with directional ones w = np.multiply(w,1+weight_directional) elif (adw_type=='CDD'): diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 5d990c8..e64c8d2 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -45,6 +45,26 @@ def test_interpolation_create_scipy_invdist(self): assert shape_target == values_resampled.shape os.unlink('tests/data/tbl_pf10tp_550800_scipy_invdist.npy.gz') + @pytest.mark.slow + def test_interpolation_create_scipy_adw(self): + d = deepcopy(config_dict) + d['interpolation.create'] = True + d['interpolation.parallel'] = True + d['interpolation.mode'] = 'adw' + file = d['input.file'] + reader = GRIBReader(file) + messages = reader.select_messages(shortName='2t') + grid_id = messages.grid_id + missing = messages.missing_value + ctx = MockedExecutionContext(d, False) + interpolator = Interpolator(ctx, missing) + values_in = messages.values_first_or_single_res[messages.first_step_range] + lats, lons = messages.latlons + values_resampled = interpolator.interpolate_scipy(lats, lons, values_in, grid_id, messages.grid_details) + shape_target = PCRasterReader(d['interpolation.latMap']).values.shape + assert shape_target == values_resampled.shape + os.unlink('tests/data/tbl_pf10tp_550800_scipy_adw.npy.gz') + @pytest.mark.slow def test_interpolation_create_scipy_bilinear(self): d = deepcopy(config_dict) From c797204d5abc9b7738f8904e4165ff0d73830414 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Fri, 7 Jul 2023 10:45:16 +0200 Subject: [PATCH 13/17] Modified ADW to exclude identity matrix points in angular weights evaluation. Added flag to enable full broadcasting in adw. Added some debugging. --- README.md | 14 +++ configuration/global/intertables.json | 11 +++ src/pyg2p/main/api.py | 1 + src/pyg2p/main/context.py | 1 + src/pyg2p/main/interpolation/__init__.py | 33 +++++-- .../interpolation/scipy_interpolation_lib.py | 86 ++++++++++++++++--- tests/test_interpolation.py | 1 + 7 files changed, 129 insertions(+), 18 deletions(-) diff --git a/README.md b/README.md index 2c55c35..8fe5c2c 100644 --- a/README.md +++ b/README.md @@ -669,6 +669,20 @@ Attributes p, leafsize and eps for the kd tree algorithm are default in scipy li | eps | 0 | | 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 + +```json +{ +"Interpolation": { + "@latMap": "/dataset/maps/europe5km/lat.map", + "@lonMap": "/dataset/maps/europe5km/long.map", + "@mode": "adw", + "@adw_broadcasting": False} +} +``` + #### bilinear It's the bilinear interpolation algorithm applyied on regular and irregular grids. On irregular grids, it tries to get the best quatrilateral around each target point, but at the same time tries to use the best stable and grid-like shape from starting points. To do so, evaluates interpolation looking at point on similar latitude, thus on projected grib files may show some irregular results. diff --git a/configuration/global/intertables.json b/configuration/global/intertables.json index e164b50..2fb617d 100644 --- a/configuration/global/intertables.json +++ b/configuration/global/intertables.json @@ -21,6 +21,17 @@ 680 ] }, + "I-15.75_16.125_511_415_212065_rotated_ll_-1700000_2700000_5000_-5000_34.97_71.07_-1700000_2700000_5000_-5000_-24.35_51.35scipy_adw": { + "filename": "tbl_pf10tp_550800_scipy_adw.npy", + "method": "adw", + "source_shape": [ + 212065 + ], + "target_shape": [ + 810, + 680 + ] + }, "I-15.75_16.125_511_415_212065_rotated_ll_-1700000_2700000_5000_-5000_34.97_71.07_-1700000_2700000_5000_-5000_-24.35_51.35scipy_nearest": { "filename": "tbl_pf10slhf_550800_scipy_nearest.npy", "method": "nearest", diff --git a/src/pyg2p/main/api.py b/src/pyg2p/main/api.py index 70b9341..5d4b3f8 100644 --- a/src/pyg2p/main/api.py +++ b/src/pyg2p/main/api.py @@ -168,6 +168,7 @@ 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.use_broadcasting'] = interpolation_conf.get('use_broadcasting', False) 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'] diff --git a/src/pyg2p/main/context.py b/src/pyg2p/main/context.py index 14681d3..92bc6e2 100644 --- a/src/pyg2p/main/context.py +++ b/src/pyg2p/main/context.py @@ -360,6 +360,7 @@ 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.rotated_target'] = interpolation_conf.get('@rotated_target', False) if not self._vars['interpolation.dir'] and interpolation_conf.get('@intertableDir'): # get from JSON diff --git a/src/pyg2p/main/interpolation/__init__.py b/src/pyg2p/main/interpolation/__init__.py index 4904e52..5759d80 100644 --- a/src/pyg2p/main/interpolation/__init__.py +++ b/src/pyg2p/main/interpolation/__init__.py @@ -9,7 +9,7 @@ from . import grib_interpolation_lib from .latlong import LatLong -from .scipy_interpolation_lib import ScipyInterpolation, DEBUG_BILINEAR_INTERPOLATION, \ +from .scipy_interpolation_lib import ScipyInterpolation, DEBUG_BILINEAR_INTERPOLATION, DEBUG_ADW_INTERPOLATION, \ DEBUG_MIN_LAT, DEBUG_MIN_LON, DEBUG_MAX_LAT, DEBUG_MAX_LON from ...exceptions import ApplicationException, NO_INTERTABLE_CREATED @@ -31,6 +31,7 @@ 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._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') @@ -177,11 +178,33 @@ def interpolate_scipy(self, latgrib, longrib, v, grid_id, grid_details=None): else: intertable_id, intertable_name = self._intertable_filename(grid_id) + if DEBUG_ADW_INTERPOLATION: + # to debug create a limited controlled set of input values + # + # np.random.seed(0) + # latgrib = np.random.uniform(low=3, high=11, size=10) + # longrib = np.random.uniform(low=46, high=50, size=10) + # v = np.random.uniform(low=100, high=200, size=10) + # latgrib = np.array([ 7.39050803, 8.72151493, 7.82210701, 7.35906546, 6.38923839, + # 8.1671529, 6.50069769, 10.13418401, 10.70930208, 6.06753215]) + # longrib = np.array([49.16690015, 48.11557968, 48.27217824, 49.70238655, 46.28414423, + # 46.3485172, 46.08087359, 49.33047938, 49.112627, 49.48004859]) + # v = np.array([197.86183422, 179.91585642, 146.14793623, 178.05291763, 111.82744259, + # 163.99210213, 114.33532874, 194.4668917, 152.18483218, 141.466194 ]) + # 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 ]) + intertable_id, intertable_name = 'DEBUG_ADW','DEBUG_ADW.npy' + nnear = self.scipy_modes_nnear[self._mode] - if pyg2p.util.files.exists(intertable_name) or pyg2p.util.files.exists(intertable_name + '.gz'): - indexes, weights = self._read_intertable(intertable_name) - result = self._interpolate_scipy_append_mv(v, weights, indexes, nnear) + if (not DEBUG_ADW_INTERPOLATION) and \ + (pyg2p.util.files.exists(intertable_name) or pyg2p.util.files.exists(intertable_name + '.gz')): + indexes, weights = self._read_intertable(intertable_name) + result = self._interpolate_scipy_append_mv(v, weights, indexes, nnear) elif self.create_if_missing: self.intertables_config.check_write() @@ -192,7 +215,7 @@ def interpolate_scipy(self, latgrib, longrib, v, grid_id, grid_details=None): 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) + parallel=self.parallel, mode=self._mode, use_broadcasting=self._adw_broadcasting) _, weights, indexes = scipy_interpolation.interpolate(lonefas, latefas) result = self._interpolate_scipy_append_mv(v, weights, indexes, nnear) diff --git a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py index a8a8cf5..46671af 100644 --- a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py +++ b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py @@ -19,6 +19,7 @@ DEBUG_BILINEAR_INTERPOLATION = False +DEBUG_ADW_INTERPOLATION = False # DEBUG_MIN_LAT = 88 # DEBUG_MIN_LON = -180 # DEBUG_MAX_LAT = 90 @@ -31,10 +32,14 @@ # DEBUG_MIN_LON = -24 # DEBUG_MAX_LAT = 70 # DEBUG_MAX_LON = -22 -DEBUG_MIN_LAT = -10 -DEBUG_MIN_LON = -100 -DEBUG_MAX_LAT = 25 -DEBUG_MAX_LON = -50 +# DEBUG_MIN_LAT = -10 +# DEBUG_MIN_LON = -100 +# DEBUG_MAX_LAT = 25 +# DEBUG_MAX_LON = -50 +DEBUG_MIN_LAT = 7 +DEBUG_MIN_LON = 45 +DEBUG_MAX_LAT = 9 +DEBUG_MAX_LON = 50 #DEBUG_NN = 15410182 @@ -283,7 +288,7 @@ class ScipyInterpolation(object): def __init__(self, longrib, latgrib, grid_details, source_values, nnear, mv_target, mv_source, target_is_rotated=False, parallel=False, - mode='nearest'): + mode='nearest', use_broadcasting = False): stdout.write('Start scipy interpolation: {}\n'.format(now_string())) self.geodetic_info = grid_details self.source_grid_is_rotated = 'rotated' in grid_details.get('gridType') @@ -294,6 +299,7 @@ def __init__(self, longrib, latgrib, grid_details, source_values, nnear, self.latgrib = latgrib self.nnear = nnear self.mode = mode + self.use_broadcasting = use_broadcasting self._mv_target = mv_target self._mv_source = mv_source self.z = source_values @@ -339,9 +345,9 @@ def interpolate(self, target_lons, target_lats): # return results, distances, indexes result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear) elif self.mode == 'adw': - result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='Shepard') + result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='Shepard', use_broadcasting=self.use_broadcasting) elif self.mode == 'cdd': - result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='CDD') + result, weights, indexes = self._build_weights_invdist(distances, indexes, self.nnear, adw_type='CDD', use_broadcasting=self.use_broadcasting) elif self.mode == 'bilinear' and self.nnear == 4: # bilinear interpolation only supported with nnear = 4 # BILINEAR INTERPOLATION result, weights, indexes = self._build_weights_bilinear(distances, indexes, efas_locations, self.nnear) @@ -466,6 +472,12 @@ def _build_nn(self, distances, indexes): # if adw_type == None -> simple IDW # if adw_type == Shepard -> Shepard 1968 original formulation def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use_broadcasting = False): + if DEBUG_ADW_INTERPOLATION: + self.min_upper_bound = 1000000000 + if DEBUG_BILINEAR_INTERPOLATION: + n_debug=1940 + else: + n_debug=11805340 z = self.z result = mask_it(np.empty((len(distances),) + np.shape(z[0])), self._mv_target, 1) weights = empty((len(distances),) + (nnear,)) @@ -533,12 +545,27 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use # cos_theta = [(x-xi)*(x-xj) + (y-yi)*(y-yj)] / di*dj. #cos_theta = (xi_diff[:,np.newaxis] * xj_diff + yi_diff[:,np.newaxis] * yj_diff) / (Di[:,np.newaxis] * Dj) - cos_theta = (xj_diff[:,i,np.newaxis] * xj_diff + yj_diff[:,i,np.newaxis] * yj_diff) / (Dj[:,i,np.newaxis] * Dj) - - numerator = np.sum((1 - cos_theta) * s, axis=1) - denominator = np.sum(s, axis=1) + # cos_theta = (xj_diff[:,i,np.newaxis] * xj_diff + yj_diff[:,i,np.newaxis] * yj_diff) / (Dj[:,i,np.newaxis] * Dj) + # numerator = np.sum((1 - cos_theta) * s, axis=1) + # denominator = np.sum(s, axis=1) + cos_theta = (xj_diff[:,i,np.newaxis] * xj_diff + yj_diff[:,i,np.newaxis] * yj_diff) / (Dj[:,i,np.newaxis] * Dj) + # skip when i==j, so that directional weight of i is evaluated on all points j where j!=i + # TODO: tip: since cos_theta = 1 for i==j, to speed up we can skip this passage and apply i!=j only on denominator + cos_theta = cos_theta[:,np.concatenate([np.arange(i), np.arange(i+1, cos_theta.shape[1])])] + sj = s[:,np.concatenate([np.arange(i), np.arange(i+1, s.shape[1])])] + numerator = np.sum((1 - cos_theta) * sj, axis=1) + denominator = np.sum(sj, axis=1) + weight_directional[dist_leq_min_upper_bound, i] = numerator[dist_leq_min_upper_bound] / denominator[dist_leq_min_upper_bound] + + # DEBUG: test the point at n_debug 11805340=lat 8.025 lon 47.0249999 + if DEBUG_ADW_INTERPOLATION: + print("i: {}".format(i)) + print("cos_theta: {}".format(cos_theta[n_debug])) + print("s: {}".format(s[n_debug])) + print("numerator: {}".format(numerator[n_debug])) + print("denominator: {}".format(denominator[n_debug])) else: # All in broadcasting: # this algorithm uses huge amount of memory and deas not speed up much on standard Virtual Machine @@ -553,8 +580,21 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use Di = np.sqrt(xi_diff**2 + yi_diff**2) # Compute cos_theta for all elements cos_theta = (xi_diff[:, :, np.newaxis] * xi_diff[:, np.newaxis, :] + yi_diff[:, :, np.newaxis] * yi_diff[:, np.newaxis, :]) / (Di[:, :, np.newaxis] * Di[:, np.newaxis, :]) - numerator = np.sum((1 - cos_theta) * s[:, np.newaxis, :], axis=2) - denominator = np.sum(s[:, np.newaxis, :], axis=2) + # skip when i==j, so that directional weight of i is evaluated on all points j where j!=i + # TODO: tip: since cos_theta = 1 for i==j, to speed up we can skip this passage and apply i!=j only on denominator + # Delete the diagonal elements from cos_theta + # Reshape cos_theta to (n, nnear, nnear-1) + n = cos_theta.shape[0] + m = cos_theta.shape[1] + strided = np.lib.stride_tricks.as_strided + s0,s1,s2 = cos_theta[:].strides + cos_theta = strided(cos_theta.ravel()[1:], shape=(n,m-1,m), strides=(s0,s1+s2,s2)).reshape(n,m,-1) + sj = np.tile(s[:, np.newaxis, :], (1, 4, 1)) + s0,s1,s2 = sj[:].strides + sj = strided(sj.ravel()[1:], shape=(n,m-1,m), strides=(s0,s1+s2,s2)).reshape(n,m,-1) + + numerator = np.sum((1 - cos_theta) * sj, axis=2) + denominator = np.sum(sj, axis=2) # Compute weight_directional for all elements weight_directional[dist_leq_min_upper_bound, :] = numerator[dist_leq_min_upper_bound, :] / denominator[dist_leq_min_upper_bound, :] @@ -563,7 +603,12 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use # print(f"Elapsed time for computation: {elapsed_time:.6f} seconds") # update weights with directional ones + if DEBUG_ADW_INTERPOLATION: + print("weight_directional: {}".format(weight_directional[n_debug])) + print("w (before adding angular weights): {}".format(w[n_debug])) w = np.multiply(w,1+weight_directional) + if DEBUG_ADW_INTERPOLATION: + print("w (after adding angular weights): {}".format(w[n_debug])) elif (adw_type=='CDD'): raise ApplicationException.get_exc(INVALID_INTERPOL_METHOD, f"interpolation method not implemented yet (mode = {self.mode}, nnear = {self.nnear}, adw_type = {adw_type})") @@ -576,6 +621,12 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 2, 5, 2*100/5)) stdout.flush() weights[dist_leq_min_upper_bound] = w[dist_leq_min_upper_bound] / sums + if DEBUG_ADW_INTERPOLATION: + dist_smalltest = distances[:, 0] <= 10000 + onlyfirst_array_test = np.zeros(nnear) + onlyfirst_array_test[0] = 1000 + weights[dist_smalltest]=onlyfirst_array_test + print("weights (after normalization): {}".format(weights[n_debug])) stdout.write('{}Building coeffs: {}/{} ({:.2f}%)\n'.format(back_char, 3, 5, 3*100/5)) stdout.flush() @@ -584,6 +635,15 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use stdout.flush() idxs[dist_leq_min_upper_bound] = indexes[dist_leq_min_upper_bound] result[dist_leq_min_upper_bound] = wz[dist_leq_min_upper_bound] + if DEBUG_ADW_INTERPOLATION: + print("result: {}".format(result[n_debug])) + print("lat: {}".format(self.lat_inALL[n_debug])) + print("lon: {}".format(self.lon_inALL[n_debug])) + print("idxs: {}".format(idxs[n_debug])) + print("distances: {}".format(distances[n_debug])) + print("latgrib: {}".format(self.latgrib[idxs[n_debug]])) + print("longrib: {}".format(self.longrib[idxs[n_debug]])) + print("value: {}".format(self.z[idxs[n_debug]])) stdout.write('{}Building coeffs: {}/{} (100%)\n'.format(back_char, 5, 5)) stdout.flush() diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index e64c8d2..67009b2 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -51,6 +51,7 @@ def test_interpolation_create_scipy_adw(self): d['interpolation.create'] = True d['interpolation.parallel'] = True d['interpolation.mode'] = 'adw' + d['interpolation.adw_broadcasting'] = True file = d['input.file'] reader = GRIBReader(file) messages = reader.select_messages(shortName='2t') From 384a2a790ffd1327471ce318979a4c77d1059b93 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 12 Jul 2023 17:14:31 +0200 Subject: [PATCH 14/17] Fixed key time shift with multiple resolution inputs in average computation --- src/pyg2p/main/manipulation/aggregator.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/pyg2p/main/manipulation/aggregator.py b/src/pyg2p/main/manipulation/aggregator.py index 6035e63..ea305a1 100644 --- a/src/pyg2p/main/manipulation/aggregator.py +++ b/src/pyg2p/main/manipulation/aggregator.py @@ -165,27 +165,20 @@ def _average(self, values): shape_iter = values[first_key].shape v_ord = collections.OrderedDict(sorted(dict((k.end_step, v_) for (k, v_) in values.items()).items(), key=lambda k: k)) - iter_end = self._end - self._aggregation_step + 2 - if self._start > 0 and not self._second_t_res: - if self._aggregation_halfweights: - iter_start = self._start - self._aggregation_step - iter_end = self._end - self._aggregation_step + 1 - else: - iter_start = self._start - self._aggregation_step + 1 - elif self._second_t_res: + if self._start == 0 or self._second_t_res: iter_start = self._start + iter_end = self._end - self._aggregation_step + 2 else: - iter_start = 0 + iter_start = self._start - self._aggregation_step + iter_end = self._end - self._aggregation_step + 1 for iter_ in range(iter_start, iter_end, self._aggregation_step): - if self._aggregation_halfweights == False and self._start == 0: - iter_from = iter_ + 1 + if self._aggregation_halfweights or self._start == 0: + iter_from = iter_ iter_to = iter_ + self._aggregation_step + 1 else: - iter_from = iter_ - iter_to = iter_ + self._aggregation_step - if self._aggregation_halfweights: - iter_to += 1 + iter_from = iter_ + 1 + iter_to = iter_ + self._aggregation_step + 1 temp_sum = np.zeros(shape_iter) v_ord_keys = list(v_ord.keys()) From 0cc593f7264deb47bbfdfb30bb8e4c99589f744b Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Wed, 12 Jul 2023 17:16:30 +0200 Subject: [PATCH 15/17] Updated version number --- src/pyg2p/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyg2p/VERSION b/src/pyg2p/VERSION index 448ada3..c4a602d 100644 --- a/src/pyg2p/VERSION +++ b/src/pyg2p/VERSION @@ -1 +1 @@ -3.2.5 \ No newline at end of file +3.2.6 \ No newline at end of file From 0e53ef05f6aa1d2feceeff5f9dac91a2c97c189b Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Mon, 17 Jul 2023 09:06:26 +0200 Subject: [PATCH 16/17] Fixed read of halfweights parameter from json file --- README.md | 12 ++++++------ src/pyg2p/VERSION | 2 +- src/pyg2p/main/context.py | 1 + 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 8fe5c2c..54fe4be 100644 --- a/README.md +++ b/README.md @@ -410,7 +410,7 @@ defining the step numbers to skip when writing PCRaster maps. Same as old grib2p grib2pcraster). It can be average or accumulation. -  halfweightsIf set to True and type is "average", it will avaluate the average by using half weights for the first and the last step +  halfweightsIf set to true and type is "average", it will avaluate the average by using half weights for the first and the last step  forceZeroArrayOptional. In case of “accumulation”, and only @@ -748,7 +748,7 @@ The JSON configuration in the execution file will look like: { "Aggregation": { "@type": "average", - "@halfweights": False} + "@halfweights": false} } ``` @@ -766,7 +766,7 @@ Temperatures are often extracted as averages on 24 hours or 6 hours. Here's a ty "Aggregation": { "@step": 24, "@type": "average", -"@halfweights": False +"@halfweights": false }, "OutMaps": { "@cloneMap": "/dataset/maps/europe/dem.map", @@ -812,7 +812,7 @@ This is needed because we performed 24 hours average over 6 hourly steps. The to evaluate the average, the following steps are executed: -- when "halfweights" is False, the results of the function is the sum of all the values from "start_step-aggregation_step+1" to end_step, taking for each step the value corresponding to the next available value in the grib file. E.g: +- when "halfweights" is false, the results of the function is the sum of all the values from "start_step-aggregation_step+1" to end_step, taking for each step the value corresponding to the next available value in the grib file. E.g: INPUT: start_step=24, end_step=, aggregation_step=24 GRIB File: contains data starting from step 0 to 48 every 6 hours: 0,6,12,18,24,30,.... @@ -821,7 +821,7 @@ GRIB File: contains data starting from step 0 to 48 every 6 hours: 0,6,12,18,24, Day 2: same as Day 1 starting from (24+24)-24+1=25... -- when "halfweights" is True, the results of the function is the sum of all the values from "start_step-aggregation_step" to end_step, taking for each step the value corresponding to the next available value in the grib file but using half of the weights for the first and the last step in each aggregation_step cicle. E.g: +- when "halfweights" is true, the results of the function is the sum of all the values from "start_step-aggregation_step" to end_step, taking for each step the value corresponding to the next available value in the grib file but using half of the weights for the first and the last step in each aggregation_step cicle. E.g: INPUT: start_step=24, end_step=, aggregation_step=24 GRIB File: contains data starting from step 0 to 72 every 6 hours: 0,6,12,18,24,30,36,.... @@ -830,7 +830,7 @@ GRIB File: contains data starting from step 0 to 72 every 6 hours: 0,6,12,18,24, Day 2: same as Day 1 starting from (24+24)-24=24: the step 24 will have a weight of 3, while steps 30,36 and 42 will be counted 6 times, and finally the step 48 will have a weight of 3. -- if start_step is zero or is not specified, the aggregation will start from 1 +- if start_step is zero or is not specified, the aggregation will start from 0 ### Accumulation For precipitation values, accumulation over 6 or 24 hours is often performed. Here's an example of configuration and execution output in DEBUG mode. diff --git a/src/pyg2p/VERSION b/src/pyg2p/VERSION index c4a602d..448ada3 100644 --- a/src/pyg2p/VERSION +++ b/src/pyg2p/VERSION @@ -1 +1 @@ -3.2.6 \ No newline at end of file +3.2.5 \ No newline at end of file diff --git a/src/pyg2p/main/context.py b/src/pyg2p/main/context.py index 92bc6e2..90416b3 100644 --- a/src/pyg2p/main/context.py +++ b/src/pyg2p/main/context.py @@ -413,6 +413,7 @@ def _define_exec_params(self): if exec_conf.get('Aggregation'): self._vars['aggregation.step'] = exec_conf['Aggregation'].get('@step') self._vars['aggregation.type'] = exec_conf['Aggregation'].get('@type') + self._vars['aggregation.halfweights'] = bool(exec_conf['Aggregation'].get('@halfweights')) self._vars['execution.doAggregation'] = bool(self._vars.get('aggregation.step')) \ and bool(self._vars.get('aggregation.type')) From 01495b312efe56966875d117ca8a4935247d271a Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Mon, 17 Jul 2023 10:51:13 +0200 Subject: [PATCH 17/17] Fixed typo in readme file --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 54fe4be..502a174 100644 --- a/README.md +++ b/README.md @@ -671,7 +671,7 @@ Attributes p, leafsize and eps for the kd tree algorithm are default in scipy li #### 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 +If @adw_broadcasting is set to true, computations will run in full broadcasting mode but requires more memory ```json { @@ -679,7 +679,7 @@ If @adw_broadcasting is set to True, computations will run in full broadcasting "@latMap": "/dataset/maps/europe5km/lat.map", "@lonMap": "/dataset/maps/europe5km/long.map", "@mode": "adw", - "@adw_broadcasting": False} + "@adw_broadcasting": false} } ```