diff --git a/nes/interpolation/horizontal_interpolation.py b/nes/interpolation/horizontal_interpolation.py index 972c990a949003fa9498780c2259aa50510f63e5..716681220df78b5b679f434e5a69a6fccc398557 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) @@ -414,7 +418,9 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): weight_matrix.set_communicator(MPI.COMM_SELF) # 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])) @@ -425,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 @@ -447,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]) diff --git a/nes/interpolation/vertical_interpolation.py b/nes/interpolation/vertical_interpolation.py index 88e2a86fc5597e8f5c4448bc8a8510b7a24b2d5a..88f8fb524c7069dd33c7f0194f878b5a18f541d6 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 33d86b139cfd3f29bc80c7f26bf60fb20f8584fe..bfd93bd2218cdc62809abecf0690e3f37d487f87 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']