Skip to content

Commit

Permalink
Merge pull request #434 from w-k-jones/fix_interp_fp_bug
Browse files Browse the repository at this point in the history 
Fix floating point error in 2d/3d coordinate interpolation
  • Loading branch information
@w-k-jones
w-k-jones authored Sep 11, 2024
2 parents dadad41 + 294ced5 commit 4fa1efa
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 37 deletions.
22 changes: 13 additions & 9 deletions tobac/feature_detection.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
"""

from __future__ import annotations
from typing import Union, Callable
from typing import Optional, Union, Callable
import warnings
import logging

Expand Down Expand Up @@ -199,8 +199,8 @@ def feature_position(

elif position_threshold == "weighted_diff":
# get position as centre of identified region, weighted by difference from the threshold:
weights = abs(track_data_region[region_small] - threshold_i)
if sum(weights) == 0:
weights = np.abs(track_data_region[region_small] - threshold_i)
if np.sum(weights) == 0:
weights = None
hdim1_weights = weights
hdim2_weights = weights
Expand All @@ -211,8 +211,8 @@ def feature_position(

elif position_threshold == "weighted_abs":
# get position as centre of identified region, weighted by absolute values if the field:
weights = abs(track_data_region[region_small])
if sum(weights) == 0:
weights = np.abs(track_data_region[region_small])
if np.sum(weights) == 0:
weights = None
hdim1_weights = weights
hdim2_weights = weights
Expand All @@ -230,14 +230,18 @@ def feature_position(
hdim1_index = pbc_utils.weighted_circmean(
hdim1_indices, weights=hdim1_weights, high=hdim1_max + 1, low=hdim1_min
)
hdim1_index = np.clip(hdim1_index, 0, hdim1_max + 1)
else:
hdim1_index = np.average(hdim1_indices, weights=hdim1_weights)
hdim1_index = np.clip(hdim1_index, 0, hdim1_max)
if PBC_flag in ("hdim_2", "both"):
hdim2_index = pbc_utils.weighted_circmean(
hdim2_indices, weights=hdim2_weights, high=hdim2_max + 1, low=hdim2_min
)
hdim2_index = np.clip(hdim2_index, 0, hdim2_max + 1)
else:
hdim2_index = np.average(hdim2_indices, weights=hdim2_weights)
hdim2_index = np.clip(hdim2_index, 0, hdim2_max)
if is_3D:
vdim_index = np.average(vdim_indices, weights=vdim_weights)

Expand Down Expand Up @@ -1147,10 +1151,10 @@ def feature_detection_multithreshold(
min_distance: float = 0,
feature_number_start: int = 1,
PBC_flag: Literal["none", "hdim_1", "hdim_2", "both"] = "none",
vertical_coord: str = None,
vertical_axis: int = None,
detect_subset: dict = None,
wavelength_filtering: tuple = None,
vertical_coord: Optional[str] = None,
vertical_axis: Optional[int] = None,
detect_subset: Optional[dict] = None,
wavelength_filtering: Optional[tuple] = None,
dz: Union[float, None] = None,
strict_thresholding: bool = False,
statistic: Union[dict[str, Union[Callable, tuple[Callable, dict]]], None] = None,
Expand Down
61 changes: 33 additions & 28 deletions tobac/utils/general.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@

import copy
import logging
from typing_extensions import Literal
import iris
import pandas as pd

from . import internal as internal_utils
Expand All @@ -16,13 +18,16 @@
import warnings


def add_coordinates(t, variable_cube):
def add_coordinates(
features: pd.DataFrame,
variable_cube: iris.cube.Cube,
):
"""Add coordinates from the input cube of the feature detection
to the trajectories/features.
Parameters
----------
t : pandas.DataFrame
features : pandas.DataFrame
Trajectories/features from feature detection or linking step.
variable_cube : iris.cube.Cube
Expand All @@ -32,7 +37,7 @@ def add_coordinates(t, variable_cube):
Returns
-------
t : pandas.DataFrame
pandas.DataFrame
Trajectories with added coordinates.
"""
Expand All @@ -42,17 +47,17 @@ def add_coordinates(t, variable_cube):
logging.debug("start adding coordinates from cube")

# pull time as datetime object and timestr from input data and add it to DataFrame:
t["time"] = None
t["timestr"] = None
features["time"] = None
features["timestr"] = None

logging.debug("adding time coordinate")

time_in = variable_cube.coord("time")
time_in_datetime = time_in.units.num2date(time_in.points)

t["time"] = time_in_datetime[t["frame"]]
t["timestr"] = [
x.strftime("%Y-%m-%d %H:%M:%S") for x in time_in_datetime[t["frame"]]
features["time"] = time_in_datetime[features["frame"]]
features["timestr"] = [
x.strftime("%Y-%m-%d %H:%M:%S") for x in time_in_datetime[features["frame"]]
]

# Get list of all coordinates in input cube except for time (already treated):
Expand Down Expand Up @@ -88,29 +93,29 @@ def add_coordinates(t, variable_cube):
variable_cube.coord(coord).points,
fill_value="extrapolate",
)
coordinate_points = f(t["hdim_1"])
coordinate_points = f(features["hdim_1"])

if variable_cube.coord_dims(coord) == (hdim_2,):
f = interp1d(
dimvec_2,
variable_cube.coord(coord).points,
fill_value="extrapolate",
)
coordinate_points = f(t["hdim_2"])
coordinate_points = f(features["hdim_2"])

# interpolate 2D coordinates:
elif variable_cube.coord(coord).ndim == 2:
if variable_cube.coord_dims(coord) == (hdim_1, hdim_2):
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)
xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (hdim_2, hdim_1):
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)
xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

# interpolate 3D coordinates:
# mainly workaround for wrf latitude and longitude (to be fixed in future)
Expand All @@ -119,44 +124,44 @@ def add_coordinates(t, variable_cube):
if variable_cube.coord_dims(coord) == (ndim_time, hdim_1, hdim_2):
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)
xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (ndim_time, hdim_2, hdim_1):
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)
xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (hdim_1, ndim_time, hdim_2):
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)
xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (hdim_1, hdim_2, ndim_time):
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)
xi = np.column_stack((features["hdim_1"], features["hdim_2"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (hdim_2, ndim_time, hdim_1):
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)
xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

if variable_cube.coord_dims(coord) == (hdim_2, hdim_1, ndim_time):
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)
xi = np.column_stack((features["hdim_2"], features["hdim_1"]))
coordinate_points = interpn(points, values, xi, bounds_error=False)

# write resulting array or list into DataFrame:
t[coord] = coordinate_points
features[coord] = coordinate_points

logging.debug("added coord: " + coord)
return t
return features


def add_coordinates_3D(
Expand Down

0 comments on commit 4fa1efa

Please sign in to comment.