Skip to content

Commit

Permalink
Multiple bugfixes for CARQ correction of OFCL track (#68)
Browse files Browse the repository at this point in the history 
* Fix multiple bugs

* Corrections for adding forecast background pressure using CARQ background pressure, correcting comments

* Bugfix for background and central pressure correction on write

* Test for carq fix of ofcl track

* Update outdated test references

---------

Co-authored-by: William Pringle <wpringle@anl.gov>
  • Loading branch information
@SorooshMani-NOAA @WPringle
SorooshMani-NOAA and WPringle authored Mar 24, 2023
1 parent 2895957 commit a3e9c47
Show file tree
Hide file tree
Showing 16 changed files with 21,873 additions and 21,822 deletions.
164 changes: 100 additions & 64 deletions stormevents/nhc/track.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -564,7 +564,6 @@ def atcf(self, advisory: ATCF_Advisory = None) -> DataFrame:
atcf["max_sustained_wind_speed"] = (
atcf["max_sustained_wind_speed"].astype("string").str.pad(5)
)
atcf["central_pressure"] = atcf["central_pressure"].astype("string").str.pad(5)
atcf["development_level"] = atcf["development_level"].str.pad(3)
atcf["isotach_radius"] = atcf["isotach_radius"].astype("string").str.pad(4)
atcf["isotach_quadrant_code"] = atcf["isotach_quadrant_code"].str.pad(4)
Expand All @@ -582,22 +581,22 @@ def atcf(self, advisory: ATCF_Advisory = None) -> DataFrame:
)

atcf["background_pressure"].fillna(method="ffill", inplace=True)
atcf.loc[
~pandas.isna(self.data["central_pressure"])
& (self.data["background_pressure"] <= self.data["central_pressure"])
& (self.data["central_pressure"] < 1013),
"background_pressure",
] = "1013"
atcf.loc[
~pandas.isna(self.data["central_pressure"])
& (self.data["background_pressure"] <= self.data["central_pressure"])
& (self.data["central_pressure"] < 1013),
"background_pressure",
] = (
self.data["central_pressure"] + 1
atcf["background_pressure"] = atcf["background_pressure"].astype(int)
atcf["central_pressure"] = atcf["central_pressure"].astype(int)

press_cond_nobg = ~atcf["central_pressure"].isna() & (
(atcf["background_pressure"] <= atcf["central_pressure"])
| (atcf["background_pressure"].isna())
)
atcf.loc[press_cond_nobg, "background_pressure"] = 1013

press_cond_nobg_hieye = press_cond_nobg & (atcf["central_pressure"] >= 1013)
atcf.loc[press_cond_nobg_hieye, "background_pressure"] = (
atcf.loc[press_cond_nobg_hieye, "central_pressure"] + 1
)
atcf["central_pressure"] = atcf["central_pressure"].astype("string").str.pad(5)
atcf["background_pressure"] = (
atcf["background_pressure"].astype(int).astype("string").str.pad(5)
atcf["background_pressure"].astype("string").str.pad(5)
)

atcf["radius_of_last_closed_isobar"] = (
Expand Down Expand Up @@ -981,55 +980,9 @@ def unfiltered_data(self, dataframe: DataFrame):
# fill missing values of MRD and MSLP in the OFCL advisory
if "OFCL" in self.advisories:
tracks = separate_tracks(dataframe)

if "OFCL" in tracks:
ofcl_tracks = tracks["OFCL"]
carq_tracks = tracks["CARQ"]

for initial_time, forecast in ofcl_tracks.items():
if initial_time in carq_tracks:
carq_forecast = carq_tracks[initial_time]
else:
carq_forecast = carq_tracks[list(carq_tracks)[0]]

relation = HollandBRelation()
holland_b = relation.holland_b(
max_sustained_wind_speed=carq_forecast[
"max_sustained_wind_speed"
],
background_pressure=carq_forecast["background_pressure"],
central_pressure=carq_forecast["central_pressure"],
)

holland_b[holland_b == numpy.inf] = numpy.nan
holland_b = numpy.nanmean(holland_b)

mrd_missing = pandas.isna(forecast["radius_of_maximum_winds"])
mslp_missing = pandas.isna(forecast["central_pressure"])
radp_missing = pandas.isna(forecast["background_pressure"])

# fill OFCL maximum wind radius with the first entry from the CARQ advisory
forecast.loc[
mrd_missing, "radius_of_maximum_winds"
] = carq_forecast["radius_of_maximum_winds"].iloc[0]

# fill OFCL background pressure with the first entry from the CARQ advisory central pressure (at sea level)
forecast.loc[radp_missing, "background_pressure"] = carq_forecast[
"central_pressure"
].iloc[0]

# fill OFCL central pressure (at sea level) with the 3rd hour entry, preserving Holland B
forecast.loc[
mslp_missing, "central_pressure"
] = relation.central_pressure(
max_sustained_wind_speed=forecast.loc[
mslp_missing, "max_sustained_wind_speed"
],
background_pressure=forecast.loc[
mslp_missing, "background_pressure"
],
holland_b=holland_b,
)
if all(adv in tracks for adv in ["OFCL", "CARQ"]):
tracks = correct_ofcl_based_on_carq_n_hollandb(tracks)
dataframe = combine_tracks(tracks)

if len(self.__advisories_to_remove) > 0:
dataframe = dataframe[
Expand Down Expand Up @@ -1175,6 +1128,7 @@ def separate_tracks(data: DataFrame) -> Dict[str, Dict[str, DataFrame]]:
tracks = {}
for advisory in pandas.unique(data["advisory"]):
advisory_data = data[data["advisory"] == advisory]
advisory_data["forecast_hours"] = advisory_data.forecast_hours.astype(int)

if advisory == "BEST":
advisory_data = advisory_data.sort_values("datetime")
Expand All @@ -1196,3 +1150,85 @@ def separate_tracks(data: DataFrame) -> Dict[str, Dict[str, DataFrame]]:
] = track_data

return tracks


def combine_tracks(tracks: Dict[str, Dict[str, DataFrame]]) -> DataFrame:
"""
combine tracks separated using `separate_tracks`
:param tracks: dictionary of forecasts for each advisory (aside from best track ``BEST``, which only has one hindcast)
:return: data frame of track
"""

return pandas.concat([df for adv_trk in tracks.values() for df in adv_trk.values()])


def correct_ofcl_based_on_carq_n_hollandb(
tracks: Dict[str, Dict[str, DataFrame]]
) -> Dict[str, Dict[str, DataFrame]]:
"""
Correct official forecast using consensus track along with holland-b
relation
:param tracks: dictionary of forecasts for each advisory (aside from best track ``BEST``, which only has one hindcast)
:return: dictionary of forecasts for each advisory (aside from best track ``BEST``, which only has one hindcast) with corrected OFCL
"""

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

corr_ofcl_tracks = dict()

for initial_time, forecast in ofcl_tracks.items():
if initial_time in carq_tracks:
carq_forecast = carq_tracks[initial_time]
else:
carq_forecast = carq_tracks[list(carq_tracks)[0]]

relation = HollandBRelation()
holland_b = relation.holland_b(
max_sustained_wind_speed=carq_forecast["max_sustained_wind_speed"],
background_pressure=carq_forecast["background_pressure"],
central_pressure=carq_forecast["central_pressure"],
)

holland_b[holland_b == numpy.inf] = numpy.nan
holland_b = numpy.nanmean(holland_b)

# Get CARQ from forecast hour 0 and isotach 34kt (i.e. the first item)
carq_ref = carq_forecast.loc[carq_forecast.forecast_hours == 0].iloc[0]

columns_of_interest = forecast[
["radius_of_maximum_winds", "central_pressure", "background_pressure"]
]
columns_of_interest[columns_of_interest == 0] = pandas.NA
missing = columns_of_interest.isna()
# Order of columns is the same as columns_of_interest
mrd_missing = missing.iloc[:, 0]
mslp_missing = missing.iloc[:, 1]
radp_missing = missing.iloc[:, 2]

# fill OFCL maximum wind radius with the first entry from 0-hr CARQ
forecast.loc[mrd_missing, "radius_of_maximum_winds"] = carq_ref[
"radius_of_maximum_winds"
]

# 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[
"background_pressure"
]

# fill OFCL central pressure (at sea level), preserving Holland B from 0-hr CARQ
forecast.loc[mslp_missing, "central_pressure"] = relation.central_pressure(
max_sustained_wind_speed=forecast.loc[
mslp_missing, "max_sustained_wind_speed"
],
background_pressure=forecast.loc[mslp_missing, "background_pressure"],
holland_b=holland_b,
)

corr_ofcl_tracks[initial_time] = forecast

tracks["OFCL"] = corr_ofcl_tracks

return tracks
4 changes: 2 additions & 2 deletions tests/data/reference/test_vortex_track/harvey2017.fort.22
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
AL, 09, 2017081606, , BEST, 0, 137N,458W, 25, 1013, LO, 0, , 0, 0, 0, 0, 1013, 150, 80, 0, 0, ,,, 0, 0,, 1
AL, 09, 2017081606, , BEST, 0, 137N,458W, 25, 1013, LO, 0, , 0, 0, 0, 0, 1014, 150, 80, 0, 0, ,,, 0, 0,, 1
AL, 09, 2017081612, , BEST, 0, 137N,474W, 25, 1010, LO, 0, , 0, 0, 0, 0, 1013, 150, 80, 0, 0, L, 0, ,270, 8, INVEST, 2
AL, 09, 2017081618, , BEST, 0, 136N,490W, 25, 1009, LO, 0, , 0, 0, 0, 0, 1013, 150, 80, 0, 0, L, 0, ,267, 8, INVEST, 3
AL, 09, 2017081700, , BEST, 0, 136N,506W, 25, 1010, LO, 0, , 0, 0, 0, 0, 1012, 120, 80, 0, 0, L, 0, ,270, 8, INVEST, 4
Expand Down Expand Up @@ -92,4 +92,4 @@ AL, 09, 2017090112, , BEST, 0, 360N,871W, 25, 1002, EX, 0, , 0,
AL, 09, 2017090118, , BEST, 0, 365N,864W, 20, 1004, EX, 0, , 0, 0, 0, 0, 1010, 220, 80, 0, 0, , 0, , 48, 4, HARVEY, 71
AL, 09, 2017090200, , BEST, 0, 372N,856W, 20, 1007, EX, 0, , 0, 0, 0, 0, 1010, 220, 80, 0, 0, , 0, , 42, 5, HARVEY, 72
AL, 09, 2017090206, , BEST, 0, 379N,849W, 20, 1009, EX, 0, , 0, 0, 0, 0, 1010, 220, 80, 0, 0, , 0, , 38, 5, HARVEY, 73
AL, 09, 2017090212, , BEST, 0, 382N,847W, 15, 1013, EX, 0, , 0, 0, 0, 0, 1010, 220, 80, 0, 0, , 0, , 28, 2, HARVEY, 74
AL, 09, 2017090212, , BEST, 0, 382N,847W, 15, 1013, EX, 0, , 0, 0, 0, 0, 1014, 220, 80, 0, 0, , 0, , 28, 2, HARVEY, 74
Loading

0 comments on commit a3e9c47

Please sign in to comment.