Skip to content

Commit

Permalink
Merge pull request #21 from FESOM/cdo_interpolation
Browse files Browse the repository at this point in the history 
Cdo interpolation
  • Loading branch information
@koldunovn
koldunovn authored Oct 2, 2023
2 parents af0023a + d0a89f9 commit e1e1047
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 8 deletions.
130 changes: 129 additions & 1 deletion src/fint/fint.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
import numpy as np
import pandas as pd
import xarray as xr
import subprocess
from scipy.interpolate import (
CloughTocher2DInterpolator,
LinearNDInterpolator,
Expand Down Expand Up @@ -264,6 +265,40 @@ 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):
"""
Interpolate a variable in a file using CDO (Climate Data Operators).
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.
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).
Returns:
np.ndarray: Interpolated variable data as a NumPy array.
"""

command = [
"cdo",
f"-{interpolation.split('_')[1]},{target_grid}",
f"-setgrid,{gridfile}",
f"{original_file}",
f"{output_file}"
]
if mask_zero:
command.insert(1, "-setctomiss,0")


# Execute the command
subprocess.run(command)

interpolated = xr.open_dataset(output_file)[variable_name].values
os.remove(output_file)
return interpolated

def parse_depths(depths, depths_from_file):
"""
Parses the selected depths from the available depth values and returns the corresponding depth indices and values.
Expand Down Expand Up @@ -549,7 +584,7 @@ def fint(args=None):
)
parser.add_argument(
"--interp",
choices=["nn", "mtri_linear", "linear_scipy"], # "idist", "linear", "cubic"],
choices=["nn", "mtri_linear", "linear_scipy", "cdo_remapcon","cdo_remaplaf","cdo_remapnn", "cdo_remapdis"], # "idist", "linear", "cubic"],
default="nn",
help="Interpolation method. Options are \
nn - nearest neighbor (KDTree implementation, fast), \
Expand Down Expand Up @@ -732,6 +767,76 @@ 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"]:
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":
endswith = "_elements_IFS.nc"
for root, dirs, files in os.walk(args.meshpath):
for file in files:
if file.endswith(endswith):
gridfile = os.path.join(root, file)


# we will fill this array with interpolated values
if not args.oneout:
Expand Down Expand Up @@ -816,6 +921,27 @@ def fint(args=None):
interpolated = interpolate_linear_scipy(data_in, x2, y2, lon, lat)
if args.rotate:
interpolated2 = interpolate_linear_scipy(data_in2, x2, y2, lon, lat)

elif interpolation in ["cdo_remapcon", "cdo_remaplaf", "cdo_remapnn", "cdo_remapdis"]:
input_data = xr.Dataset({variable_name: (["nod2"], data_in)})
if args.rotate:
input_data = xr.Dataset({variable_name: (["nod2"], data_in2)})

output_file_path = out_path.replace(".nc", "output_cdo_file.nc")
input_file_path = args.data.replace(".nc","cdo_interpolation.nc")
input_data.to_netcdf(input_file_path,encoding={
variable_name: {"dtype": np.dtype("double")},
},
)
interpolated = interpolate_cdo(target_grid_path,
gridfile,
input_file_path,
output_file_path,
variable_name,
interpolation,
mask_zero=args.no_mask_zero
)
os.remove(input_file_path)

# masking of the data
if mask_file is not None:
Expand Down Expand Up @@ -870,6 +996,8 @@ def fint(args=None):
lat,
out_path_one2,
)
if interpolation in ["cdo_remapcon", "cdo_remaplaf", "cdo_remapnn", "cdo_remapdis"]:
os.remove(target_grid_path)

# save data (always 4D array)
if not args.oneout:
Expand Down
18 changes: 12 additions & 6 deletions src/fint/regions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,8 @@ def define_region(box, res, projection=None):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 360, 170
res = (360,170)
lonNumber, latNumber = res

x = np.linspace(left, right, lonNumber)
y = np.linspace(down, up, latNumber)
Expand Down Expand Up @@ -135,7 +136,8 @@ def region_gs(res):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 300, 250
res = (300,250)
lonNumber, latNumber = res
left = -80
right = -30
bottom = 20
Expand Down Expand Up @@ -166,7 +168,8 @@ def region_trop(res):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 751, 400
res = (751,400)
lonNumber, latNumber = res
left = -60
right = 15
bottom = -9.95
Expand Down Expand Up @@ -200,7 +203,8 @@ def region_arctic(res):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 500, 500
res = (500,500)
lonNumber, latNumber = res
left = -180
right = 180
bottom = 60
Expand Down Expand Up @@ -235,7 +239,8 @@ def region_gulf(res):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 1000, 500
res = (1000, 500)
lonNumber, latNumber = res
left = -80
right = -30
bottom = 20
Expand Down Expand Up @@ -280,7 +285,8 @@ def region_cartopy(box, res, projection="mer"):
if not res is None:
lonNumber, latNumber = res
else:
lonNumber, latNumber = 500, 500
res = (500,500)
lonNumber, latNumber = res
left, right, down, up = list(map(float, box.split(",")))
print(left, right, down, up)
# left = -80
Expand Down
14 changes: 13 additions & 1 deletion test/tests.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@ set -o xtrace

wget https://swift.dkrz.de/v1/dkrz_c719fbc3-98ea-446c-8e01-356dac22ed90/fint/test_fint.tar
tar -xvf test_fint.tar

wget -O ./test/mesh/pi/pi_griddes_elements_IFS.nc https://gitlab.awi.de/fesom/pi/-/raw/master/pi_griddes_elements_IFS.nc
wget -O ./test/mesh/pi/pi_griddes_nodes_IFS.nc https://gitlab.awi.de/fesom/pi/-/raw/master/pi_griddes_nodes_IFS.nc
export FILE="./test/data/temp.fesom.1948.nc"
export MESH="./test/mesh/pi/"
export INFL="--influence 500000"
Expand Down Expand Up @@ -36,6 +37,17 @@ fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 --interp mtri_linear
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp mtri_linear
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b gulf --interp mtri_linear

#cdo_remapcon
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 --interp cdo_remapcon
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp cdo_remapcon
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b gulf --interp cdo_remapcon

#cdo_remaplaf
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 --interp cdo_remaplaf
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b "-150, 150, -50, 70" --interp cdo_remaplaf
fint ${FILE} ${MESH} ${INFL} -t 0 -d 0 -b arctic --interp cdo_remapcon


# create mask
fint ${FILE} ${MESH} ${INFL} -t 0 -d -1 --interp mtri_linear -o mask.nc
fint ${FILE} ${MESH} ${INFL} -t 0 -d -1 --interp mtri_linear --mask mask.nc
Expand Down

0 comments on commit e1e1047

Please sign in to comment.