diff --git a/environment.yml b/environment.yml index 157b849..a2b85e9 100644 --- a/environment.yml +++ b/environment.yml @@ -14,6 +14,8 @@ dependencies: - cartopy - shapely - cdo + - xesmf + - sparse - pip: - smmregrid diff --git a/src/fint/fint.py b/src/fint/fint.py index 64c8a55..8f5a4d4 100644 --- a/src/fint/fint.py +++ b/src/fint/fint.py @@ -6,6 +6,8 @@ import numpy as np import pandas as pd import xarray as xr +import xesmf as xe +import sparse import subprocess from smmregrid import Regridder from scipy.interpolate import ( @@ -266,25 +268,26 @@ def interpolate_linear_scipy(data_in, x2, y2, lon2, lat2): return interpolated -def interpolate_cdo(target_grid,gridfile,original_file,output_file,variable_name,interpolation, mask_zero=True): +def interpolate_cdo(target_grid,gridfile,original_file,output_file,variable_name,interpolation,weights_file_path, mask_zero=True): """ - Interpolate a variable in a file using CDO (Climate Data Operators). + Interpolate a climate variable in a NetCDF file using Climate Data Operators (CDO). Args: - target_grid (str): Path to the target grid file. - gridfile (str): Path to the grid file associated with the original file. - original_file (str): Path to the original file containing the variable to be interpolated. - output_file (str): Path to the output file where the interpolated variable will be saved. + target_grid (str): Path to the target grid file (NetCDF format) for interpolation. + gridfile (str): Path to the grid file (NetCDF format) associated with the original data. + original_file (str): Path to the original NetCDF file containing the variable to be interpolated. + output_file (str): Path to the output NetCDF file where the interpolated variable will be saved. variable_name (str): Name of the variable to be interpolated. - interpolation (str): Interpolation method to be used (cdo_remapcon,cdo_remaplaf,cdo_remapnn, cdo_remapdis). + interpolation (str): Interpolation method to be used (e.g., 'remapcon', 'remaplaf', 'remapnn', 'remapdis'). + weights_file_path (str): Path to the weights file generated by CDO for interpolation. + mask_zero (bool, optional): Whether to mask zero values in the output. Default is True. Returns: np.ndarray: Interpolated variable data as a NumPy array. """ - command = [ "cdo", - f"-{interpolation.split('_')[1]},{target_grid}", + f"-remap,{target_grid},{weights_file_path}", f"-setgrid,{gridfile}", f"{original_file}", f"{output_file}" @@ -299,34 +302,88 @@ def interpolate_cdo(target_grid,gridfile,original_file,output_file,variable_name os.remove(output_file) return interpolated -def generate_cdo_weights(target_grid,gridfile,original_file,output_file,interpolation): +def generate_cdo_weights(target_grid,gridfile,original_file,output_file,interpolation, save = False): """ - Generate CDO weights for interpolation using ssmregrid. + Generate CDO interpolation weights for smmregrid and cdo interpolation. Args: - target_grid (str): Path to the target grid file. - gridfile (str): Path to the grid file associated with the original file. - original_file (str): Path to the original file containing the data to be remapped. - output_file (str): Path to the output file where the weights will be saved. + target_grid (str): Path to the target grid file (NetCDF format). + gridfile (str): Path to the grid file (NetCDF format) associated with the original data. + original_file (str): Path to the original NetCDF file containing the data to be remapped. + output_file (str): Path to the output NetCDF file where the weights will be saved. + interpolation (str): Interpolation method to be used. + save (bool, optional): Whether to save the weights file. Default is False. Returns: - xr.Dataset: Generated weights as an xarray Dataset. + xr.Dataset: Generated interpolation weights as an xarray Dataset. """ + + int_method = interpolation.split('_')[-1][5:] + command = [ "cdo", - f"-gen{interpolation.split('_')[-1]},{target_grid}", + f"-gen{int_method},{target_grid}", f"-setgrid,{gridfile}", f"{original_file}", f"{output_file}" ] - # Execute the command subprocess.run(command) - + weights = xr.open_dataset(output_file) - os.remove(output_file) + if save == False: + os.remove(output_file) + return weights +def xesmf_weights_to_xarray(regridder): + """ + Converts xESMF regridder weights to an xarray Dataset. + + This function takes a regridder object from xESMF, extracts the weights data, + and converts it into an xarray Dataset with relevant dimensions and attributes. + + Args: + regridder (xesmf.Regridder): A regridder object created using xESMF, which contains the weights to be converted. + + Returns: + xr.Dataset: An xarray Dataset containing the weights data with dimensions 'n_s', 'col', and 'row'. + """ + w = regridder.weights.data + dim = 'n_s' + ds = xr.Dataset( + { + 'S': (dim, w.data), + 'col': (dim, w.coords[1, :] + 1), + 'row': (dim, w.coords[0, :] + 1), + } + ) + ds.attrs = {'n_in': regridder.n_in, 'n_out': regridder.n_out} + return ds + +def reconstruct_xesmf_weights(ds_w): + """ + Reconstruct weights into a format that xESMF understands. + + This function takes a dataset with weights in a specific format and converts + it into a format that can be used by xESMF for regridding. + + Args: + ds_w (xarray.Dataset): The input dataset containing weights data. + It should have 'S', 'col', 'row', and appropriate attributes 'n_out' and 'n_in'. + + Returns: + xarray.DataArray: A DataArray containing reconstructed weights in COO format suitable for use with xESMF. + """ + col = ds_w['col'].values - 1 + row = ds_w['row'].values - 1 + s = ds_w['S'].values + n_out, n_in = ds_w.attrs['n_out'], ds_w.attrs['n_in'] + crds = np.stack([row, col]) + return xr.DataArray( + sparse.COO(crds, s, (n_out, n_in)), dims=('out_dim', 'in_dim'), name='weights' + ) + def parse_depths(depths, depths_from_file): """ Parses the selected depths from the available depth values and returns the corresponding depth indices and values. @@ -443,6 +500,7 @@ def save_data( lon, lat, out_path, + unstructured = False ): """ Saves the interpolated data to a NetCDF file. @@ -459,6 +517,7 @@ def save_data( lon (np.ndarray): The longitudes of the original dataset. lat (np.ndarray): The latitudes of the original dataset. out_path (str): The path to save the output NetCDF file. + unstructured (bool, optional): Interpolation to unstructured mesh. Default is False. Returns: None: This function does not return any value. @@ -476,16 +535,26 @@ def save_data( out_time = np.atleast_1d(shifted_time) else: out_time = np.atleast_1d(data.time.data[timesteps]) - - out1 = xr.Dataset( - {variable_name: (["time", "depth", "lat", "lon"], interpolated3d)}, + if unstructured is False: + out1 = xr.Dataset( + {variable_name: (["time", "depth", "lat", "lon"], interpolated3d)}, + coords={ + "time": out_time, + "depth": realdepths, + "lon": (["lon"], x), + "lat": (["lat"], y), + "longitude": (["lat", "lon"], lon), + "latitude": (["lat", "lon"], lat), + }, + attrs=data.attrs, + ) + else: + out1 = xr.Dataset( + {variable_name: (["time", "depth", "ncells"], interpolated3d)}, coords={ "time": out_time, "depth": realdepths, - "lon": (["lon"], x), - "lat": (["lat"], y), - "longitude": (["lat", "lon"], lon), - "latitude": (["lat", "lon"], lat), + }, attrs=data.attrs, ) @@ -504,18 +573,28 @@ def save_data( # ) # out1.to_netcdf(out_path, encoding={variable_name: {"zlib": True, "complevel": 9}}) - out1.to_netcdf( - out_path, - encoding={ - "time": {"dtype": np.dtype("double")}, - "depth": {"dtype": np.dtype("double")}, - "lat": {"dtype": np.dtype("double")}, - "lon": {"dtype": np.dtype("double")}, - "longitude": {"dtype": np.dtype("double")}, - "latitude": {"dtype": np.dtype("double")}, - variable_name: {"zlib": True, "complevel": 1, "dtype": np.dtype("single")}, - }, - ) + if unstructured is False: + out1.to_netcdf( + out_path, + encoding={ + "time": {"dtype": np.dtype("double")}, + "depth": {"dtype": np.dtype("double")}, + "lat": {"dtype": np.dtype("double")}, + "lon": {"dtype": np.dtype("double")}, + "longitude": {"dtype": np.dtype("double")}, + "latitude": {"dtype": np.dtype("double")}, + variable_name: {"zlib": True, "complevel": 1, "dtype": np.dtype("single")}, + }, + ) + else: + out1.to_netcdf( + out_path, + encoding={ + "time": {"dtype": np.dtype("double")}, + "depth": {"dtype": np.dtype("double")}, + variable_name: {"zlib": True, "complevel": 1, "dtype": np.dtype("single")}, + }, + ) # if args.rotate: # out2.to_netcdf( # out_path2, @@ -614,7 +693,7 @@ def fint(args=None): "--interp", choices=["nn", "mtri_linear", "linear_scipy", "cdo_remapcon","cdo_remaplaf","cdo_remapnn", "cdo_remapdis", - "smm_con","smm_laf","smm_nn","smm_dis"], # "idist", "linear", "cubic"], + "smm_remapcon","smm_remaplaf","smm_remapnn", "smm_remapdis", "xesmf_nearest_s2d"], # "idist", "linear", "cubic"], default="nn", help="Interpolation method. Options are \ nn - nearest neighbor (KDTree implementation, fast), \ @@ -678,6 +757,19 @@ def fint(args=None): Valid units are 'D' (days), 'h' (hours), 'm' (minutes), 's' (seconds). \ To substract timedelta, put argument in quotes, and prepend ' -', so SPACE and then -, e.g. ' -10D'.", ) + parser.add_argument( + "--weightspath", + type=str, + help="File with CDO weiths for smm interpolation.") + parser.add_argument( + "--save_weights", + action = "store_true", + help = "Save CDO weights ac netcdf file at the output directory.") + parser.add_argument( + "--to_fesom_mesh", + action = "store_true", + help = "Interpolate data to another fesom mesh. Works only with the target flag.") + args = parser.parse_args() @@ -782,12 +874,27 @@ def fint(args=None): # define region of interpolation if args.target is None: x, y, lon, lat = define_region(args.box, args.res, projection) + if args.to_fesom_mesh: + raise ValueError("target mesh for interpolation is not provided") else: - x, y, lon, lat = define_region_from_file(args.target) + if args.to_fesom_mesh: + lon, lat = load_mesh(args.target)[:2] + x = 0 + y = 0 + if placement == "elements": + endswith = "_elements_IFS.nc" + for root, dirs, files in os.walk(args.target): + for file in files: + if file.endswith(endswith): + mesh_xr = xr.open_dataset(os.path.join(root, file)) + lon = mesh_xr.grid_center_lon.values + lat = mesh_xr.grid_center_lat.values + else: + x, y, lon, lat = define_region_from_file(args.target) x2, y2 = match_longitude_format(x2, y2, lon, lat) # if we want to use shapelly mask, load it - if args.no_shape_mask is False: + if args.no_shape_mask is False and args.to_fesom_mesh is False: m2 = mask_ne(lon, lat) # additional variables, that we need for different interplations @@ -797,67 +904,99 @@ def fint(args=None): trifinder = triang2.get_trifinder() elif interpolation == "nn": distances, inds = create_indexes_and_distances(x2, y2, lon, lat, k=1, workers=4) - elif interpolation in ["cdo_remapcon", "cdo_remaplaf", "cdo_remapnn", "cdo_remapdis", "smm_con", "smm_laf", "smm_nn", "smm_dis"]: - gridtype = 'latlon' - gridsize = x.size*y.size - xsize = x.size - ysize = y.size - xname = 'longitude' - xlongname = 'longitude' - xunits = 'degrees_east' - yname = 'latitude' - ylongname = 'latitude' - yunits = 'degrees_north' - xfirst = float(lon[0,0]) - xinc = float(lon[0,1]-lon[0,0]) - yfirst = float(lat[0,0]) - yinc = float(lat[1,0]-lat[0,0]) - grid_mapping = [] - grid_mapping_name = [] - straight_vertical_longitude_from_pole = [] - latitude_of_projection_origin = [] - standard_parallel = [] - if projection == "np": - gridtype = 'projection' - xlongname = 'x coordinate of projection' - xunits = 'meters' - ylongname = 'y coordinate of projection' - yunits = 'meters' - xfirst = float(x[0]) - xinc = float(x[1]-x[0]) - yfirst = float(y[0]) - yinc = float(y[1]-y[0]) - grid_mapping = 'crs' - grid_mapping_name = 'polar_stereographic' - straight_vertical_longitude_from_pole = 0.0 - latitude_of_projection_origin = 90.0 - standard_parallel = 71.0 - - - formatted_content = f"""\ - gridtype = {gridtype} - gridsize = {gridsize} - xsize = {xsize} - ysize = {ysize} - xname = {xname} - xlongname = "{xlongname}" - xunits = "{xunits}" - yname = {yname} - ylongname = "{ylongname}" - yunits = "{yunits}" - xfirst = {xfirst} - xinc = {xinc} - yfirst = {yfirst} - yinc = {yinc} - grid_mapping = {grid_mapping} - grid_mapping_name = {grid_mapping_name} - straight_vertical_longitude_from_pole = {straight_vertical_longitude_from_pole} - latitude_of_projection_origin = {latitude_of_projection_origin} - standard_parallel = {standard_parallel}""" - - target_grid_path = out_path.replace(".nc", "target_grid.txt") - with open(target_grid_path, 'w') as file: - file.write(formatted_content) + elif interpolation in ['xesmf_nearest_s2d']: + ds_in = xr.open_dataset(args.data) + ds_in = ds_in.assign_coords(lat=('nod2',y2), lon=('nod2',x2)) + ds_in.lat.attrs = {'units': 'degrees', 'standard_name': 'latitude'} + ds_in.lon.attrs = {'units': 'degrees', 'standard_name': 'longitude'} + ds_out = xr.Dataset( + { + 'x': xr.DataArray(x, dims=['x']), + 'y': xr.DataArray(y, dims=['y']), + 'lat': xr.DataArray(lat, dims=['y', 'x']), + 'lon': xr.DataArray(lon, dims=['y', 'x']), + }) + if args.weightspath is not None: + xesmf_weights_path = args.weightspath + ds_w = xr.open_dataset(xesmf_weights_path) + wegiths_xesmf = reconstruct_xesmf_weights(ds_w) + regridder = xe.Regridder(ds_in,ds_out, method='nearest_s2d', weights=wegiths_xesmf,locstream_in=True) + else: + regridder = xe.Regridder(ds_in, ds_out, method='nearest_s2d', locstream_in=True) + if args.save_weights is True: + ds_w = xesmf_weights_to_xarray(regridder) + xesmf_weights_path = out_path.replace(".nc", "xesmf_weights.nc") + ds_w.to_netcdf(xesmf_weights_path) + elif interpolation in ["cdo_remapcon", "cdo_remaplaf", "cdo_remapnn", "cdo_remapdis", "smm_remapcon", "smm_remaplaf", "smm_remapnn", "smm_remapdis"]: + if args.to_fesom_mesh: + endswith = "_nodes_IFS.nc" + if placement == "elements": + endswith = "_elements_IFS.nc" + for root, dirs, files in os.walk(args.target): + for file in files: + if file.endswith(endswith): + target_grid_path = os.path.join(root, file) + else: + gridtype = 'latlon' + gridsize = x.size*y.size + xsize = x.size + ysize = y.size + xname = 'longitude' + xlongname = 'longitude' + xunits = 'degrees_east' + yname = 'latitude' + ylongname = 'latitude' + yunits = 'degrees_north' + xfirst = float(lon[0,0]) + xinc = float(lon[0,1]-lon[0,0]) + yfirst = float(lat[0,0]) + yinc = float(lat[1,0]-lat[0,0]) + grid_mapping = [] + grid_mapping_name = [] + straight_vertical_longitude_from_pole = [] + latitude_of_projection_origin = [] + standard_parallel = [] + if projection == "np": + gridtype = 'projection' + xlongname = 'x coordinate of projection' + xunits = 'meters' + ylongname = 'y coordinate of projection' + yunits = 'meters' + xfirst = float(x[0]) + xinc = float(x[1]-x[0]) + yfirst = float(y[0]) + yinc = float(y[1]-y[0]) + grid_mapping = 'crs' + grid_mapping_name = 'polar_stereographic' + straight_vertical_longitude_from_pole = 0.0 + latitude_of_projection_origin = 90.0 + standard_parallel = 71.0 + + + formatted_content = f"""\ + gridtype = {gridtype} + gridsize = {gridsize} + xsize = {xsize} + ysize = {ysize} + xname = {xname} + xlongname = "{xlongname}" + xunits = "{xunits}" + yname = {yname} + ylongname = "{ylongname}" + yunits = "{yunits}" + xfirst = {xfirst} + xinc = {xinc} + yfirst = {yfirst} + yinc = {yinc} + grid_mapping = {grid_mapping} + grid_mapping_name = {grid_mapping_name} + straight_vertical_longitude_from_pole = {straight_vertical_longitude_from_pole} + latitude_of_projection_origin = {latitude_of_projection_origin} + standard_parallel = {standard_parallel}""" + + target_grid_path = out_path.replace(".nc", "target_grid.txt") + with open(target_grid_path, 'w') as file: + file.write(formatted_content) endswith = "_nodes_IFS.nc" if placement == "elements": @@ -870,15 +1009,27 @@ def fint(args=None): # we will fill this array with interpolated values if not args.oneout: - interpolated3d = np.zeros((len(timesteps), len(realdepths), len(y), len(x))) - if args.rotate: - interpolated3d2 = np.zeros( - (len(timesteps), len(realdepths), len(y), len(x)) - ) + if args.to_fesom_mesh is False: + interpolated3d = np.zeros((len(timesteps), len(realdepths), len(y), len(x))) + if args.rotate: + interpolated3d2 = np.zeros( + (len(timesteps), len(realdepths), len(y), len(x)) + ) + else: + interpolated3d = np.zeros((len(timesteps), len(realdepths), len(lat))) + if args.rotate: + interpolated3d2 = np.zeros( + (len(timesteps), len(realdepths), len(lat)) + ) else: - interpolated3d = np.zeros((1, len(realdepths), len(y), len(x))) - if args.rotate: - interpolated3d2 = np.zeros((1, len(realdepths), len(y), len(x))) + if args.to_fesom_mesh is False: + interpolated3d = np.zeros((1, len(realdepths), len(y), len(x))) + if args.rotate: + interpolated3d2 = np.zeros((1, len(realdepths), len(y), len(x))) + else: + interpolated3d = np.zeros((1, len(realdepths), len(lat))) + if args.rotate: + interpolated3d2 = np.zeros((1, len(realdepths), len(lat))) # main loop for t_index, ttime in enumerate(timesteps): @@ -962,38 +1113,63 @@ def fint(args=None): variable_name: {"dtype": np.dtype("double")}, }, ) + if args.weightspath is not None: + weights_file_path = args.weightspath + else: + if t_index == 0 and d_index == 0: + weights_file_path = out_path.replace(".nc", "weighs_cdo.nc") + weights = generate_cdo_weights(target_grid_path, + gridfile, + input_file_path, + weights_file_path, + interpolation, + save = True) interpolated = interpolate_cdo(target_grid_path, gridfile, input_file_path, output_file_path, variable_name, interpolation, + weights_file_path, mask_zero=args.no_mask_zero ) os.remove(input_file_path) - elif interpolation in ["smm_con", "smm_laf", "smm_nn", "smm_dis"]: + elif interpolation in ["smm_remapcon", "smm_remaplaf", "smm_remapnn", "smm_remapdis"]: input_data = xr.Dataset({variable_name: (["nod2"], data_in)}) if args.rotate: input_data = xr.Dataset({variable_name: (["nod2"], data_in2)}) - input_file_path = args.data.replace(".nc","cdo_interpolation.nc") - input_data.to_netcdf(input_file_path,encoding={ - variable_name: {"dtype": np.dtype("double")}, - }, - ) - output_file_path = out_path.replace(".nc", "weighs_cdo.nc") - weights = generate_cdo_weights(target_grid_path, - gridfile, - input_file_path, - output_file_path, - interpolation) - os.remove(input_file_path) - interpolator = Regridder(weights=weights) + if t_index == 0 and d_index == 0: + input_file_path = args.data.replace(".nc","cdo_interpolation.nc") + input_data.to_netcdf(input_file_path,encoding={ + variable_name: {"dtype": np.dtype("double")}, + }, + ) + output_file_path = out_path.replace(".nc", "weighs_cdo.nc") + if args.weightspath is not None: + weights_file = args.weightspath + weights = xr.open_dataset(weights_file) + else: + weights = generate_cdo_weights(target_grid_path, + gridfile, + input_file_path, + output_file_path, + interpolation, + save =args.save_weights) + os.remove(input_file_path) + interpolator = Regridder(weights=weights) interpolated = interpolator.regrid(input_data) interpolated = interpolated[variable_name].values mask_zero=args.no_mask_zero if mask_zero: interpolated[interpolated == 0] = np.nan + + elif interpolation in ['xesmf_nearest_s2d']: + ds_in = xr.Dataset({variable_name: (["nod2"], data_in)}) + ds_in = ds_in.assign_coords(lat=('nod2',y2), lon=('nod2',x2)) + ds_in.lat.attrs = {'units': 'degrees', 'standard_name': 'latitude'} + ds_in.lon.attrs = {'units': 'degrees', 'standard_name': 'longitude'} + interpolated = regridder(ds_in)[variable_name].values # masking of the data if mask_file is not None: @@ -1002,19 +1178,29 @@ def fint(args=None): interpolated[mask] = np.nan if args.rotate: interpolated2[mask] = np.nan - elif args.no_shape_mask is False: + elif args.no_shape_mask is False and args.to_fesom_mesh is False: interpolated[m2] = np.nan if args.rotate: interpolated2[m2] = np.nan if args.oneout: - interpolated3d[0, d_index, :, :] = interpolated - if args.rotate: - interpolated3d2[0, d_index, :, :] = interpolated2 + if args.to_fesom_mesh is False: + interpolated3d[0, d_index, :, :] = interpolated + if args.rotate: + interpolated3d2[0, d_index, :, :] = interpolated2 + else: + interpolated3d[0, d_index, :] = interpolated + if args.rotate: + interpolated3d2[0, d_index, :] = interpolated2 else: - interpolated3d[t_index, d_index, :, :] = interpolated - if args.rotate: - interpolated3d2[t_index, d_index, :, :] = interpolated2 + if args.to_fesom_mesh is False: + interpolated3d[t_index, d_index, :, :] = interpolated + if args.rotate: + interpolated3d2[t_index, d_index, :, :] = interpolated2 + else: + interpolated3d[t_index, d_index, :] = interpolated + if args.rotate: + interpolated3d2[t_index, d_index, :] = interpolated2 if args.oneout: out_path_one = out_path.replace(".nc", f"_{str(t_index).zfill(10)}.nc") @@ -1030,6 +1216,7 @@ def fint(args=None): lon, lat, out_path_one, + args.to_fesom_mesh ) if args.rotate: out_path_one2 = out_path2.replace( @@ -1047,9 +1234,16 @@ def fint(args=None): lon, lat, out_path_one2, + args.to_fesom_mesh ) - if interpolation in ["cdo_remapcon","cdo_remaplaf","cdo_remapnn","cdo_remapdis","smm_con","smm_nn","smm_laf","smm_dis"]: - os.remove(target_grid_path) + if interpolation in ["cdo_remapcon","cdo_remaplaf","cdo_remapnn","cdo_remapdis"]: + if args.to_fesom_mesh is False: + os.remove(target_grid_path) + if args.save_weights is False and args.weightspath is None and weights_file_path is not None: + os.remove(weights_file_path) + elif interpolation in ["smm_remapcon","smm_remapnn","smm_remaplaf","smm_remapdis"]: + if args.to_fesom_mesh is False: + os.remove(target_grid_path) # save data (always 4D array) if not args.oneout: @@ -1065,6 +1259,7 @@ def fint(args=None): lon, lat, out_path, + args.to_fesom_mesh ) if args.rotate: save_data( @@ -1079,6 +1274,7 @@ def fint(args=None): lon, lat, out_path2, + args.to_fesom_mesh ) diff --git a/src/fint/ut.py b/src/fint/ut.py index e3baac4..a295fe8 100644 --- a/src/fint/ut.py +++ b/src/fint/ut.py @@ -42,9 +42,9 @@ def nodes_or_ements(data, variable_name, node_num, elem_num): """ - if data[variable_name].shape[-1] == node_num: + if node_num in data[variable_name].shape: return "nodes" - elif data[variable_name].shape[-1] == elem_num: + elif elem_num in data[variable_name].shape : return "elements" diff --git a/test/tests.sh b/test/tests.sh index 80cd73f..ae4b575 100755 --- a/test/tests.sh +++ b/test/tests.sh @@ -48,11 +48,24 @@ fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp cdo_rema fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b arctic --interp cdo_remapcon #smm_regrid -fint ${FILE} ${MESH} --influence 500000 -t 0 -d -1 --interp smm_con --no_shape_mask -fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp smm_laf -fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b arctic --interp smm_nn -fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 --interp smm_dis +fint ${FILE} ${MESH} ${INFL} -t 0 -d -1 --interp smm_remapcon --no_shape_mask +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp smm_remaplaf +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b arctic --interp smm_remapnn +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 --interp smm_remapdis +#xesmf_regrid +fint ${FILE} ${MESH} ${INFL} -t -1 -d -1 -b "-150, 150, -50, 70" --interp xesmf_nearest_s2d +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b arctic --interp xesmf_nearest_s2d +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b gulf --interp xesmf_nearest_s2d + +#saving weights and reuse it +fint ${FILE} ${MESH} ${INFL} -t 1:5 -d -1 -b "-150, 150, -50, 70" --interp smm_remapcon --save_weights +export WEIGHTS="--weightspath ./temp.fesom.1948_interpolated_-150_150_-50_70_2.5_6125.0_1_4weighs_cdo.nc" +fint ${FILE} ${MESH} ${INFL} -t 1:5 -d -1 -b "-150, 150, -50, 70" --interp cdo_remapcon ${WEIGHTS} + +fint ${FILE} ${MESH} ${INFL} -t 1:5 -d -1 -b "-150, 150, -50, 70" --interp xesmf_nearest_s2d --save_weights +export WEIGHTS="--weightspath ./temp.fesom.1948_interpolated_-150_150_-50_70_2.5_6125.0_1_4xesmf_weights.nc" +fint ${FILE} ${MESH} ${INFL} -t -1 -d -1 -b "-150, 150, -50, 70" --interp xesmf_nearest_s2d ${WEIGHTS} # create mask fint ${FILE} ${MESH} ${INFL} -t 0 -d -1 --interp mtri_linear -o mask.nc @@ -67,6 +80,25 @@ fint ${FILE} ${MESH} ${INFL} -t 0 -d 500:3000 --interp nn --target mask.nc # Don't apply shapely mask fint ${FILE} ${MESH} ${INFL} -t 0 -d 500:3000 --interp nn --target mask.nc --no_shape_mask +# interpolate to a fesom (unstructured) grid +if [ ! -d "./test/mesh/core2" ]; then + mkdir "./test/mesh/core2" + wget -O ./test/mesh/core2/core2_old_griddes_elements_IFS.nc https://gitlab.awi.de/fesom/core2_old/-/raw/master/core2_old_griddes_elements_IFS.nc + wget -O ./test/mesh/core2/core2_old_griddes_nodes_IFS.nc https://gitlab.awi.de/fesom/core2_old/-/raw/master/core2_old_griddes_nodes_IFS.nc + wget -O ./test/mesh/core2/nod2d.out https://gitlab.awi.de/fesom/core2_old/-/raw/master/nod2d.out + wget -O ./test/mesh/core2/elem2d.out https://gitlab.awi.de/fesom/core2_old/-/raw/master/elem2d.out +fi +export TARGET="./test/mesh/core2/" +fint ${FILE} ${MESH} ${INFL} -d -1 --target ${TARGET} --to_fesom_mesh +fint ${FILE} ${MESH} ${INFL} -t 0 -d 0:10 --target ${TARGET} --to_fesom_mesh --interp mtri_linear +fint ${FILE} ${MESH} ${INFL} -d 0:50 --target ${TARGET} --to_fesom_mesh --interp cdo_remapcon +fint ${FILE} ${MESH} ${INFL} -d 0:100 --target ${TARGET} --to_fesom_mesh --interp cdo_remaplaf + +export FILE="./test/data/v.fesom.1948.nc" +fint ${FILE} ${MESH} ${INFL} -d -1 --target ${TARGET} --to_fesom_mesh +fint ${FILE} ${MESH} ${INFL} -d 0:50 --target ${TARGET} --to_fesom_mesh --interp cdo_remapcon +fint ${FILE} ${MESH} ${INFL} -d 0:100 --target ${TARGET} --to_fesom_mesh --interp cdo_remaplaf + # Different variables FILE="./test/data/u.fesom.1948.nc" fint ${FILE} ${MESH} ${INFL} -t 0:3 -d -1