From 37c5721f53077f426c15a49f4c113d0b92a8ef19 Mon Sep 17 00:00:00 2001 From: ctena Date: Mon, 12 Dec 2022 16:09:11 +0100 Subject: [PATCH 1/3] Fixed bug on horizontal interpolation when distance is zero. --- nes/interpolation/horizontal_interpolation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nes/interpolation/horizontal_interpolation.py b/nes/interpolation/horizontal_interpolation.py index 972c990..7b7ba65 100644 --- a/nes/interpolation/horizontal_interpolation.py +++ b/nes/interpolation/horizontal_interpolation.py @@ -413,6 +413,7 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): weight_matrix.hours_end = 0 weight_matrix.set_communicator(MPI.COMM_SELF) + dists[dists < 1] = 1 # take the reciprocals of the nearest neighbours distances inverse_dists = np.reciprocal(dists) inverse_dists_transf = inverse_dists.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1])) -- GitLab From f3546b2e3f5302a6f305ff630cbda0a49396a1f2 Mon Sep 17 00:00:00 2001 From: ctena Date: Thu, 15 Dec 2022 17:37:42 +0100 Subject: [PATCH 2/3] Horiz Interp: Changed the method to convert lat_lon to cartesian. Now providentia interpolation and NES uses the same one. Previous method provides slight differences between machines (Nord3v2 vs MN4) and that one produces exactly the same values. --- nes/interpolation/horizontal_interpolation.py | 42 +++++++++++++++++-- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/nes/interpolation/horizontal_interpolation.py b/nes/interpolation/horizontal_interpolation.py index 7b7ba65..7166812 100644 --- a/nes/interpolation/horizontal_interpolation.py +++ b/nes/interpolation/horizontal_interpolation.py @@ -10,6 +10,7 @@ from scipy import spatial from filelock import FileLock from datetime import datetime from warnings import warn +import pyproj def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, @@ -390,8 +391,11 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): # to cartesian ECEF (Earth Centred, Earth Fixed) coordinates, before # calculating distances. - src_mod_xy = lon_lat_to_cartesian(src_lon.flatten(), src_lat.flatten()) - dst_mod_xy = lon_lat_to_cartesian(dst_lon.flatten(), dst_lat.flatten()) + # src_mod_xy = lon_lat_to_cartesian(src_lon.flatten(), src_lat.flatten()) + # dst_mod_xy = lon_lat_to_cartesian(dst_lon.flatten(), dst_lat.flatten()) + + src_mod_xy = lon_lat_to_cartesian_ecef(src_lon.flatten(), src_lat.flatten()) + dst_mod_xy = lon_lat_to_cartesian_ecef(dst_lon.flatten(), dst_lat.flatten()) # generate KDtree using model 1 coordinates (i.e. the model grid you are # interpolating from) @@ -413,9 +417,10 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): weight_matrix.hours_end = 0 weight_matrix.set_communicator(MPI.COMM_SELF) - dists[dists < 1] = 1 # take the reciprocals of the nearest neighbours distances + dists[dists < 1] = 1 inverse_dists = np.reciprocal(dists) + inverse_dists_transf = inverse_dists.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1])) weight_matrix.variables['inverse_dists'] = {'data': inverse_dists_transf, 'units': 'm'} idx_transf = idx.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1])) @@ -426,10 +431,12 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): return weight_matrix -def lon_lat_to_cartesian(lon, lat, radius=1.0): +def lon_lat_to_cartesian(lon, lat, radius=6378137.0): """ Calculate lon, lat coordinates of a point on a sphere. + DEPRECATED!!!! + Parameters ---------- lon : np.array @@ -448,3 +455,30 @@ def lon_lat_to_cartesian(lon, lat, radius=1.0): z = radius * np.sin(lat_r) return np.column_stack([x, y, z]) + + +def lon_lat_to_cartesian_ecef(lon, lat): + """ + # convert observational/model geographic longitude/latitude coordinates to cartesian ECEF (Earth Centred, + # Earth Fixed) coordinates, assuming WGS84 datum and ellipsoid, and that all heights = 0. + # ECEF coordiantes represent positions (in meters) as X, Y, Z coordinates, approximating the earth surface + # as an ellipsoid of revolution. + # This conversion is for the subsequent calculation of euclidean distances of the model gridcell centres + # from each observational station. + # Defining the distance between two points on the earth's surface as simply the euclidean distance + # between the two lat/lon pairs could lead to inaccurate results depending on the distance + # between two points (i.e. 1 deg. of longitude varies with latitude). + + Parameters + ---------- + lon : np.array + Longitude values. + lat : np.array + Latitude values. + """ + lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') + ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') + + x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False) + + return np.column_stack([x, y, z]) -- GitLab From 0a95c1768a3662e00d95c70475666787e4203a6c Mon Sep 17 00:00:00 2001 From: ctena Date: Fri, 16 Dec 2022 15:29:39 +0100 Subject: [PATCH 3/3] Corrected the way to use level.positive metadata. added functionalities to play with that parameter --- nes/interpolation/vertical_interpolation.py | 9 +++++++-- nes/nc_projections/default_nes.py | 22 ++++++++++++++++++++- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/nes/interpolation/vertical_interpolation.py b/nes/interpolation/vertical_interpolation.py index 88e2a86..88f8fb5 100644 --- a/nes/interpolation/vertical_interpolation.py +++ b/nes/interpolation/vertical_interpolation.py @@ -173,11 +173,16 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', # Update level information new_lev_info = {'data': new_levels} + if 'positive' in self._lev.keys(): + # Vertical level direction + if flip: + self.reverse_level_direction() + new_lev_info['positive'] = self._lev['positive'] for var_attr, attr_info in self.variables[self.vertical_var_name].items(): if var_attr not in ['data', 'dimensions', 'crs', 'grid_mapping']: new_lev_info[var_attr] = copy(attr_info) - self.lev = new_lev_info - self._lev = new_lev_info + self.set_levels(new_lev_info) + self.free_vars(self.vertical_var_name) self.vertical_var_name = None diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 33d86b1..bfd93bd 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -374,6 +374,24 @@ class Nes(object): def get_full_levels(self): return self._lev + def set_level_direction(self, new_direction): + if new_direction not in ['up', 'down']: + raise ValueError("Level direction mus be up or down. '{0}' is not a valid option".format(new_direction)) + self._lev['positive'] = new_direction + self.lev['positive'] = new_direction + + return True + + def reverse_level_direction(self): + if 'positive' in self._lev.keys(): + if self._lev['positive'] == 'up': + self._lev['positive'] = 'down' + self.lev['positive'] = 'down' + else: + self._lev['positive'] = 'up' + self.lev['positive'] = 'up' + return True + def clear_communicator(self): """ Erase the communicator and the parallelization indexes. @@ -1938,7 +1956,9 @@ class Nes(object): lev.units = Units(self._lev['units'], formatted=True).units else: lev.units = '' - lev.positive = 'up' + if 'positive' in self._lev.keys(): + lev.positive = self._lev['positive'] + if self.size > 1: lev.set_collective(True) lev[:] = self._lev['data'] -- GitLab