Skip to content

Feature/rmax fill nosmooth #114

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
SorooshMani-NOAA merged 6 commits into from
Jan 27, 2025
Merged

Feature/rmax fill nosmooth #114

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
edd0771
changing format of hours and minutes in moving mean smoother call
WPringle Jan 22, 2025
d505a90
moving the moving mean smoother for the RMW fill to a function for cl…
WPringle Jan 22, 2025
daea5f7
adding option for the penny et al (2023) rmw regression fill without …
WPringle Jan 22, 2025
0b523a0
adding regression fill option with_smoothing that is same as the defa…
WPringle Jan 22, 2025
0f0d693
changing to make the penny_2023 flag actually just equal to the one e…
WPringle Jan 23, 2025
c6ca722
remove macos from test matrix
SorooshMani-NOAA Jan 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions stormevents/nhc/const.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ class RMWFillMethod(Enum):
none = None
persistent = auto()
regression_penny_2023 = auto()
regression_penny_2023_with_smoothing = auto()
Copy link
Collaborator

@SorooshMani-NOAA SorooshMani-NOAA Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@WPringle you don't need to define a new number, you can just make sure both enum names are the same, e.g.:

In [1]: import enum

In [2]: class E(enum.Enum):
   ...:     a = enum.auto()
   ...:     b = enum.auto()
   ...:     c = b
   ...:     d = enum.auto()
   ...:

In [3]: E
Out[3]: <enum 'E'>

In [4]: E.a
Out[4]: <E.a: 1>

In [5]: E.b
Out[5]: <E.b: 2>

In [6]: E.c
Out[6]: <E.b: 2>

In [7]: E.d
Out[7]: <E.d: 3>

Copy link
Collaborator

@SorooshMani-NOAA SorooshMani-NOAA Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually in this case we probably have to set penny_regression_2023 = penny_regression_2023_with_smoothing rather than the other way around, but I think it's still a good idea to keep the resulting enum value the same.

Copy link
Contributor Author

@WPringle WPringle Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was wondering about that, so you can just set c=b as you did in that class ok.

regression_penny_2023_no_smoothing = auto()


# Bias correction values for the Rmax forecast
Expand Down
106 changes: 59 additions & 47 deletions stormevents/nhc/track.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1287,6 +1287,55 @@ def correct_ofcl_based_on_carq_n_hollandb(
def clamp(n, minn, maxn):
return max(min(maxn, n), minn)

def movingmean(dff):
fcsthr_index = dff["forecast_hours"].drop_duplicates().index
df_temp = dff.loc[fcsthr_index].copy()
# make sure 60, 84, and 108 are added
fcsthrs_12hr = numpy.unique(
numpy.append(df_temp["forecast_hours"].values, [60, 84, 108])
)
rmw_12hr = numpy.interp(
fcsthrs_12hr,
df_temp["forecast_hours"],
df_temp["radius_of_maximum_winds"],
)
dt_12hr = pandas.to_datetime(
fcsthrs_12hr, unit="h", origin=df_temp["datetime"].iloc[0]
)
df_temp = DataFrame(
data={
"forecast_hours": fcsthrs_12hr,
"radius_of_maximum_winds": rmw_12hr,
},
index=dt_12hr,
)
rmw_rolling = df_temp.rolling(window="24h1min", center=True, min_periods=1)[
"radius_of_maximum_winds"
].mean()
for valid_time, rmw in rmw_rolling.items():
valid_index = dff["datetime"] == valid_time
if (
valid_index.sum() == 0
or dff.loc[valid_index, "forecast_hours"].iloc[0] == 0
):
continue
# make sure rolling rmw is not larger than the maximum radii of the strongest isotach
# this problem usually comes from the rolling average
max_isotach_radii = isotach_radii.loc[valid_index].iloc[-1].max()
if rmw < max_isotach_radii or numpy.isnan(max_isotach_radii):
dff.loc[valid_index, "radius_of_maximum_winds"] = rmw
# in case it does not come from rolling average just set to be Vr/Vmax ratio of max_isotach_radii
if (
dff.loc[valid_index, "radius_of_maximum_winds"].iloc[-1]
> max_isotach_radii
):
dff.loc[valid_index, "radius_of_maximum_winds"] = (
max_isotach_radii
* dff.loc[valid_index, "isotach_radius"].iloc[-1]
/ dff.loc[valid_index, "max_sustained_wind_speed"].iloc[-1]
)
return dff

ofcl_tracks = tracks["OFCL"]
carq_tracks = tracks["CARQ"]

Expand Down Expand Up @@ -1328,7 +1377,11 @@ def clamp(n, minn, maxn):
"radius_of_maximum_winds"
]

elif rmw_fill == RMWFillMethod.regression_penny_2023:
elif (
rmw_fill == RMWFillMethod.regression_penny_2023
or rmw_fill == RMWFillMethod.regression_penny_2023_with_smoothing
or rmw_fill == RMWFillMethod.regression_penny_2023_no_smoothing
):
# fill OFCL maximum wind radius based on regression method from
# Penny et al. (2023). https://doi.org/10.1175/WAF-D-22-0209.1
isotach_radii = forecast[
Expand Down Expand Up @@ -1384,52 +1437,11 @@ def clamp(n, minn, maxn):
rmw_, 5.0, max(120.0, rmw0)
)
# apply 24-HR moving mean to unique datetimes
fcsthr_index = forecast["forecast_hours"].drop_duplicates().index
df_temp = forecast.loc[fcsthr_index].copy()
# make sure 60, 84, and 108 are added
fcsthrs_12hr = numpy.unique(
numpy.append(df_temp["forecast_hours"].values, [60, 84, 108])
)
rmw_12hr = numpy.interp(
fcsthrs_12hr,
df_temp["forecast_hours"],
df_temp["radius_of_maximum_winds"],
)
dt_12hr = pandas.to_datetime(
fcsthrs_12hr, unit="h", origin=df_temp["datetime"].iloc[0]
)
df_temp = DataFrame(
data={
"forecast_hours": fcsthrs_12hr,
"radius_of_maximum_winds": rmw_12hr,
},
index=dt_12hr,
)
rmw_rolling = df_temp.rolling(window="24.01 h", center=True, min_periods=1)[
"radius_of_maximum_winds"
].mean()
for valid_time, rmw in rmw_rolling.items():
valid_index = forecast["datetime"] == valid_time
if (
valid_index.sum() == 0
or forecast.loc[valid_index, "forecast_hours"].iloc[0] == 0
):
continue
# make sure rolling rmw is not larger than the maximum radii of the strongest isotach
# this problem usually comes from the rolling average
max_isotach_radii = isotach_radii.loc[valid_index].iloc[-1].max()
if rmw < max_isotach_radii or numpy.isnan(max_isotach_radii):
forecast.loc[valid_index, "radius_of_maximum_winds"] = rmw
# in case it does not come from rolling average just set to be Vr/Vmax ratio of max_isotach_radii
if (
forecast.loc[valid_index, "radius_of_maximum_winds"].iloc[-1]
> max_isotach_radii
):
forecast.loc[valid_index, "radius_of_maximum_winds"] = (
max_isotach_radii
* forecast.loc[valid_index, "isotach_radius"].iloc[-1]
/ forecast.loc[valid_index, "max_sustained_wind_speed"].iloc[-1]
)
if (
rmw_fill == RMWFillMethod.regression_penny_2023
or rmw_fill == RMWFillMethod.regression_penny_2023_with_smoothing
):
forecast = movingmean(forecast)

# fill OFCL background pressure with the first entry from 0-hr CARQ background pressure (at sea level)
forecast.loc[radp_missing, "background_pressure"] = carq_ref[
Expand Down
25 changes: 21 additions & 4 deletions tests/test_nhc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -468,21 +468,38 @@ def test_rmw_fill_method_persistent():
assert rmw.unique() == 10


def test_rmw_fill_method_regression_penny_2023():
def test_rmw_fill_method_regression_penny_2023_with_smoothing():
tr_florence2018 = VortexTrack.from_storm_name(
"Florence",
2018,
file_deck="a",
advisories=["OFCL"],
rmw_fill=RMWFillMethod.regression_penny_2023,
rmw_fill=RMWFillMethod.regression_penny_2023_with_smoothing,
)
assert tr_florence2018.rmw_fill == RMWFillMethod.regression_penny_2023
assert tr_florence2018.rmw_fill == RMWFillMethod.regression_penny_2023_with_smoothing
data = tr_florence2018.data
i_uq_row = 40
rmw = data.loc[data.track_start_time == data.track_start_time.unique()[i_uq_row]][
"radius_of_maximum_winds"
]
assert len(rmw.unique()) == 11


def test_rmw_fill_method_regression_penny_2023_no_smoothing():
tr_florence2018 = VortexTrack.from_storm_name(
"Florence",
2018,
file_deck="a",
advisories=["OFCL"],
rmw_fill=RMWFillMethod.regression_penny_2023_no_smoothing,
)
assert tr_florence2018.rmw_fill == RMWFillMethod.regression_penny_2023_no_smoothing
data = tr_florence2018.data
i_uq_row = 40
rmw = data.loc[data.track_start_time == data.track_start_time.unique()[i_uq_row]][
"radius_of_maximum_winds"
]
assert len(rmw.unique()) > 1
assert len(rmw.unique()) == 10


def test_rmw_fill_method_setvalue_invalid():
Expand Down
Loading