From 3862747773386aafe1b100fbe03bce7c67511e81 Mon Sep 17 00:00:00 2001 From: Carlo Russo Date: Thu, 19 Oct 2023 17:52:38 +0200 Subject: [PATCH] Bugfix in ADW Shepard (zero distance) --- .../main/interpolation/scipy_interpolation_lib.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py index 2df4263..27b3ba8 100644 --- a/src/pyg2p/main/interpolation/scipy_interpolation_lib.py +++ b/src/pyg2p/main/interpolation/scipy_interpolation_lib.py @@ -981,6 +981,18 @@ def _build_weights_invdist(self, distances, indexes, nnear, adw_type = None, use if (adw_type=='Shepard') or (adw_type=='CDD'): idxs = indexes result = wz + + if (adw_type=='Shepard'): + # in Shepard, we still have NaN values because of the denominator when the distance is equal to 0 + # so we just take the exact value when the point is exacly on the station coordinates + dist_leq_1e_10 = distances[:, 0] <= 1e-10 + + # 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]] else: idxs[dist_leq_min_upper_bound] = indexes[dist_leq_min_upper_bound] result[dist_leq_min_upper_bound] = wz[dist_leq_min_upper_bound]