diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 8eaf68b..016eddd 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -40,7 +40,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [ ubuntu-latest, macos-latest ] + os: [ ubuntu-latest ] python: [ '3.9', '3.10', '3.11' ] steps: - uses: actions/checkout@main diff --git a/stormevents/nhc/const.py b/stormevents/nhc/const.py index 2dbab48..8c69a94 100644 --- a/stormevents/nhc/const.py +++ b/stormevents/nhc/const.py @@ -7,7 +7,9 @@ class RMWFillMethod(Enum): none = None persistent = auto() - regression_penny_2023 = auto() + regression_penny_2023_with_smoothing = auto() + regression_penny_2023 = regression_penny_2023_with_smoothing + regression_penny_2023_no_smoothing = auto() # Bias correction values for the Rmax forecast diff --git a/stormevents/nhc/track.py b/stormevents/nhc/track.py index 3ef4ccc..2b0d4fb 100644 --- a/stormevents/nhc/track.py +++ b/stormevents/nhc/track.py @@ -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"] @@ -1328,7 +1377,10 @@ def clamp(n, minn, maxn): "radius_of_maximum_winds" ] - elif rmw_fill == RMWFillMethod.regression_penny_2023: + elif ( + 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[ @@ -1384,52 +1436,8 @@ 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_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[ diff --git a/tests/test_nhc.py b/tests/test_nhc.py index 549cf3c..6d4fb3b 100644 --- a/tests/test_nhc.py +++ b/tests/test_nhc.py @@ -468,7 +468,7 @@ 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, @@ -476,13 +476,32 @@ def test_rmw_fill_method_regression_penny_2023(): advisories=["OFCL"], rmw_fill=RMWFillMethod.regression_penny_2023, ) - 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():