Skip to content

Commit

Permalink
Merge pull request #131 from zxdawn/interp_kind
Browse files Browse the repository at this point in the history 
Support different interpolation kinds
  • Loading branch information
@freemansw1
freemansw1 committed May 13, 2022
2 parents 4b9c798 + dc2a655 commit dc9166c
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 10 deletions.
10 changes: 9 additions & 1 deletion tobac/themes/tobac_v1/feature_detection.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ def feature_detection_multithreshold(
min_num=0,
target="maximum",
position_threshold="center",
coord_interp_kind='linear',
sigma_threshold=0.5,
n_erosion_threshold=0,
n_min_threshold=0,
Expand Down Expand Up @@ -63,6 +64,13 @@ def feature_detection_multithreshold(
Flag choosing method used for the position of the tracked
feature. Default is 'center'.
coord_interp_kind : str, optional
The kind of interpolation for coordinates. Default is 'linear'.
For 1d interp, {'linear', 'nearest', 'nearest-up', 'zero',
'slinear', 'quadratic', 'cubic',
'previous', 'next'}.
For 2d interp, {'linear', 'cubic', 'quintic'}.
sigma_threshold: float, optional
Standard deviation for intial filtering step. Default is 0.5.
Expand Down Expand Up @@ -146,7 +154,7 @@ def feature_detection_multithreshold(
features["feature"] = features.index + feature_number_start
# features_filtered = features.drop(features[features['num'] < min_num].index)
# features_filtered.drop(columns=['idx','num','threshold_value'],inplace=True)
features = add_coordinates(features, field_iris)
features = add_coordinates(features, field_iris, coord_interp_kind)
# convert pandas DataFrame to xarray
list_coords = [key for key in field_in.coords.keys()]
print(features)
Expand Down
48 changes: 39 additions & 9 deletions tobac/utils/general.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@


@xarray_to_iris
def add_coordinates(t, variable_cube):
def add_coordinates(t, variable_cube, coord_interp_kind):
"""Add coordinates from the tracking cube to the trajectories.
Coordinates: time, longitude&latitude, x&y dimensions.
Expand All @@ -20,6 +20,9 @@ def add_coordinates(t, variable_cube):
variable_cube : iris.cube.Cube
Input data used for the tracking to transfer coodinate information to resulting DataFrame
coord_interp_kind: str
The kind of interpolation for coordinates.
Returns
-------
t : pandas.DataFrame
Expand Down Expand Up @@ -73,11 +76,21 @@ def add_coordinates(t, variable_cube):
logging.debug("adding coord: " + coord)
# interpolate 2D coordinates:
if variable_cube.coord(coord).ndim == 1:
# check the interpolation method is valid
if coord_interp_kind not in ["linear", "nearest", "nearest-up",
"zero", "slinear", "quadratic",
"cubic", "previous", "next"]:
raise ValueError(
"coord_interp_kind must be linear, nearest, nearest-up, \
zero, slinear, quadratic, \
cubic, previous, or next for 1D coord"
)

if variable_cube.coord_dims(coord) == (hdim_1,):
f = interp1d(
dimvec_1,
variable_cube.coord(coord).points,
kind=coord_interp_kind,
fill_value="extrapolate",
)
coordinate_points = f(t["hdim_1"])
Expand All @@ -86,59 +99,76 @@ def add_coordinates(t, variable_cube):
f = interp1d(
dimvec_2,
variable_cube.coord(coord).points,
kind=coord_interp_kind,
fill_value="extrapolate",
)
coordinate_points = f(t["hdim_2"])

# interpolate 2D coordinates:
elif variable_cube.coord(coord).ndim == 2:
# check the interpolation method is valid
if coord_interp_kind not in ["linear", "cubic", "quintic"]:
raise ValueError(
"coord_interp_kind must be linear, cubic, or quintic for 2D coord"
)

if variable_cube.coord_dims(coord) == (hdim_1, hdim_2):
f = interp2d(dimvec_2, dimvec_1, variable_cube.coord(coord).points)
f = interp2d(dimvec_2, dimvec_1, variable_cube.coord(coord).points, kind=coord_interp_kind)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])]

if variable_cube.coord_dims(coord) == (hdim_2, hdim_1):
f = interp2d(dimvec_1, dimvec_2, variable_cube.coord(coord).points)
f = interp2d(dimvec_1, dimvec_2, variable_cube.coord(coord).points, kind=coord_interp_kind)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])]

# interpolate 3D coordinates:
# mainly workaround for wrf latitude and longitude (to be fixed in future)

elif variable_cube.coord(coord).ndim == 3:
# check the interpolation method is valid
if coord_interp_kind not in ["linear", "cubic", "quintic"]:
raise ValueError(
"coord_interp_kind must be linear, cubic, or quintic for 3D coord"
)

if variable_cube.coord_dims(coord) == (ndim_time, hdim_1, hdim_2):
f = interp2d(
dimvec_2, dimvec_1, variable_cube[0, :, :].coord(coord).points
dimvec_2, dimvec_1, variable_cube[0, :, :].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])]

if variable_cube.coord_dims(coord) == (ndim_time, hdim_2, hdim_1):
f = interp2d(
dimvec_1, dimvec_2, variable_cube[0, :, :].coord(coord).points
dimvec_1, dimvec_2, variable_cube[0, :, :].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])]

if variable_cube.coord_dims(coord) == (hdim_1, ndim_time, hdim_2):
f = interp2d(
dimvec_2, dimvec_1, variable_cube[:, 0, :].coord(coord).points
dimvec_2, dimvec_1, variable_cube[:, 0, :].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim_1"])]

if variable_cube.coord_dims(coord) == (hdim_1, hdim_2, ndim_time):
f = interp2d(
dimvec_2, dimvec_1, variable_cube[:, :, 0].coord(coord).points
dimvec_2, dimvec_1, variable_cube[:, :, 0].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_2"], t["hdim1"])]

if variable_cube.coord_dims(coord) == (hdim_2, ndim_time, hdim_1):
f = interp2d(
dimvec_1, dimvec_2, variable_cube[:, 0, :].coord(coord).points
dimvec_1, dimvec_2, variable_cube[:, 0, :].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])]

if variable_cube.coord_dims(coord) == (hdim_2, hdim_1, ndim_time):
f = interp2d(
dimvec_1, dimvec_2, variable_cube[:, :, 0].coord(coord).points
dimvec_1, dimvec_2, variable_cube[:, :, 0].coord(coord).points,
kind=coord_interp_kind
)
coordinate_points = [f(a, b) for a, b in zip(t["hdim_1"], t["hdim_2"])]

Expand Down

0 comments on commit dc9166c

Please sign in to comment.