From 4267b6b678c022a1948248b15f7f1fb6f621578a Mon Sep 17 00:00:00 2001
From: Boris Shapkin <shap99boris@gmail.com>
Date: Mon, 9 Oct 2023 14:44:25 +0200
Subject: [PATCH 1/7] add reusable weights

---
 src/fint/fint.py | 112 +++++++++++++++++++++++++++++++----------------
 test/tests.sh    |  13 ++++--
 2 files changed, 84 insertions(+), 41 deletions(-)

diff --git a/src/fint/fint.py b/src/fint/fint.py
index 64c8a55..a8e9a34 100644
--- a/src/fint/fint.py
+++ b/src/fint/fint.py
@@ -266,25 +266,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,32 +300,38 @@ 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 parse_depths(depths, depths_from_file):
@@ -614,7 +621,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"],  # "idist", "linear", "cubic"],
         default="nn",
         help="Interpolation method. Options are \
             nn - nearest neighbor (KDTree implementation, fast), \
@@ -678,6 +685,15 @@ 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.")
+    
 
     args = parser.parse_args()
 
@@ -797,7 +813,7 @@ 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"]:
+    elif interpolation in ["cdo_remapcon", "cdo_remaplaf", "cdo_remapnn", "cdo_remapdis", "smm_remapcon", "smm_remaplaf", "smm_remapnn", "smm_remapdis"]:
         gridtype = 'latlon'
         gridsize = x.size*y.size
         xsize = x.size
@@ -962,33 +978,51 @@ 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
@@ -1048,7 +1082,11 @@ def fint(args=None):
                     lat,
                     out_path_one2,
                 )
-    if interpolation in ["cdo_remapcon","cdo_remaplaf","cdo_remapnn","cdo_remapdis","smm_con","smm_nn","smm_laf","smm_dis"]:    
+    if interpolation in ["cdo_remapcon","cdo_remaplaf","cdo_remapnn","cdo_remapdis"]:    
+        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"]:
         os.remove(target_grid_path)
 
     # save data (always 4D array)
diff --git a/test/tests.sh b/test/tests.sh
index 80cd73f..eb6c69d 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -48,10 +48,15 @@ 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
+
+#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}
 
 
 # create mask

From e155735523fa5a89101abfaca041b5ce52e8ae1a Mon Sep 17 00:00:00 2001
From: Boris Shapkin <shap99boris@gmail.com>
Date: Fri, 27 Oct 2023 22:00:41 +0200
Subject: [PATCH 2/7] xesmf nearest_s2d regridding

---
 environment.yml  |  2 ++
 src/fint/fint.py | 82 +++++++++++++++++++++++++++++++++++++++++++++++-
 test/tests.sh    |  8 +++++
 3 files changed, 91 insertions(+), 1 deletion(-)

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 a8e9a34..5c039f6 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 (
@@ -334,6 +336,54 @@ def generate_cdo_weights(target_grid,gridfile,original_file,output_file,interpol
 
     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.
@@ -621,7 +671,7 @@ def fint(args=None):
         "--interp",
         choices=["nn", "mtri_linear", "linear_scipy",
                  "cdo_remapcon","cdo_remaplaf","cdo_remapnn", "cdo_remapdis",
-                 "smm_remapcon","smm_remaplaf","smm_remapnn", "smm_remapdis"],  # "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), \
@@ -813,6 +863,29 @@ 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 ['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"]:
         gridtype = 'latlon'
         gridsize = x.size*y.size
@@ -1028,6 +1101,13 @@ def fint(args=None):
                 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:
diff --git a/test/tests.sh b/test/tests.sh
index eb6c69d..3ef9a22 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -53,11 +53,19 @@ fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp smm_rema
 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

From 628512aa03417378e65c7949f7e5a4f960f133af Mon Sep 17 00:00:00 2001
From: Boris Shapkin <shap99boris@gmail.com>
Date: Thu, 21 Dec 2023 19:13:30 +0100
Subject: [PATCH 3/7] Interpolation to unstructured meshes

---
 src/fint/fint.py | 276 ++++++++++++++++++++++++++++++-----------------
 src/fint/ut.py   |   4 +-
 2 files changed, 179 insertions(+), 101 deletions(-)

diff --git a/src/fint/fint.py b/src/fint/fint.py
index 5c039f6..8f5a4d4 100644
--- a/src/fint/fint.py
+++ b/src/fint/fint.py
@@ -500,6 +500,7 @@ def save_data(
     lon,
     lat,
     out_path,
+    unstructured = False 
 ):
     """
     Saves the interpolated data to a NetCDF file.
@@ -516,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.
@@ -533,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,
     )
@@ -561,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,
@@ -743,6 +765,10 @@ def fint(args=None):
         "--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()
@@ -848,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
@@ -887,66 +928,75 @@ def fint(args=None):
             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"]:
-        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)
+        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":
@@ -959,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):
@@ -1116,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")
@@ -1144,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(
@@ -1161,13 +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"]:    
-        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"]:
-        os.remove(target_grid_path)
+        if args.to_fesom_mesh is False:
+            os.remove(target_grid_path)
 
     # save data (always 4D array)
     if not args.oneout:
@@ -1183,6 +1259,7 @@ def fint(args=None):
             lon,
             lat,
             out_path,
+            args.to_fesom_mesh
         )
         if args.rotate:
             save_data(
@@ -1197,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"
 
 

From 870f6e949945c1f393d9d39bc2dda82a6761c52f Mon Sep 17 00:00:00 2001
From: boryasbora <80640421+boryasbora@users.noreply.github.com>
Date: Fri, 22 Dec 2023 11:11:40 +0100
Subject: [PATCH 4/7] added test for unstructured meshes

---
 test/tests.sh | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/test/tests.sh b/test/tests.sh
index 3ef9a22..70709ce 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -80,6 +80,18 @@ 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
+export TARGET="/work/ab0995/a270226/swift.dkrz.de/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

From 155289cdfbc32f0bd454edf7c3baa68beaa499a2 Mon Sep 17 00:00:00 2001
From: boryasbora <80640421+boryasbora@users.noreply.github.com>
Date: Fri, 22 Dec 2023 11:27:40 +0100
Subject: [PATCH 5/7] download core2 grid files

---
 test/tests.sh | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/tests.sh b/test/tests.sh
index 70709ce..3b250c2 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -81,7 +81,12 @@ fint ${FILE} ${MESH} ${INFL} -t 0 -d 500:3000  --interp nn --target mask.nc
 fint ${FILE} ${MESH} ${INFL} -t 0 -d 500:3000  --interp nn --target mask.nc --no_shape_mask
 
 # interpolate to a fesom (unstructured) grid
-export TARGET="/work/ab0995/a270226/swift.dkrz.de/core2/"
+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
+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

From 0ed49c7f928fa0c8335c97eb86eb2aa96ca2d506 Mon Sep 17 00:00:00 2001
From: boryasbora <80640421+boryasbora@users.noreply.github.com>
Date: Fri, 22 Dec 2023 11:33:57 +0100
Subject: [PATCH 6/7] fix typo

---
 test/tests.sh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/tests.sh b/test/tests.sh
index 3b250c2..873aba5 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -86,7 +86,7 @@ if [ ! -d "./test/mesh/core2" ]; then
     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
 fi
-export TARGET="/test/mesh/core2/"
+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

From 8b2af9b4e2e15a8cebd1a203ece082f8fb62ab2b Mon Sep 17 00:00:00 2001
From: boryasbora <80640421+boryasbora@users.noreply.github.com>
Date: Fri, 22 Dec 2023 11:41:59 +0100
Subject: [PATCH 7/7] more files for core2 mesh

---
 test/tests.sh | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/tests.sh b/test/tests.sh
index 873aba5..ae4b575 100755
--- a/test/tests.sh
+++ b/test/tests.sh
@@ -85,6 +85,8 @@ 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