From 456f2a0d91edccefdab73726fdff1c3e6de3a9f0 Mon Sep 17 00:00:00 2001 From: Kolya Lettl Date: Thu, 11 May 2023 09:44:31 +0200 Subject: [PATCH 1/5] chore: interpolate 2D coordinates with interpn --- tobac/utils/general.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tobac/utils/general.py b/tobac/utils/general.py index 00aa8f7d..0cac32a0 100644 --- a/tobac/utils/general.py +++ b/tobac/utils/general.py @@ -27,7 +27,7 @@ def add_coordinates(t, variable_cube): """ - from scipy.interpolate import interp2d, interp1d + from scipy.interpolate import interp2d, interp1d, interpn logging.debug("start adding coordinates from cube") @@ -91,16 +91,16 @@ def add_coordinates(t, variable_cube): # interpolate 2D coordinates: elif variable_cube.coord(coord).ndim == 2: if variable_cube.coord_dims(coord) == (hdim_1, hdim_2): - f = interp2d(dimvec_2, dimvec_1, variable_cube.coord(coord).points) - coordinate_points = np.asarray( - [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])] - ) + points = (dimvec_1, dimvec_2) + values = variable_cube.coord(coord).points + xi = np.column_stack((t["hdim_1"], t["hdim_2"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (hdim_2, hdim_1): - f = interp2d(dimvec_1, dimvec_2, variable_cube.coord(coord).points) - coordinate_points = np.asarray( - [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])] - ) + points = (dimvec_2, dimvec_1) + values = variable_cube.coord(coord).points + xi = np.column_stack((t["hdim_2"], t["hdim_1"])) + coordinate_points = interpn(points, values, xi) # interpolate 3D coordinates: # mainly workaround for wrf latitude and longitude (to be fixed in future) From 122f50bdfd5ec5431ae60126ac4abda8ab7f0f41 Mon Sep 17 00:00:00 2001 From: Kolya Lettl Date: Thu, 11 May 2023 11:31:53 +0200 Subject: [PATCH 2/5] chore: interpolate 3D coordinates with interpn --- tobac/utils/general.py | 48 +++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/tobac/utils/general.py b/tobac/utils/general.py index 0cac32a0..565411dd 100644 --- a/tobac/utils/general.py +++ b/tobac/utils/general.py @@ -107,40 +107,40 @@ def add_coordinates(t, variable_cube): elif variable_cube.coord(coord).ndim == 3: if variable_cube.coord_dims(coord) == (ndim_time, hdim_1, hdim_2): - f = interp2d( - dimvec_2, dimvec_1, variable_cube[0, :, :].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])] + points = (dimvec_1, dimvec_2) + values = variable_cube[0, :, :].coord(coord).points + xi = np.column_stack((t["hdim_1"], t["hdim_2"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (ndim_time, hdim_2, hdim_1): - f = interp2d( - dimvec_1, dimvec_2, variable_cube[0, :, :].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])] + points = (dimvec_2, dimvec_1) + values = variable_cube[0, :, :].coord(coord).points + xi = np.column_stack((t["hdim_2"], t["hdim_1"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (hdim_1, ndim_time, hdim_2): - f = interp2d( - dimvec_2, dimvec_1, variable_cube[:, 0, :].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])] + points = (dimvec_1, dimvec_2) + values = variable_cube[:, 0, :].coord(coord).points + xi = np.column_stack((t["hdim_1"], t["hdim_2"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (hdim_1, hdim_2, ndim_time): - f = interp2d( - dimvec_2, dimvec_1, variable_cube[:, :, 0].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim1"])] + points = (dimvec_1, dimvec_2) + values = variable_cube[:, :, 0].coord(coord).points + xi = np.column_stack((t["hdim_1"], t["hdim_2"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (hdim_2, ndim_time, hdim_1): - f = interp2d( - dimvec_1, dimvec_2, variable_cube[:, 0, :].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])] + points = (dimvec_2, dimvec_1) + values = variable_cube[:, 0, :].coord(coord).points + xi = np.column_stack((t["hdim_2"], t["hdim_1"])) + coordinate_points = interpn(points, values, xi) if variable_cube.coord_dims(coord) == (hdim_2, hdim_1, ndim_time): - f = interp2d( - dimvec_1, dimvec_2, variable_cube[:, :, 0].coord(coord).points - ) - coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])] + points = (dimvec_2, dimvec_1) + values = variable_cube[:, :, 0].coord(coord).points + xi = np.column_stack((t["hdim_2"], t["hdim_1"])) + coordinate_points = interpn(points, values, xi) # write resulting array or list into DataFrame: t[coord] = coordinate_points From 3e1ba1a1a2373013ae319edfbb69504755742ac0 Mon Sep 17 00:00:00 2001 From: Kolya Lettl Date: Thu, 11 May 2023 11:33:53 +0200 Subject: [PATCH 3/5] chore: removed obsolete import of interp2d --- tobac/utils/general.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tobac/utils/general.py b/tobac/utils/general.py index 565411dd..fb319dcb 100644 --- a/tobac/utils/general.py +++ b/tobac/utils/general.py @@ -27,7 +27,7 @@ def add_coordinates(t, variable_cube): """ - from scipy.interpolate import interp2d, interp1d, interpn + from scipy.interpolate import interp1d, interpn logging.debug("start adding coordinates from cube") From 8af25cfe064c9d985280b2fb08de0b78ff761244 Mon Sep 17 00:00:00 2001 From: Kolya Lettl Date: Thu, 11 May 2023 12:03:18 +0200 Subject: [PATCH 4/5] chore: replaced interp2d with interpn in add_coordinates_3D --- tobac/utils/general.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tobac/utils/general.py b/tobac/utils/general.py index fb319dcb..bafe6c8c 100644 --- a/tobac/utils/general.py +++ b/tobac/utils/general.py @@ -258,10 +258,10 @@ def add_coordinates_3D( elif var_coord.ndim == 2: first_dim = coord_to_ax[variable_cube.coord_dims(coord)[1]] second_dim = coord_to_ax[variable_cube.coord_dims(coord)[0]] - f = interp2d(first_dim[0], second_dim[0], var_coord.points) - coordinate_points = [ - f(a, b) for a, b in zip(t[first_dim[1]], t[second_dim[1]]) - ] + points = (second_dim[0], first_dim[0]) + values = var_coord.points + xi = np.column_stack((t[second_dim[1]], t[first_dim[1]])) + coordinate_points = interpn(points, values, xi) # Deal with the special case where the coordinate is 3D but # one of the dimensions is time and we assume the coordinates @@ -276,10 +276,10 @@ def add_coordinates_3D( hdim2_pos = 1 if time_pos == 2 else 2 first_dim = coord_to_ax[variable_cube.coord_dims(coord)[hdim2_pos]] second_dim = coord_to_ax[variable_cube.coord_dims(coord)[hdim1_pos]] - f = interp2d(first_dim[0], second_dim[0], var_coord.points) - coordinate_points = [ - f(a, b) for a, b in zip(t[first_dim[1]], t[second_dim[1]]) - ] + points = (second_dim[0], second_dim[1]) + values = var_coord.points + xi = np.column_stack((t[second_dim[1]], t[first_dim[1]])) + coordinate_points = interpn(points, values, xi) # interpolate 3D coordinates: elif var_coord.ndim == 3: From 28eb1735215687c1c5b1b37efecdb84229843984 Mon Sep 17 00:00:00 2001 From: Kolya Lettl Date: Wed, 17 May 2023 12:54:05 +0200 Subject: [PATCH 5/5] fix: fixed incorrect definition of points --- tobac/utils/general.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tobac/utils/general.py b/tobac/utils/general.py index bafe6c8c..c542a875 100644 --- a/tobac/utils/general.py +++ b/tobac/utils/general.py @@ -276,7 +276,7 @@ def add_coordinates_3D( hdim2_pos = 1 if time_pos == 2 else 2 first_dim = coord_to_ax[variable_cube.coord_dims(coord)[hdim2_pos]] second_dim = coord_to_ax[variable_cube.coord_dims(coord)[hdim1_pos]] - points = (second_dim[0], second_dim[1]) + points = (second_dim[0], first_dim[0]) values = var_coord.points xi = np.column_stack((t[second_dim[1]], t[first_dim[1]])) coordinate_points = interpn(points, values, xi)