From 5d37204779874ab86925a5165f67e49ecdd8ba1d Mon Sep 17 00:00:00 2001 From: ctena Date: Tue, 24 Jan 2023 16:48:07 +0100 Subject: [PATCH 1/5] Conservative horizontal interpolation: working in serial --- nes/methods/horizontal_interpolation.py | 202 +++++++++++++++++++++--- tests/sbatch_test_MN4.cmd | 35 ++++ tests/test_interp_conservative.py | 151 ++++++++++++++++++ 3 files changed, 370 insertions(+), 18 deletions(-) create mode 100644 tests/sbatch_test_MN4.cmd create mode 100644 tests/test_interp_conservative.py diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py index b5d5d18..f658c28 100644 --- a/nes/methods/horizontal_interpolation.py +++ b/nes/methods/horizontal_interpolation.py @@ -3,6 +3,7 @@ import sys import warnings import numpy as np +import pandas as pd import os import nes from mpi4py import MPI @@ -10,8 +11,13 @@ from scipy import spatial from filelock import FileLock from datetime import datetime from warnings import warn +import copy import pyproj +# CONSTANTS +NEAREST_OPTS = ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN'] +CONSERVATIVE_OPTS = ['Conservative', 'Area_Conservative', 'cons', 'conservative', 'area'] + def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, info=False, to_providentia=False): @@ -35,7 +41,8 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares to_providentia : bool Indicates if we want the interpolated grid in Providentia format. """ - + if info and self.master: + print("Creating Weight Matrix") # Obtain weight matrix if self.parallel_method == 'T': weights, idx = get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours) @@ -44,6 +51,13 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares else: raise NotImplemented("Parallel method {0} is not implemented yet for horizontal interpolations. Use 'T'".format( self.parallel_method)) + if info and self.master: + print("Weight Matrix done!") + + # idx[idx < 0] = np.nan + idx = np.ma.masked_array(idx, mask=idx == -999) + # weights[weights < 0] = np.nan + weights = np.ma.masked_array(weights, mask=weights == -999) # Copy NES final_dst = dst_grid.copy() @@ -62,6 +76,8 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares final_dst.hours_start = self.hours_start final_dst.hours_end = self.hours_end + if info and self.master: + print("Applying weights") # Apply weights for var_name, var_info in self.variables.items(): if info and self.master: @@ -95,6 +111,9 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares final_dst.global_attrs = self.global_attrs + if info and self.master: + print("Formatting") + if to_providentia: # self = experiment to interpolate (regular, rotated, etc.) # final_dst = interpolated experiment (points) @@ -192,13 +211,17 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour if len(weight_matrix.lev['data']) != n_neighbours: warn("The selected weight matrix does not have the same number of nearest neighbours." + "Re-calculating again but not saving it.") - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: raise NotImplementedError(kind) else: - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: raise NotImplementedError(kind) if weight_matrix_path is not None: @@ -206,19 +229,30 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour else: weight_matrix = None else: - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: + if kind in NEAREST_OPTS: if self.master: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) else: weight_matrix = None + elif kind in CONSERVATIVE_OPTS: + if self.master: + if self.master: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) + else: + weight_matrix = None + else: + weight_matrix = None else: raise NotImplementedError(kind) - # Normalize to 1 if self.master: - weights = np.array(np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64) / - np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64).sum(axis=1), - dtype=np.float64) + if kind in NEAREST_OPTS: + # Normalize to 1 + weights = np.array(np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) / + np.array(weight_matrix.variables['weight']['data'], dtype=np.float64).sum(axis=1), + dtype=np.float64) + else: + weights = np.array(weight_matrix.variables['weight']['data'][0], dtype=np.float64) idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int) else: weights = None @@ -267,13 +301,17 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou if len(weight_matrix.lev['data']) != n_neighbours: warn("The selected weight matrix does not have the same number of nearest neighbours." + "Re-calculating again but not saving it.") - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: raise NotImplementedError(kind) else: - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: raise NotImplementedError(kind) if weight_matrix_path is not None: @@ -281,19 +319,24 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou else: weight_matrix = None else: - if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']: - if self.master: + if self.master: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: - weight_matrix = None + raise NotImplementedError(kind) else: - raise NotImplementedError(kind) + weight_matrix = None # Normalize to 1 if self.master: - weights = np.array(np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64) / - np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64).sum(axis=1), - dtype=np.float64) + if kind in NEAREST_OPTS: + weights = np.array(np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) / + np.array(weight_matrix.variables['weight']['data'], dtype=np.float64).sum(axis=1), + dtype=np.float64) + else: + weights = np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int) else: weights = None @@ -339,6 +382,10 @@ def read_weight_matrix(weight_matrix_path, comm=None, parallel_method='T'): weight_matrix = nes.open_netcdf(path=weight_matrix_path, comm=comm, parallel_method=parallel_method, balanced=True) weight_matrix.load() + # In previous versions of NES weight was called inverse_dists + if 'inverse_dists' in weight_matrix.variables.keys(): + weight_matrix.variables['weight'] = weight_matrix.variables['inverse_dists'] + return weight_matrix @@ -412,6 +459,8 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): weight_matrix = dst_grid.copy() weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] + weight_matrix._time_bnds = None + weight_matrix.time_bnds = None weight_matrix.last_level = None weight_matrix.first_level = 0 weight_matrix.hours_start = 0 @@ -423,7 +472,7 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): 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'} + weight_matrix.variables['weight'] = {'data': inverse_dists_transf, 'units': 'm'} idx_transf = idx.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1])) weight_matrix.variables['idx'] = {'data': idx_transf, 'units': ''} weight_matrix.lev = {'data': np.arange(inverse_dists_transf.shape[1]), 'units': ''} @@ -432,6 +481,123 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): return weight_matrix +def create_area_conservative_weight_matrix(self, dst_nes, info=False): + """ + To create the weight matrix with the area conservative method. + + Parameters + ---------- + self : nes.Nes + Source projection Nes Object. + dst_nes : nes.Nes + Final projection Nes object. + + info: bool + Indicates if you want to print extra info during the methods process. + + Returns + ------- + nes.Nes + Weight matrix. + """ + + if info and self.master: + print("\tCreating area conservative Weight Matrix") + sys.stdout.flush() + # Get a portion of the destiny grid + if dst_nes.shapefile is None: + dst_nes.create_shapefile() + dst_grid = copy.deepcopy(dst_nes.shapefile) + + # Get the complete source grid + if self.shapefile is None: + self.create_shapefile() + src_grid = copy.deepcopy(self.shapefile) + if self.parallel_method == 'T': + # All process has the same shapefile + pass + else: + src_grid = self.comm.gather(src_grid, root=0) + if self.master: + src_grid = pd.concat(src_grid) + src_grid = self.comm.bcast(src_grid) + + # Normalizing projections + dst_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) + dst_grid.index.name = 'FID_dst' + src_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) + src_grid.index.name = 'FID_src' + + # Get intersected areas + inp, res = dst_grid.sindex.query_bulk(src_grid.geometry, predicate='intersects') + + # Calculate intersected areas and fractions + intersection = pd.concat([src_grid.geometry[inp].reset_index(), dst_grid.geometry[res].reset_index()], axis=1, + ignore_index=True) + intersection.columns = (list(src_grid.geometry[inp].reset_index().rename(columns={'geometry': 'geometry_src'}).columns) + list(dst_grid.geometry[res].reset_index().rename(columns={'geometry': 'geometry_dst'}).columns)) + intersection["src_area"] = src_grid.loc[intersection['FID_src'], 'geometry'].area.values + + intersection['weight'] = intersection.apply( + lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area / x['src_area'], axis=1) + + # Format & Clean + intersection.drop(columns=["geometry_src", "geometry_dst", "src_area"], inplace=True) + intersection = intersection.set_index(['FID_dst', intersection.groupby('FID_dst').cumcount()]).rename_axis( + ('FID', 'level')).sort_index() + intersection.rename(columns={"FID_src": "idx"}, inplace=True) + + # Initialising weight matrix + weight_matrix = dst_nes.copy() + weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] + weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] + weight_matrix._time_bnds = None + weight_matrix.time_bnds = None + weight_matrix.last_level = None + weight_matrix.first_level = 0 + weight_matrix.hours_start = 0 + weight_matrix.hours_end = 0 + + weight_matrix.set_communicator(MPI.COMM_SELF) + + weight_matrix.set_levels({'data': np.arange(intersection.index.get_level_values('level').max() + 1), + 'dimensions': ('lev',), + 'units': '', + 'positive': 'up'}) + + # Creating Weight matrix empty variables + if len(weight_matrix.lat['data'].shape) == 1: + shape = (1, len(weight_matrix.lev['data']), + weight_matrix.lat['data'].shape[0], weight_matrix.lon['data'].shape[0],) + shape_flat = (1, len(weight_matrix.lev['data']), + weight_matrix.lat['data'].shape[0] * weight_matrix.lon['data'].shape[0],) + else: + shape = (1, len(weight_matrix.lev['data']), + weight_matrix.lat['data'].shape[0], weight_matrix.lat['data'].shape[1],) + shape_flat = (1, len(weight_matrix.lev['data']), + weight_matrix.lat['data'].shape[0] * weight_matrix.lat['data'].shape[1],) + + weight_matrix.variables['weight'] = {'data': np.empty(shape_flat), 'units': '-'} + weight_matrix.variables['weight']['data'][:] = -999 + weight_matrix.variables['idx'] = {'data': np.empty(shape_flat), 'units': '-'} + weight_matrix.variables['idx']['data'][:] = -999 + + # Filling Weight matrix variables + for aux_lev in weight_matrix.lev['data']: + aux_data = intersection.xs(level='level', key=aux_lev) + weight_matrix.variables['weight']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'weight'] + weight_matrix.variables['idx']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'idx'] + + # Re-shaping + weight_matrix.variables['weight']['data'] = weight_matrix.variables['weight']['data'].reshape(shape) + weight_matrix.variables['idx']['data'] = weight_matrix.variables['idx']['data'].reshape(shape) + + # TODO: Transpose do to FID od shapefile are transposed. + for aux_lev in weight_matrix.lev['data']: + weight_matrix.variables['weight']['data'][0, aux_lev] = weight_matrix.variables['weight']['data'][0, aux_lev].T + weight_matrix.variables['idx']['data'][0, aux_lev] = weight_matrix.variables['idx']['data'][0, aux_lev].T + return weight_matrix + + def lon_lat_to_cartesian(lon, lat, radius=6378137.0): """ Calculate lon, lat coordinates of a point on a sphere. diff --git a/tests/sbatch_test_MN4.cmd b/tests/sbatch_test_MN4.cmd new file mode 100644 index 0000000..0901115 --- /dev/null +++ b/tests/sbatch_test_MN4.cmd @@ -0,0 +1,35 @@ +#!/bin/bash + +#SBATCH --qos=debug +#SBATCH -A bsc32 +#SBATCH --cpus-per-task=1 +#SBATCH -n 1 +#SBATCH -t 00:30:00 +#SBATCH -J NES_mn4_tests +#SBATCH --output=log_mn4_NES_%j.out +#SBATCH --error=log_mn4_NES_%j.err +#SBATCH --exclusive + +### ulimit -s 128000 + +module purge +module use /gpfs/projects/bsc32/software/suselinux/11/modules/all + +module load Python/3.7.4-GCCcore-8.3.0 +module load netcdf4-python/1.5.3-foss-2019b-Python-3.7.4 +module load cfunits/1.8-foss-2019b-Python-3.7.4 +module load xarray/0.17.0-foss-2019b-Python-3.7.4 +module load OpenMPI/4.0.5-GCC-8.3.0-mn4 +module load filelock/3.7.1-foss-2019b-Python-3.7.4 +module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4 +module load pyproj/2.5.0-foss-2019b-Python-3.7.4 +module load geopandas/0.8.1-foss-2019b-Python-3.7.4 +module load Shapely/1.7.1-foss-2019b-Python-3.7.4 + +# export PYTHONPATH=/gpfs/projects/bsc32/models/NES_master:${PYTHONPATH} +# cd /gpfs/projects/bsc32/models/NES_master/tests + +export PYTHONPATH=/gpfs/scratch/bsc32/bsc32538/NES_tests/NES:${PYTHONPATH} +cd /gpfs/scratch/bsc32/bsc32538/NES_tests/NES/tests + +mpirun --mca mpi_warn_on_fork 0 -np 1 python test_interp_conservative.py diff --git a/tests/test_interp_conservative.py b/tests/test_interp_conservative.py new file mode 100644 index 0000000..3514000 --- /dev/null +++ b/tests/test_interp_conservative.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python +import sys +import timeit +import pandas as pd +from mpi4py import MPI +from nes import * +import os, sys + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +parallel_method = 'Y' + +result_path = "Times_test_interp_conservative_{0}_{1:03d}.csv".format(parallel_method, size) +result = pd.DataFrame(index=['read', 'calcul', 'write'], + columns=['Only interp', 'Create_WM', "Read_WM"]) + +# NAMEE +src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc" +var_list = ['O3'] + +# ====================================================================================================================== +# ====================================== Only interp ===================================================== +# ====================================================================================================================== +test_name = "Only interp" +print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative') +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +# ====================================================================================================================== +# ====================================== Create_WM ===================================================== +# ====================================================================================================================== +test_name = "Create_WM" +print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +# Cleaning WM +if os.path.exists("NAMEE_to_IP.nc"): + os.remove("NAMEE_to_IP.nc") + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="NAMEE_to_IP.nc") +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +# ====================================================================================================================== +# ====================================== Read_WM ===================================================== +# ====================================================================================================================== +test_name = "Read_WM" +print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 + +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="NAMEE_to_IP.nc") +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +result.to_csv(result_path) \ No newline at end of file -- GitLab From 69b2a3e9e35ba39418117d033a4b690287d579a7 Mon Sep 17 00:00:00 2001 From: ctena Date: Fri, 27 Jan 2023 13:38:04 +0100 Subject: [PATCH 2/5] Time: Fixed bug when loading single time step NetCDF With NaN data inside --- nes/nc_projections/default_nes.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 5e62696..da964bf 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -1340,7 +1340,7 @@ class Nes(object): CF compliant time. """ - units = time.units + units = self.__parse_time_unit(time.units) if not hasattr(time, 'calendar'): calendar = 'standard' else: @@ -1350,7 +1350,12 @@ class Nes(object): units = 'days since ' + time.units.split('since')[1].lstrip() time = self.__get_dates_from_months(time, units, calendar) - return time, units, calendar + time_data = time[:] + + if len(time_data) == 1 and np.isnan(time_data[0]): + time_data[0] = 0 + + return time_data, units, calendar @staticmethod def __parse_time_unit(t_units): @@ -1389,8 +1394,9 @@ class Nes(object): else: if self.master: nc_var = self.netcdf.variables['time'] - nc_var, units, calendar = self.__parse_time(nc_var) - time = num2date(nc_var[:], self.__parse_time_unit(units), calendar=calendar) + time_data, units, calendar = self.__parse_time(nc_var) + + time = num2date(time_data, units, calendar=calendar) time = [aux.replace(second=0, microsecond=0) for aux in time] else: time = None -- GitLab From d7bbed7dfa22d74b93a99498d9a1a89ac89838f0 Mon Sep 17 00:00:00 2001 From: ctena Date: Thu, 2 Feb 2023 15:02:47 +0100 Subject: [PATCH 3/5] Working on conservative interpolation --- nes/methods/horizontal_interpolation.py | 82 ++++++++++++------- nes/nc_projections/default_nes.py | 10 ++- ... => 5.3-test_horiz_interp_conservative.py} | 20 +++-- tests/test_bash_mn4.cmd | 3 +- tests/test_bash_nord3v2.cmd | 1 + 5 files changed, 75 insertions(+), 41 deletions(-) rename tests/{test_interp_conservative.py => 5.3-test_horiz_interp_conservative.py} (90%) diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py index f658c28..8425bbb 100644 --- a/nes/methods/horizontal_interpolation.py +++ b/nes/methods/horizontal_interpolation.py @@ -20,7 +20,7 @@ CONSERVATIVE_OPTS = ['Conservative', 'Area_Conservative', 'cons', 'conservative' def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, - info=False, to_providentia=False): + info=False, to_providentia=False, only_create_wm=False): """ Horizontal methods from one grid to another one. @@ -33,13 +33,15 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares weight_matrix_path : str, None Path to the weight matrix to read/create. kind : str - Kind of horizontal methods. Accepted values: ['NearestNeighbour']. + Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative']. n_neighbours : int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. info : bool Indicates if you want to print extra info during the methods process. to_providentia : bool Indicates if we want the interpolated grid in Providentia format. + only_create_wm : bool + Indicates if you want to only create the Weight Matrix. """ if info and self.master: print("Creating Weight Matrix") @@ -53,11 +55,17 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares self.parallel_method)) if info and self.master: print("Weight Matrix done!") + if only_create_wm: + return None # idx[idx < 0] = np.nan idx = np.ma.masked_array(idx, mask=idx == -999) + # idx = np.array(idx, dtype=float) + # idx[idx < 0] = np.nan # weights[weights < 0] = np.nan weights = np.ma.masked_array(weights, mask=weights == -999) + # weights = np.array(weights, dtype=float) + # weights[weights < 0] = np.nan # Copy NES final_dst = dst_grid.copy() @@ -174,7 +182,8 @@ def get_src_data(comm, var_data, idx, parallel_method): var_data = comm.bcast(var_data) - var_data = np.take(var_data, idx) + var_data = np.pad(var_data, [1, 1], 'constant', constant_values=np.nan).take(idx + 1, mode='clip') + #var_data = np.take(var_data, idx) return var_data @@ -193,7 +202,7 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour weight_matrix_path : str, None Path to the weight matrix to read/create. kind : str - Kind of horizontal methods. Accepted values: ['NearestNeighbour']. + Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative']. n_neighbours : int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. @@ -208,15 +217,12 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour if self.master: if os.path.isfile(weight_matrix_path): weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF) - if len(weight_matrix.lev['data']) != n_neighbours: - warn("The selected weight matrix does not have the same number of nearest neighbours." + - "Re-calculating again but not saving it.") - if kind in NEAREST_OPTS: + if kind in NEAREST_OPTS: + if len(weight_matrix.lev['data']) != n_neighbours: + warn("The selected weight matrix does not have the same number of nearest neighbours." + + "Re-calculating again but not saving it.") weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) - elif kind in CONSERVATIVE_OPTS: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) - else: - raise NotImplementedError(kind) + else: if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) @@ -228,6 +234,9 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour weight_matrix.to_netcdf(weight_matrix_path) else: weight_matrix = None + if self.master: + if os.path.exists(weight_matrix_path.replace('.nc', '.lock')): + os.remove(weight_matrix_path.replace('.nc', '.lock')) else: if kind in NEAREST_OPTS: if self.master: @@ -278,7 +287,7 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou weight_matrix_path : str, None Path to the weight matrix to read/create. kind : str - Kind of horizontal methods. Accepted values: ['NearestNeighbour'] + Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative']. n_neighbours : int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. @@ -298,15 +307,11 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou if self.master: if os.path.isfile(weight_matrix_path): weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF) - if len(weight_matrix.lev['data']) != n_neighbours: - warn("The selected weight matrix does not have the same number of nearest neighbours." + - "Re-calculating again but not saving it.") - if kind in NEAREST_OPTS: + if kind in NEAREST_OPTS: + if len(weight_matrix.lev['data']) != n_neighbours: + warn("The selected weight matrix does not have the same number of nearest neighbours." + + "Re-calculating again but not saving it.") weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) - elif kind in CONSERVATIVE_OPTS: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) - else: - raise NotImplementedError(kind) else: if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) @@ -318,6 +323,9 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou weight_matrix.to_netcdf(weight_matrix_path) else: weight_matrix = None + if self.master: + if os.path.exists(weight_matrix_path.replace('.nc', '.lock')): + os.remove(weight_matrix_path.replace('.nc', '.lock')) else: if self.master: if kind in NEAREST_OPTS: @@ -522,10 +530,12 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): src_grid = pd.concat(src_grid) src_grid = self.comm.bcast(src_grid) + my_crs = pyproj.CRS.from_proj4("+proj=latlon") # Normalizing projections - dst_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) + dst_grid.to_crs(crs=my_crs, inplace=True) dst_grid.index.name = 'FID_dst' - src_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) + # src_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) + src_grid.to_crs(crs=my_crs, inplace=True) src_grid.index.name = 'FID_src' # Get intersected areas @@ -534,14 +544,28 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): # Calculate intersected areas and fractions intersection = pd.concat([src_grid.geometry[inp].reset_index(), dst_grid.geometry[res].reset_index()], axis=1, ignore_index=True) - intersection.columns = (list(src_grid.geometry[inp].reset_index().rename(columns={'geometry': 'geometry_src'}).columns) + list(dst_grid.geometry[res].reset_index().rename(columns={'geometry': 'geometry_dst'}).columns)) + intersection.columns = (list(src_grid.geometry[inp].reset_index().rename( + columns={'geometry': 'geometry_src'}).columns) + list(dst_grid.geometry[res].reset_index().rename( + columns={'geometry': 'geometry_dst'}).columns)) + + warnings.filterwarnings('ignore') + # NO WARNING ZONE intersection["src_area"] = src_grid.loc[intersection['FID_src'], 'geometry'].area.values + # intersection["src_area"] = gpd_geographic_area(src_grid.loc[intersection['FID_src'], 'geometry'].reset_index()) + + intersection['intersect_area'] = intersection.apply( + lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area, axis=1) - intersection['weight'] = intersection.apply( - lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area / x['src_area'], axis=1) + # intersection['geometry'] = intersection.apply(lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0), axis=1) + # intersection['intersect_area'] = gpd_geographic_area(intersection['geometry']) + + intersection.drop(intersection[intersection['intersect_area'] <= 0].index, inplace=True) + intersection['weight'] = intersection['intersect_area'] / intersection["src_area"] + + warnings.filterwarnings('default') # Format & Clean - intersection.drop(columns=["geometry_src", "geometry_dst", "src_area"], inplace=True) + intersection.drop(columns=["geometry_src", "geometry_dst", "src_area", "intersect_area"], inplace=True) intersection = intersection.set_index(['FID_dst', intersection.groupby('FID_dst').cumcount()]).rename_axis( ('FID', 'level')).sort_index() intersection.rename(columns={"FID_src": "idx"}, inplace=True) @@ -591,10 +615,6 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): weight_matrix.variables['weight']['data'] = weight_matrix.variables['weight']['data'].reshape(shape) weight_matrix.variables['idx']['data'] = weight_matrix.variables['idx']['data'].reshape(shape) - # TODO: Transpose do to FID od shapefile are transposed. - for aux_lev in weight_matrix.lev['data']: - weight_matrix.variables['weight']['data'][0, aux_lev] = weight_matrix.variables['weight']['data'][0, aux_lev].T - weight_matrix.variables['idx']['data'][0, aux_lev] = weight_matrix.variables['idx']['data'][0, aux_lev].T return weight_matrix diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index da964bf..65d5c9f 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -2577,7 +2577,7 @@ class Nes(object): intersection.columns = (list(self.shapefile.geometry[inp].reset_index().columns) + list(shapefile.geometry[res].reset_index().rename(columns={'geometry': 'geometry_shapefile', 'index': 'index_shapefile'}).columns)) - intersection['area'] = intersection.apply(lambda x: x['geometry'].intersection(x['geometry_shapefile']).buffer(0).area, + intersection['area'] = intersection.apply(lambda x: x['geometry'].intersection(x['geometry_shapefile']).buffer(0).area, axis=1) intersection['fraction'] = intersection.apply(lambda x: x['area'] / x['geometry'].area, axis=1) @@ -2806,7 +2806,7 @@ class Nes(object): self, new_levels, new_src_vertical=new_src_vertical, kind=kind, extrapolate=extrapolate, info=info) def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, - info=False, to_providentia=False): + info=False, to_providentia=False, only_create_wm=False): """ Horizontal methods from the current grid to another one. @@ -2817,15 +2817,17 @@ class Nes(object): weight_matrix_path : str, None Path to the weight matrix to read/create. kind : str - Kind of horizontal methods. choices = ['NearestNeighbour']. + Kind of horizontal methods. choices = ['NearestNeighbour', 'Conservative']. n_neighbours: int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. info: bool Indicates if you want to print extra info during the methods process. to_providentia : bool Indicates if we want the interpolated grid in Providentia format. + only_create_wm : bool + Indicates if you want to only create the Weight Matrix. """ return horizontal_interpolation.interpolate_horizontal( self, dst_grid, weight_matrix_path=weight_matrix_path, kind=kind, n_neighbours=n_neighbours, info=info, - to_providentia=to_providentia) + to_providentia=to_providentia, only_create_wm=only_create_wm) diff --git a/tests/test_interp_conservative.py b/tests/5.3-test_horiz_interp_conservative.py similarity index 90% rename from tests/test_interp_conservative.py rename to tests/5.3-test_horiz_interp_conservative.py index 3514000..7ca8f27 100644 --- a/tests/test_interp_conservative.py +++ b/tests/5.3-test_horiz_interp_conservative.py @@ -12,9 +12,9 @@ size = comm.Get_size() parallel_method = 'Y' -result_path = "Times_test_interp_conservative_{0}_{1:03d}.csv".format(parallel_method, size) +result_path = "Times_test_5.3_horiz_interp_conservative.py_{0}_{1:03d}.csv".format(parallel_method, size) result = pd.DataFrame(index=['read', 'calcul', 'write'], - columns=['Only interp', 'Create_WM', "Read_WM"]) + columns=['5.3.1.Only interp', '5.3.2.Create_WM', "5.3.3.Read_WM"]) # NAMEE src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc" @@ -23,7 +23,7 @@ var_list = ['O3'] # ====================================================================================================================== # ====================================== Only interp ===================================================== # ====================================================================================================================== -test_name = "Only interp" +test_name = '5.3.1.Only interp' print(test_name) # READING @@ -61,10 +61,15 @@ interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") comm.Barrier() result.loc['write', test_name] = timeit.default_timer() - st_time +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + # ====================================================================================================================== # ====================================== Create_WM ===================================================== # ====================================================================================================================== -test_name = "Create_WM" +test_name = '5.3.2.Create_WM' print(test_name) # READING @@ -106,10 +111,15 @@ interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") comm.Barrier() result.loc['write', test_name] = timeit.default_timer() - st_time +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + # ====================================================================================================================== # ====================================== Read_WM ===================================================== # ====================================================================================================================== -test_name = "Read_WM" +test_name = "5.3.3.Read_WM" print(test_name) # READING diff --git a/tests/test_bash_mn4.cmd b/tests/test_bash_mn4.cmd index 28ff46a..5e0c5c4 100644 --- a/tests/test_bash_mn4.cmd +++ b/tests/test_bash_mn4.cmd @@ -23,4 +23,5 @@ cd /gpfs/projects/bsc32/models/NES_master/tests #mpirun --mca mpi_warn_on_fork 0 -np 4 python 1-test_read_write_size.py #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2-test_read_write_projection.py #mpirun --mca mpi_warn_on_fork 0 -np 4 python 3-test_spatial_join.py -mpirun --mca mpi_warn_on_fork 0 -np 4 python 4-test_bounds.py +#mpirun --mca mpi_warn_on_fork 0 -np 4 python 4-test_bounds.py +mpirun --mca mpi_warn_on_fork 0 -np 4 python 5.3-test_horiz_interp_conservative.py diff --git a/tests/test_bash_nord3v2.cmd b/tests/test_bash_nord3v2.cmd index c255b60..c0cfb6e 100644 --- a/tests/test_bash_nord3v2.cmd +++ b/tests/test_bash_nord3v2.cmd @@ -23,3 +23,4 @@ cd /gpfs/projects/bsc32/models/NES_master/tests #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2-test_read_write_projection.py #mpirun --mca mpi_warn_on_fork 0 -np 4 python 3-test_spatial_join.py mpirun --mca mpi_warn_on_fork 0 -np 4 python 4-test_bounds.py +mpirun --mca mpi_warn_on_fork 0 -np 4 python 5.3-test_horiz_interp_conservative.py -- GitLab From 51ed5a21bb06fbb497f1880d98d137b6d6b27202 Mon Sep 17 00:00:00 2001 From: ctena Date: Thu, 9 Feb 2023 17:07:20 +0100 Subject: [PATCH 4/5] Horizontal Interpolation Conservative: Testing with parallel method 'T' & 'Y' --- nes/methods/horizontal_interpolation.py | 317 ++++++++++-------- nes/nc_projections/default_nes.py | 10 +- requirements.txt | 3 +- ...e.py => 3.2-test_horiz_interp_bilinear.py} | 102 ++++-- tests/3.3-test_horiz_interp_conservative.py | 227 +++++++++++++ 5 files changed, 499 insertions(+), 160 deletions(-) rename tests/{5.3-test_horiz_interp_conservative.py => 3.2-test_horiz_interp_bilinear.py} (60%) create mode 100644 tests/3.3-test_horiz_interp_conservative.py diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py index 8425bbb..5fc29d8 100644 --- a/nes/methods/horizontal_interpolation.py +++ b/nes/methods/horizontal_interpolation.py @@ -20,7 +20,7 @@ CONSERVATIVE_OPTS = ['Conservative', 'Area_Conservative', 'cons', 'conservative' def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, - info=False, to_providentia=False, only_create_wm=False): + info=False, to_providentia=False, only_create_wm=False, wm=None): """ Horizontal methods from one grid to another one. @@ -42,21 +42,26 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares Indicates if we want the interpolated grid in Providentia format. only_create_wm : bool Indicates if you want to only create the Weight Matrix. + wm : Nes + Weight matrix Nes File """ if info and self.master: print("Creating Weight Matrix") # Obtain weight matrix if self.parallel_method == 'T': - weights, idx = get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours) + weights, idx = get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, + only_create_wm, wm) elif self.parallel_method in ['Y', 'X']: - weights, idx = get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours) + weights, idx = get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, + only_create_wm, wm) else: raise NotImplemented("Parallel method {0} is not implemented yet for horizontal interpolations. Use 'T'".format( self.parallel_method)) if info and self.master: print("Weight Matrix done!") if only_create_wm: - return None + # weights for only_create is the WM NES object + return weights # idx[idx < 0] = np.nan idx = np.ma.masked_array(idx, mask=idx == -999) @@ -106,8 +111,8 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares for time in range(dst_shape[0]): for lev in range(dst_shape[1]): src_aux = get_src_data(self.comm, var_info['data'][time, lev], idx, self.parallel_method) - # src_aux = np.take(src_data[time, lev], idx) final_dst.variables[var_name]['data'][time, lev] = np.sum(weights * src_aux, axis=1) + if isinstance(dst_grid, nes.PointsNes): # Removing level axis if src_shape[1] != 1: @@ -189,7 +194,7 @@ def get_src_data(comm, var_data, idx, parallel_method): # noinspection DuplicatedCode -def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours): +def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, only_create, wm): """ To obtain the weights and source data index through the T axis. @@ -205,54 +210,63 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative']. n_neighbours : int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. + only_create : bool + Indicates if you want to only create the Weight Matrix. + wm : Nes + Weight matrix Nes File Returns ------- tuple Weights and source data index. """ + if wm is not None: + weight_matrix = wm - if weight_matrix_path is not None: - with FileLock(weight_matrix_path.replace('.nc', '.lock')): - if self.master: - if os.path.isfile(weight_matrix_path): + elif weight_matrix_path is not None: + with FileLock(weight_matrix_path + "{0:03d}.lock".format(self.rank)): + if os.path.isfile(weight_matrix_path): + if self.master: weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF) - if kind in NEAREST_OPTS: + else: + weight_matrix = True + if kind in NEAREST_OPTS: + if self.master: if len(weight_matrix.lev['data']) != n_neighbours: warn("The selected weight matrix does not have the same number of nearest neighbours." + "Re-calculating again but not saving it.") weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + else: + weight_matrix = True - else: + else: + if self.master: if kind in NEAREST_OPTS: - weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours, + wm_path=weight_matrix_path) elif kind in CONSERVATIVE_OPTS: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid, + wm_path=weight_matrix_path) else: raise NotImplementedError(kind) - if weight_matrix_path is not None: - weight_matrix.to_netcdf(weight_matrix_path) - else: - weight_matrix = None - if self.master: - if os.path.exists(weight_matrix_path.replace('.nc', '.lock')): - os.remove(weight_matrix_path.replace('.nc', '.lock')) + else: + weight_matrix = True + + if os.path.exists(weight_matrix_path + "{0:03d}.lock".format(self.rank)): + os.remove(weight_matrix_path + "{0:03d}.lock".format(self.rank)) else: - if kind in NEAREST_OPTS: - if self.master: + if self.master: + if kind in NEAREST_OPTS: weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: - weight_matrix = None - elif kind in CONSERVATIVE_OPTS: - if self.master: - if self.master: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) - else: - weight_matrix = None - else: - weight_matrix = None + raise NotImplementedError(kind) else: - raise NotImplementedError(kind) + weight_matrix = True + + if only_create: + return weight_matrix, None if self.master: if kind in NEAREST_OPTS: @@ -261,7 +275,7 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour np.array(weight_matrix.variables['weight']['data'], dtype=np.float64).sum(axis=1), dtype=np.float64) else: - weights = np.array(weight_matrix.variables['weight']['data'][0], dtype=np.float64) + weights = np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int) else: weights = None @@ -274,7 +288,7 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour # noinspection DuplicatedCode -def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours): +def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, only_create, wm): """ To obtain the weights and source data index through the X or Y axis. @@ -290,52 +304,63 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative']. n_neighbours : int Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4. + only_create : bool + Indicates if you want to only create the Weight Matrix. + wm : Nes + Weight matrix Nes File Returns ------- tuple Weights and source data index. """ - if isinstance(dst_grid, nes.PointsNes) and weight_matrix_path is not None: if self.master: warn("To point weight matrix cannot be saved.") weight_matrix_path = None - if weight_matrix_path is not None: - with FileLock(weight_matrix_path.replace('.nc', '.lock')): - if self.master: - if os.path.isfile(weight_matrix_path): + if wm is not None: + weight_matrix = wm + + elif weight_matrix_path is not None: + with FileLock(weight_matrix_path + "{0:03d}.lock".format(self.rank)): + if os.path.isfile(weight_matrix_path): + if self.master: weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF) - if kind in NEAREST_OPTS: - if len(weight_matrix.lev['data']) != n_neighbours: - warn("The selected weight matrix does not have the same number of nearest neighbours." + - "Re-calculating again but not saving it.") - weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) else: - if kind in NEAREST_OPTS: - weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) - elif kind in CONSERVATIVE_OPTS: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) - else: - raise NotImplementedError(kind) - if weight_matrix_path is not None: - weight_matrix.to_netcdf(weight_matrix_path) + weight_matrix = True + if kind in NEAREST_OPTS: + if len(weight_matrix.lev['data']) != n_neighbours: + warn("The selected weight matrix does not have the same number of nearest neighbours." + + "Re-calculating again but not saving it.") + if self.master: + weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + else: + weight_matrix = True else: - weight_matrix = None - if self.master: - if os.path.exists(weight_matrix_path.replace('.nc', '.lock')): - os.remove(weight_matrix_path.replace('.nc', '.lock')) + if kind in NEAREST_OPTS: + if self.master: + weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours, + wm_path=weight_matrix_path) + else: + weight_matrix = True + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid, wm_path=weight_matrix_path) + else: + raise NotImplementedError(kind) + + if os.path.exists(weight_matrix_path + "{0:03d}.lock".format(self.rank)): + os.remove(weight_matrix_path + "{0:03d}.lock".format(self.rank)) else: - if self.master: - if kind in NEAREST_OPTS: - weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) - elif kind in CONSERVATIVE_OPTS: - weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) - else: - raise NotImplementedError(kind) + if kind in NEAREST_OPTS: + weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours) + elif kind in CONSERVATIVE_OPTS: + weight_matrix = create_area_conservative_weight_matrix(self, dst_grid) else: - weight_matrix = None + raise NotImplementedError(kind) + + if only_create: + return weight_matrix, None # Normalize to 1 if self.master: @@ -345,13 +370,14 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou dtype=np.float64) else: weights = np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) - idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int) + idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=np.int64) else: weights = None idx = None weights = self.comm.bcast(weights, root=0) idx = self.comm.bcast(idx, root=0) + # if isinstance(dst_grid, nes.PointsNes): # print("weights 1 ->", weights.shape) # print("idx 1 ->", idx.shape) @@ -386,7 +412,6 @@ def read_weight_matrix(weight_matrix_path, comm=None, parallel_method='T'): nes.Nes Weight matrix. """ - weight_matrix = nes.open_netcdf(path=weight_matrix_path, comm=comm, parallel_method=parallel_method, balanced=True) weight_matrix.load() @@ -397,7 +422,7 @@ def read_weight_matrix(weight_matrix_path, comm=None, parallel_method='T'): return weight_matrix -def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): +def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, wm_path=None, info=False): """ To create the weight matrix with the nearest neighbours method. @@ -485,11 +510,13 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False): weight_matrix.variables['idx'] = {'data': idx_transf, 'units': ''} weight_matrix.lev = {'data': np.arange(inverse_dists_transf.shape[1]), 'units': ''} weight_matrix._lev = {'data': np.arange(inverse_dists_transf.shape[1]), 'units': ''} + if wm_path is not None: + weight_matrix.to_netcdf(wm_path) return weight_matrix -def create_area_conservative_weight_matrix(self, dst_nes, info=False): +def create_area_conservative_weight_matrix(self, dst_nes, wm_path=None, info=False): """ To create the weight matrix with the area conservative method. @@ -508,7 +535,6 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): nes.Nes Weight matrix. """ - if info and self.master: print("\tCreating area conservative Weight Matrix") sys.stdout.flush() @@ -521,6 +547,7 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): if self.shapefile is None: self.create_shapefile() src_grid = copy.deepcopy(self.shapefile) + if self.parallel_method == 'T': # All process has the same shapefile pass @@ -533,88 +560,110 @@ def create_area_conservative_weight_matrix(self, dst_nes, info=False): my_crs = pyproj.CRS.from_proj4("+proj=latlon") # Normalizing projections dst_grid.to_crs(crs=my_crs, inplace=True) - dst_grid.index.name = 'FID_dst' + dst_grid['FID_dst'] = dst_grid.index + dst_grid = dst_grid.reset_index() + # src_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True) src_grid.to_crs(crs=my_crs, inplace=True) - src_grid.index.name = 'FID_src' + src_grid['FID_src'] = src_grid.index + + src_grid = src_grid.reset_index() + + if info and self.master: + print("\t\tGrids created and ready to interpolate") + sys.stdout.flush() # Get intersected areas inp, res = dst_grid.sindex.query_bulk(src_grid.geometry, predicate='intersects') # Calculate intersected areas and fractions - intersection = pd.concat([src_grid.geometry[inp].reset_index(), dst_grid.geometry[res].reset_index()], axis=1, - ignore_index=True) - intersection.columns = (list(src_grid.geometry[inp].reset_index().rename( - columns={'geometry': 'geometry_src'}).columns) + list(dst_grid.geometry[res].reset_index().rename( - columns={'geometry': 'geometry_dst'}).columns)) + intersection = pd.DataFrame(columns=["FID_src", "FID_dst"]) + intersection['INP'] = np.array(inp, dtype=np.uint32) + intersection['RES'] = np.array(res, dtype=np.uint32) + intersection['FID_src'] = np.array(src_grid.loc[inp, 'FID_src'], dtype=np.uint32) + intersection['FID_dst'] = np.array(dst_grid.loc[res, 'FID_dst'], dtype=np.uint32) - warnings.filterwarnings('ignore') - # NO WARNING ZONE - intersection["src_area"] = src_grid.loc[intersection['FID_src'], 'geometry'].area.values - # intersection["src_area"] = gpd_geographic_area(src_grid.loc[intersection['FID_src'], 'geometry'].reset_index()) + intersection['geometry_src'] = src_grid.loc[inp, 'geometry'].values + intersection['geometry_dst'] = dst_grid.loc[res, 'geometry'].values - intersection['intersect_area'] = intersection.apply( - lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area, axis=1) + if True: + # No Warnings Zone + warnings.filterwarnings('ignore') + intersection['intersect_area'] = intersection.apply( + lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area, axis=1) + intersection.drop(intersection[intersection['intersect_area'] <= 0].index, inplace=True) - # intersection['geometry'] = intersection.apply(lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0), axis=1) - # intersection['intersect_area'] = gpd_geographic_area(intersection['geometry']) + intersection["src_area"] = src_grid.loc[intersection['FID_src'], 'geometry'].area.values - intersection.drop(intersection[intersection['intersect_area'] <= 0].index, inplace=True) - intersection['weight'] = intersection['intersect_area'] / intersection["src_area"] + warnings.filterwarnings('default') - warnings.filterwarnings('default') + intersection['weight'] = intersection['intersect_area'] / intersection["src_area"] # Format & Clean intersection.drop(columns=["geometry_src", "geometry_dst", "src_area", "intersect_area"], inplace=True) - intersection = intersection.set_index(['FID_dst', intersection.groupby('FID_dst').cumcount()]).rename_axis( - ('FID', 'level')).sort_index() - intersection.rename(columns={"FID_src": "idx"}, inplace=True) - # Initialising weight matrix - weight_matrix = dst_nes.copy() - weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] - weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] - weight_matrix._time_bnds = None - weight_matrix.time_bnds = None - weight_matrix.last_level = None - weight_matrix.first_level = 0 - weight_matrix.hours_start = 0 - weight_matrix.hours_end = 0 - - weight_matrix.set_communicator(MPI.COMM_SELF) + if info and self.master: + print("\t\tWeights calculated. Formatting weight matrix.") + sys.stdout.flush() - weight_matrix.set_levels({'data': np.arange(intersection.index.get_level_values('level').max() + 1), - 'dimensions': ('lev',), - 'units': '', - 'positive': 'up'}) - - # Creating Weight matrix empty variables - if len(weight_matrix.lat['data'].shape) == 1: - shape = (1, len(weight_matrix.lev['data']), - weight_matrix.lat['data'].shape[0], weight_matrix.lon['data'].shape[0],) - shape_flat = (1, len(weight_matrix.lev['data']), - weight_matrix.lat['data'].shape[0] * weight_matrix.lon['data'].shape[0],) + # Initialising weight matrix + if self.parallel_method != 'T': + intersection = self.comm.gather(intersection, root=0) + if self.master: + if self.parallel_method != 'T': + intersection = pd.concat(intersection) + intersection = intersection.set_index(['FID_dst', intersection.groupby('FID_dst').cumcount()]).rename_axis( + ('FID', 'level')).sort_index() + intersection.rename(columns={"FID_src": "idx"}, inplace=True) + weight_matrix = dst_nes.copy() + weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] + weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)] + weight_matrix._time_bnds = None + weight_matrix.time_bnds = None + weight_matrix.last_level = None + weight_matrix.first_level = 0 + weight_matrix.hours_start = 0 + weight_matrix.hours_end = 0 + + weight_matrix.set_communicator(MPI.COMM_SELF) + + weight_matrix.set_levels({'data': np.arange(intersection.index.get_level_values('level').max() + 1), + 'dimensions': ('lev',), + 'units': '', + 'positive': 'up'}) + + # Creating Weight matrix empty variables + if len(weight_matrix._lat['data'].shape) == 1: + shape = (1, len(weight_matrix.lev['data']), + weight_matrix._lat['data'].shape[0], weight_matrix._lon['data'].shape[0],) + shape_flat = (1, len(weight_matrix.lev['data']), + weight_matrix._lat['data'].shape[0] * weight_matrix._lon['data'].shape[0],) + else: + shape = (1, len(weight_matrix.lev['data']), + weight_matrix._lat['data'].shape[0], weight_matrix._lat['data'].shape[1],) + shape_flat = (1, len(weight_matrix.lev['data']), + weight_matrix._lat['data'].shape[0] * weight_matrix._lat['data'].shape[1],) + + weight_matrix.variables['weight'] = {'data': np.empty(shape_flat), 'units': '-'} + weight_matrix.variables['weight']['data'][:] = -999 + weight_matrix.variables['idx'] = {'data': np.empty(shape_flat), 'units': '-'} + weight_matrix.variables['idx']['data'][:] = -999 + + # Filling Weight matrix variables + for aux_lev in weight_matrix.lev['data']: + aux_data = intersection.xs(level='level', key=aux_lev) + weight_matrix.variables['weight']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'weight'].values + weight_matrix.variables['idx']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'idx'].values + # Re-shaping + weight_matrix.variables['weight']['data'] = weight_matrix.variables['weight']['data'].reshape(shape) + weight_matrix.variables['idx']['data'] = weight_matrix.variables['idx']['data'].reshape(shape) + if wm_path is not None: + if info and self.master: + print("\t\tWeight matrix saved at {0}".format(wm_path)) + sys.stdout.flush() + weight_matrix.to_netcdf(wm_path) else: - shape = (1, len(weight_matrix.lev['data']), - weight_matrix.lat['data'].shape[0], weight_matrix.lat['data'].shape[1],) - shape_flat = (1, len(weight_matrix.lev['data']), - weight_matrix.lat['data'].shape[0] * weight_matrix.lat['data'].shape[1],) - - weight_matrix.variables['weight'] = {'data': np.empty(shape_flat), 'units': '-'} - weight_matrix.variables['weight']['data'][:] = -999 - weight_matrix.variables['idx'] = {'data': np.empty(shape_flat), 'units': '-'} - weight_matrix.variables['idx']['data'][:] = -999 - - # Filling Weight matrix variables - for aux_lev in weight_matrix.lev['data']: - aux_data = intersection.xs(level='level', key=aux_lev) - weight_matrix.variables['weight']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'weight'] - weight_matrix.variables['idx']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'idx'] - - # Re-shaping - weight_matrix.variables['weight']['data'] = weight_matrix.variables['weight']['data'].reshape(shape) - weight_matrix.variables['idx']['data'] = weight_matrix.variables['idx']['data'].reshape(shape) - + weight_matrix = True return weight_matrix diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 3e7b3be..0e91a94 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -1010,7 +1010,7 @@ class Nes(object): 't_min': None, 't_max': None} idx = self.get_idx_intervals() - + if self.parallel_method == 'Y': y_len = idx['idx_y_max'] - idx['idx_y_min'] if y_len < self.size: @@ -1048,7 +1048,7 @@ class Nes(object): axis_limits['t_max'] = idx['idx_t_max'] elif self.parallel_method == 'T': - t_len = idx['idx_t_min'] - idx['idx_t_max'] + t_len = idx['idx_t_max'] - idx['idx_t_min'] if t_len < self.size: raise IndexError('More processors (size={0}) selected than T elements (size={1})'.format( self.size, t_len)) @@ -2935,7 +2935,7 @@ class Nes(object): self, new_levels, new_src_vertical=new_src_vertical, kind=kind, extrapolate=extrapolate, info=info) def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, - info=False, to_providentia=False, only_create_wm=False): + info=False, to_providentia=False, only_create_wm=False, wm=None): """ Horizontal methods from the current grid to another one. @@ -2955,11 +2955,13 @@ class Nes(object): Indicates if we want the interpolated grid in Providentia format. only_create_wm : bool Indicates if you want to only create the Weight Matrix. + wm : Nes + Weight matrix Nes File """ return horizontal_interpolation.interpolate_horizontal( self, dst_grid, weight_matrix_path=weight_matrix_path, kind=kind, n_neighbours=n_neighbours, info=info, - to_providentia=to_providentia, only_create_wm=only_create_wm) + to_providentia=to_providentia, only_create_wm=only_create_wm, wm=wm) def calculate_grid_area(self): """ diff --git a/requirements.txt b/requirements.txt index 3e6c91a..a40856e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1,9 @@ pycodestyle>=2.10.0 geopandas>=0.10.2 +rtree>=0.9.0 pandas>=1.3.5 netcdf4>=1.6.2 -numpy==1.21.6 +numpy>=1.20.0 pyproj~=3.2.1 setuptools>=66.1.1 pytest>=7.2.1 diff --git a/tests/5.3-test_horiz_interp_conservative.py b/tests/3.2-test_horiz_interp_bilinear.py similarity index 60% rename from tests/5.3-test_horiz_interp_conservative.py rename to tests/3.2-test_horiz_interp_bilinear.py index 7ca8f27..77bddad 100644 --- a/tests/5.3-test_horiz_interp_conservative.py +++ b/tests/3.2-test_horiz_interp_bilinear.py @@ -10,11 +10,11 @@ comm = MPI.COMM_WORLD rank = comm.Get_rank() size = comm.Get_size() -parallel_method = 'Y' +parallel_method = 'T' -result_path = "Times_test_5.3_horiz_interp_conservative.py_{0}_{1:03d}.csv".format(parallel_method, size) +result_path = "Times_test_3.2_horiz_interp_bilinear_{0}_{1:03d}.csv".format(parallel_method, size) result = pd.DataFrame(index=['read', 'calcul', 'write'], - columns=['5.3.1.Only interp', '5.3.2.Create_WM', "5.3.3.Read_WM"]) + columns=['3.2.1.Only interp', '3.2.2.Create_WM', "3.2.3.Use_WM", "3.2.4.Read_WM"]) # NAMEE src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc" @@ -23,8 +23,9 @@ var_list = ['O3'] # ====================================================================================================================== # ====================================== Only interp ===================================================== # ====================================================================================================================== -test_name = '5.3.1.Only interp' -print(test_name) +test_name = '3.2.1.NN_Only interp' +if rank == 0: + print(test_name) # READING st_time = timeit.default_timer() @@ -46,18 +47,19 @@ inc_y = 4000 x_0 = -807847.688 y_0 = -797137.125 dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, - nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) comm.Barrier() result.loc['read', test_name] = timeit.default_timer() - st_time st_time = timeit.default_timer() -interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative') +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN') comm.Barrier() result.loc['calcul', test_name] = timeit.default_timer() - st_time st_time = timeit.default_timer() -interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size)) comm.Barrier() result.loc['write', test_name] = timeit.default_timer() - st_time @@ -69,8 +71,9 @@ sys.stdout.flush() # ====================================================================================================================== # ====================================== Create_WM ===================================================== # ====================================================================================================================== -test_name = '5.3.2.Create_WM' -print(test_name) +test_name = '3.2.2.NN_Create_WM' +if rank == 0: + print(test_name) # READING st_time = timeit.default_timer() @@ -92,22 +95,70 @@ inc_y = 4000 x_0 = -807847.688 y_0 = -797137.125 dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, - nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) comm.Barrier() result.loc['read', test_name] = timeit.default_timer() - st_time # Cleaning WM -if os.path.exists("NAMEE_to_IP.nc"): - os.remove("NAMEE_to_IP.nc") +if os.path.exists("NN_WM_NAMEE_to_IP.nc") and rank == 0: + os.remove("NN_WM_NAMEE_to_IP.nc") +comm.Barrier() st_time = timeit.default_timer() -interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="NAMEE_to_IP.nc") +wm_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN', info=True, + weight_matrix_path="NN_WM_NAMEE_to_IP.nc", only_create_wm=True) comm.Barrier() result.loc['calcul', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ====================================== Use_WM ===================================================== +# ====================================================================================================================== +test_name = "3.2.3.NN_Use_WM" +if rank == 0: + print(test_name) + +# READING st_time = timeit.default_timer() -interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 + +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN', wm=wm_nes) +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) comm.Barrier() result.loc['write', test_name] = timeit.default_timer() - st_time @@ -119,8 +170,9 @@ sys.stdout.flush() # ====================================================================================================================== # ====================================== Read_WM ===================================================== # ====================================================================================================================== -test_name = "5.3.3.Read_WM" -print(test_name) +test_name = "3.2.4.NN_Read_WM" +if rank == 0: + print(test_name) # READING st_time = timeit.default_timer() @@ -143,19 +195,27 @@ x_0 = -807847.688 y_0 = -797137.125 dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, - nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method) + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) comm.Barrier() result.loc['read', test_name] = timeit.default_timer() - st_time st_time = timeit.default_timer() -interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="NAMEE_to_IP.nc") +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN', + weight_matrix_path="NN_WM_NAMEE_to_IP.nc") comm.Barrier() result.loc['calcul', test_name] = timeit.default_timer() - st_time st_time = timeit.default_timer() -interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +interp_nes.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) comm.Barrier() result.loc['write', test_name] = timeit.default_timer() - st_time -result.to_csv(result_path) \ No newline at end of file +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +if rank == 0: + result.to_csv(result_path) diff --git a/tests/3.3-test_horiz_interp_conservative.py b/tests/3.3-test_horiz_interp_conservative.py new file mode 100644 index 0000000..a6af099 --- /dev/null +++ b/tests/3.3-test_horiz_interp_conservative.py @@ -0,0 +1,227 @@ +#!/usr/bin/env python +import sys +import timeit +import pandas as pd +from mpi4py import MPI +from nes import * +import os, sys + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +parallel_method = 'T' + +result_path = "Times_test_3.3_horiz_interp_conservative.py_{0}_{1:03d}.csv".format(parallel_method, size) +result = pd.DataFrame(index=['read', 'calcul', 'write'], + columns=['3.3.1.Only interp', '3.3.2.Create_WM', "3.3.3.Use_WM", "3.3.4.Read_WM"]) + +# NAMEE +src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc" +var_list = ['O3'] + +# ====================================================================================================================== +# ====================================== Only interp ===================================================== +# ====================================================================================================================== +test_name = '3.3.1.Only interp' +if rank == 0: + print(test_name) + +# READING final_dst.variables[var_name]['data'][time, lev] = np.sum(weights * src_aux, axis=1) + +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative') +# interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path='T_WM.nc') + +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size), serial=True) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ====================================== Create_WM ===================================================== +# ====================================================================================================================== +test_name = '3.3.2.Create_WM' +if rank == 0: + print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +# Cleaning WM +if os.path.exists("CONS_WM_NAMEE_to_IP.nc") and rank == 0: + os.remove("CONS_WM_NAMEE_to_IP.nc") +comm.Barrier() + +st_time = timeit.default_timer() + +wm_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', info=True, + weight_matrix_path="CONS_WM_NAMEE_to_IP.nc", only_create_wm=True) +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +# st_time = timeit.default_timer() +# interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc") +# comm.Barrier() +# result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ====================================== Use_WM ===================================================== +# ====================================================================================================================== +test_name = "3.3.3.Use_WM" +if rank == 0: + print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 + +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', wm=wm_nes) +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ====================================== Read_WM ===================================================== +# ====================================================================================================================== +test_name = "3.3.4.Read_WM" +if rank == 0: + print(test_name) + +# READING +st_time = timeit.default_timer() + +# Source data +src_nes = open_netcdf(src_path, parallel_method=parallel_method) +src_nes.keep_vars(var_list) +src_nes.load() + +# Destination Grid +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 + +dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method, + times=src_nes.get_full_times()) +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() + +interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="CONS_WM_NAMEE_to_IP.nc") +comm.Barrier() +result.loc['calcul', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +if rank == 0: + result.to_csv(result_path) -- GitLab From 7c6e6900081fb696d4896ebed4671455d7c6a25b Mon Sep 17 00:00:00 2001 From: Alba Vilanova Date: Tue, 14 Feb 2023 06:31:54 +0100 Subject: [PATCH 5/5] Add conservative interpolation test and update interpolation tests order --- .../4.1.Vertical_Interpolation.ipynb | 2 +- .../4.2.Horizontal_Interpolation.ipynb | 4 +- .../4.3.Conservative_Interpolation.ipynb | 396 ++++++++++++++++++ ...nb => 4.4.Providentia_Interpolation.ipynb} | 0 ....5.NES_vs_Providentia_Interpolation.ipynb} | 0 5 files changed, 399 insertions(+), 3 deletions(-) create mode 100644 tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb rename tutorials/4.Interpolation/{4.3.Providentia_Interpolation.ipynb => 4.4.Providentia_Interpolation.ipynb} (100%) rename tutorials/4.Interpolation/{4.4.NES_vs_Providentia_Interpolation.ipynb => 4.5.NES_vs_Providentia_Interpolation.ipynb} (100%) diff --git a/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb b/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb index 38b1ac5..d126994 100644 --- a/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb +++ b/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# How to interpolate vertically" + "# Vertical interpolation" ] }, { diff --git a/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb b/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb index 5f0379f..63264ad 100644 --- a/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb +++ b/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# How to interpolate horizontally" + "# Horizontal interpolation" ] }, { @@ -593,7 +593,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8 (default, Aug 12 2021, 07:06:15) \n[GCC 8.4.1 20200928 (Red Hat 8.4.1-1)]" + "version": "3.7.4" }, "vscode": { "interpreter": { diff --git a/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb b/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb new file mode 100644 index 0000000..6177f4c --- /dev/null +++ b/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb @@ -0,0 +1,396 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Horizontal conservative interpolation" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from nes import *" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Read data to interpolate" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Original path: /esarchive/exp/monarch/a4dd/original_files/000/2022111512/MONARCH_d01_2022111512.nc\n", + "# Rotated grid from MONARCH (NAMEE)\n", + "source_path = \"/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc\"" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "source_grid = open_netcdf(path=source_path, info=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': masked_array(\n", + " data=[[16.350338, 16.43293 , 16.515146, ..., 16.515146, 16.43293 ,\n", + " 16.350338],\n", + " [16.527426, 16.610239, 16.692677, ..., 16.692677, 16.610243,\n", + " 16.527426],\n", + " [16.704472, 16.787508, 16.870167, ..., 16.870167, 16.78751 ,\n", + " 16.704472],\n", + " ...,\n", + " [58.32095 , 58.472683, 58.62431 , ..., 58.62431 , 58.472683,\n", + " 58.32095 ],\n", + " [58.426285, 58.5782 , 58.730026, ..., 58.730026, 58.5782 ,\n", + " 58.426285],\n", + " [58.530792, 58.6829 , 58.83492 , ..., 58.83492 , 58.682903,\n", + " 58.530792]],\n", + " mask=False,\n", + " fill_value=1e+20,\n", + " dtype=float32),\n", + " 'dimensions': ('rlat', 'rlon'),\n", + " 'long_name': 'latitude',\n", + " 'units': 'degrees_north',\n", + " 'standard_name': 'latitude',\n", + " 'coordinates': 'lon lat'}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "source_grid.lat" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': masked_array(\n", + " data=[[-22.181265, -22.016672, -21.851799, ..., 41.851795, 42.016666,\n", + " 42.18126 ],\n", + " [-22.27818 , -22.113186, -21.947905, ..., 41.9479 , 42.113174,\n", + " 42.27817 ],\n", + " [-22.375267, -22.209873, -22.04419 , ..., 42.044186, 42.209873,\n", + " 42.375263],\n", + " ...,\n", + " [-67.57767 , -67.397064, -67.21535 , ..., 87.21534 , 87.39706 ,\n", + " 87.57766 ],\n", + " [-67.90188 , -67.72247 , -67.54194 , ..., 87.54194 , 87.72246 ,\n", + " 87.90187 ],\n", + " [-68.228035, -68.04982 , -67.870514, ..., 87.87051 , 88.04982 ,\n", + " 88.228035]],\n", + " mask=False,\n", + " fill_value=1e+20,\n", + " dtype=float32),\n", + " 'dimensions': ('rlat', 'rlon'),\n", + " 'long_name': 'longitude',\n", + " 'units': 'degrees_east',\n", + " 'standard_name': 'longitude',\n", + " 'coordinates': 'lon lat'}" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "source_grid.lon" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Rank 000: Loading O3 var (1/1)\n", + "Rank 000: Loaded O3 var ((37, 24, 271, 351))\n" + ] + } + ], + "source": [ + "source_grid.keep_vars('O3')\n", + "source_grid.load()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Create destination grid" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "lat_1 = 37\n", + "lat_2 = 43\n", + "lon_0 = -3\n", + "lat_0 = 40\n", + "nx = 397\n", + "ny = 397\n", + "inc_x = 4000\n", + "inc_y = 4000\n", + "x_0 = -807847.688\n", + "y_0 = -797137.125\n", + "dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,\n", + " nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, times=source_grid.get_full_times())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Interpolation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.1. Without creating weight matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "interp_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative')" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': array([[32.49460816, 32.49803615, 32.50144726, ..., 32.52460207,\n", + " 32.52130788, 32.51799681],\n", + " [32.53024662, 32.5336765 , 32.5370895 , ..., 32.5602571 ,\n", + " 32.55696109, 32.55364819],\n", + " [32.5658877 , 32.56931947, 32.57273435, ..., 32.59591474,\n", + " 32.59261692, 32.58930218],\n", + " ...,\n", + " [46.60231473, 46.60653579, 46.6107361 , ..., 46.63924881,\n", + " 46.63519229, 46.63111499],\n", + " [46.63791957, 46.64214277, 46.6463452 , ..., 46.67487234,\n", + " 46.67081377, 46.66673441],\n", + " [46.67352141, 46.67774674, 46.68195131, ..., 46.71049289,\n", + " 46.70643226, 46.70235083]])}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interp_nes.lat" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': array([[-11.56484687, -11.52259289, -11.48033507, ..., 5.18764453,\n", + " 5.22992847, 5.27220869],\n", + " [-11.56892309, -11.52664925, -11.48437156, ..., 5.1915433 ,\n", + " 5.23384715, 5.27614727],\n", + " [-11.57300318, -11.53070946, -11.48841188, ..., 5.19544578,\n", + " 5.23776955, 5.2800896 ],\n", + " ...,\n", + " [-13.53865445, -13.48682654, -13.4349915 , ..., 7.07590089,\n", + " 7.12778439, 7.17966098],\n", + " [-13.54481681, -13.49295916, -13.44109437, ..., 7.08179741,\n", + " 7.13371074, 7.18561716],\n", + " [-13.55098635, -13.49909893, -13.44720434, ..., 7.0877008 ,\n", + " 7.139644 , 7.19158028]])}" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interp_nes.lon" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.2. Creating weight matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating Weight Matrix\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/esarchive/scratch/avilanova/software/NES/nes/nc_projections/default_nes.py:2035: DeprecationWarning: tostring() is deprecated. Use tobytes() instead.\n", + " time_var.units, time_var.calendar)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Weight Matrix done!\n" + ] + } + ], + "source": [ + "wm_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', info=True,\n", + " weight_matrix_path=\"CONS_WM_NAMEE_to_IP.nc\", only_create_wm=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "interp_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', wm=wm_nes)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': array([[32.49460816, 32.49803615, 32.50144726, ..., 32.52460207,\n", + " 32.52130788, 32.51799681],\n", + " [32.53024662, 32.5336765 , 32.5370895 , ..., 32.5602571 ,\n", + " 32.55696109, 32.55364819],\n", + " [32.5658877 , 32.56931947, 32.57273435, ..., 32.59591474,\n", + " 32.59261692, 32.58930218],\n", + " ...,\n", + " [46.60231473, 46.60653579, 46.6107361 , ..., 46.63924881,\n", + " 46.63519229, 46.63111499],\n", + " [46.63791957, 46.64214277, 46.6463452 , ..., 46.67487234,\n", + " 46.67081377, 46.66673441],\n", + " [46.67352141, 46.67774674, 46.68195131, ..., 46.71049289,\n", + " 46.70643226, 46.70235083]])}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interp_nes.lat" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'data': array([[-11.56484687, -11.52259289, -11.48033507, ..., 5.18764453,\n", + " 5.22992847, 5.27220869],\n", + " [-11.56892309, -11.52664925, -11.48437156, ..., 5.1915433 ,\n", + " 5.23384715, 5.27614727],\n", + " [-11.57300318, -11.53070946, -11.48841188, ..., 5.19544578,\n", + " 5.23776955, 5.2800896 ],\n", + " ...,\n", + " [-13.53865445, -13.48682654, -13.4349915 , ..., 7.07590089,\n", + " 7.12778439, 7.17966098],\n", + " [-13.54481681, -13.49295916, -13.44109437, ..., 7.08179741,\n", + " 7.13371074, 7.18561716],\n", + " [-13.55098635, -13.49909893, -13.44720434, ..., 7.0877008 ,\n", + " 7.139644 , 7.19158028]])}" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interp_nes.lon" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb b/tutorials/4.Interpolation/4.4.Providentia_Interpolation.ipynb similarity index 100% rename from tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb rename to tutorials/4.Interpolation/4.4.Providentia_Interpolation.ipynb diff --git a/tutorials/4.Interpolation/4.4.NES_vs_Providentia_Interpolation.ipynb b/tutorials/4.Interpolation/4.5.NES_vs_Providentia_Interpolation.ipynb similarity index 100% rename from tutorials/4.Interpolation/4.4.NES_vs_Providentia_Interpolation.ipynb rename to tutorials/4.Interpolation/4.5.NES_vs_Providentia_Interpolation.ipynb -- GitLab