Merge pull request #74 from oceanmodeling/feature/rmw_option

Feature/rmw option

Author: SorooshMani-NOAA

pyproject.toml

@@ -54,7 +54,7 @@ dependencies = [
     "pytz",
     "pyyaml",
     "shapely>=2",
-    "stormevents>=2.3.2",  # rmax forecast
+    "stormevents>=2.3.4",  # rmax option
     "rasterio",
     "requests",
     "rtree",

stormworkflow/main.py

@@ -18,7 +18,7 @@
 
 _logger = logging.getLogger(__file__)
 
-CUR_INPUT_VER = Version('0.0.2')
+CUR_INPUT_VER = Version('0.0.3')
 
 
 def _handle_input_v0_0_1_to_v0_0_2(inout_conf):
@@ -43,6 +43,23 @@ def _handle_input_v0_0_1_to_v0_0_2(inout_conf):
     return Version('0.0.2')
 
 
+def _handle_input_v0_0_2_to_v0_0_3(inout_conf):
+
+    ver = Version(inout_conf['input_version'])
+
+    # Only update config if specified version matches the assumed one
+    if ver != Version('0.0.2'):
+        return ver
+
+
+    _logger.info(
+        "Adding RMW fill method default to persistent"
+    )
+    inout_conf['rmw_fill_method'] = 'persistent'
+
+    return Version('0.0.3')
+
+
 def handle_input_version(inout_conf):
 
     if 'input_version' not in inout_conf:
@@ -60,16 +77,18 @@ def handle_input_version(inout_conf):
             f"Input version not supported! Max version supported is {CUR_INPUT_VER}"
         )
 
-    ver = _handle_input_v0_0_1_to_v0_0_2(inout_conf)
+    for fn in [
+            _handle_input_v0_0_1_to_v0_0_2,
+            _handle_input_v0_0_2_to_v0_0_3,
+    ]:
+        ver = fn(inout_conf)
+        inout_conf['input_version'] = str(ver)
 
     if ver != CUR_INPUT_VER:
         raise ValueError(
             f"Could NOT update input to the latest version! Updated to {ver}"
         )
 
-    inout_conf['input_version'] = str(ver)
-
-
 def main():
 
     parser = ArgumentParser()

stormworkflow/prep/hurricane_data.py

@@ -24,6 +24,7 @@
 from shapely.geometry import box, base
 from stormevents import StormEvent
 from stormevents.nhc import VortexTrack
+from stormevents.nhc.const import RMWFillMethod
 from stormevents.nhc.track import (
     combine_tracks,
     correct_ofcl_based_on_carq_n_hollandb,
@@ -144,6 +145,7 @@ def main(args):
     hr_before_landfall = args.hours_before_landfall
     lead_times = args.lead_times
     track_dir = args.preprocessed_tracks_dir
+    rmw_fill = RMWFillMethod[args.rmw_fill.lower()]
 
     if hr_before_landfall < 0:
         hr_before_landfall = 48
@@ -183,7 +185,8 @@ def main(args):
 
         advisory = 'OFCL'
         if not local_track_file.is_file():
-            # Find and pick a single advisory based on priority
+            # Find and pick a single advisory based on priority, the
+            # track is only used to get the available advisories
             temp_track = event.track(file_deck='a')
             adv_avail = temp_track.unfiltered_data.advisory.unique()
             adv_order = ['OFCL', 'HWRF', 'HMON', 'CARQ']
@@ -193,41 +196,65 @@ def main(args):
                     advisory = adv
                     break
 
-            # TODO: THIS IS NO LONGER RELEVANT IF WE FAKE RMWP AS OFCL!
             if advisory == "OFCL" and "CARQ" not in adv_avail:
                 raise ValueError(
                     "OFCL advisory needs CARQ for fixing missing variables!"
                 )
 
-            track = VortexTrack(nhc_code, file_deck='a', advisories=[advisory])
+            track = VortexTrack(
+                nhc_code,
+                file_deck='a',
+                advisories=[advisory],
+                rmw_fill=rmw_fill,
+            )
 
         else:  # read from preprocessed file
+
             advisory = 'OFCL'
 
             # If a file exists, use the local file
             track_raw = pd.read_csv(local_track_file, header=None, dtype=str)
-            assert len(track_raw[4].unique()) == 1
+            # The provided tracks should have a single advisory type,
+            # e.g. in NHC adjusted track files the value is RMWP
+            if len(track_raw[4].unique()) != 1:
+                raise RuntimeError(
+                    "Only single advisory-name track files are supported!"
+                )
+            # Treat the existing advisory as if it's OFCL so that 
+            # stormevents supports reading it
             track_raw[4] = advisory
 
             with tempfile.NamedTemporaryFile() as tmp:
                 track_raw.to_csv(tmp.name, header=False, index=False)
 
+                # Track read from file is NOT corrected because it
+                # does NOT have CARQ advisory
                 unfixed_track = VortexTrack(
                     tmp.name, file_deck='a', advisories=[advisory]
                 )
+                # Since we're only getting CARQ, there's no need to 
+                # pass rmw fill method
                 carq_track = event.track(file_deck='a', advisories=['CARQ'])
                 unfix_dict = {
                     **separate_tracks(unfixed_track.data),
                     **separate_tracks(carq_track.data),
                 }
 
-                fix_dict = correct_ofcl_based_on_carq_n_hollandb(unfix_dict)
+                # Fix the file track with the fetched CARQ; if RMW
+                # is already filled, it fills out other missing values
+                fix_dict = correct_ofcl_based_on_carq_n_hollandb(
+                    unfix_dict, rmw_fill=rmw_fill
+                )
                 fix_track = combine_tracks(fix_dict)
 
+                # Create a new VortexTrack object from the datasets.
+                # Since the values are already filled in, there's
+                # no need to fill the rmw!
                 track = VortexTrack(
                     fix_track[fix_track.advisory == advisory],
                     file_deck='a',
-                    advisories=[advisory]
+                    advisories=[advisory],
+                    rmw_fill=RMWFillMethod.none,
                 )
 
 
@@ -418,6 +445,12 @@ def cli():
         help="Existing adjusted track directory",
     )
 
+    parser.add_argument(
+        '--rmw-fill',
+        type=str,
+        help="Method to use to fill missing RMW data for OFCL track",
+        )
+
     args = parser.parse_args()
 
     main(args)

stormworkflow/prep/setup_ensemble.py

@@ -176,6 +176,9 @@ def main(args):
     # get has unique forecast time for only the segment we want to
     # perturb, the preceeding entries are 0-hour forecasts from
     # previous forecast_times
+    # 
+    # Here we're working with NA-filled track files, so there's
+    # no need for rmw fill argument
     track_to_perturb = VortexTrack.from_file(
         track_path,
         start_date=perturb_start,

stormworkflow/refs/input.yaml

@@ -1,5 +1,5 @@
 ---
-input_version: 0.0.2
+input_version: 0.0.3
 
 storm: "florence"
 year: 2018
@@ -12,14 +12,16 @@ use_wwm: 0
 pahm_model: "gahm"
 num_perturb: 2
 sample_rule: "korobov"
+perturb_vars:
+  - "cross_track"
+  - "along_track"
+#  - "radius_of_maximum_winds"
+  - "radius_of_maximum_winds_persistent"
+  - "max_sustained_wind_speed"
+rmw_fill_method: "persistent"
+
 spinup_exec: "pschism_PAHM_TVD-VL"
 hotstart_exec: "pschism_PAHM_TVD-VL"
-perturb_vars:
-  - 'cross_track'
-  - 'along_track'
-#  - 'radius_of_maximum_winds'
-  - 'radius_of_maximum_winds_persistent'
-  - 'max_sustained_wind_speed'
 
 hpc_solver_nnodes: 3
 hpc_solver_ntasks: 108

stormworkflow/scripts/workflow.sh

@@ -77,6 +77,7 @@ hurricane_data \
     --hours-before-landfall "$hr_prelandfall" \
     --lead-times "$L_LEADTIMES_DATASET" \
     --preprocessed-tracks-dir "$L_TRACK_DIR" \
+    --rmw-fill "$rmw_fill_method" \
     $storm $year 2>&1 | tee "${run_dir}/output/head_hurricane_data.out"
 
 

stormworkflow/slurm/post.sbatch

@@ -2,6 +2,7 @@
 #SBATCH --parsable
 #SBATCH --time=05:00:00
 #SBATCH --nodes=1
+#SBATCH --exclusive
 
 set -ex
 

tests/conftest.py

@@ -0,0 +1,39 @@
+from importlib.resources import files
+
+import pytest
+import yaml
+from yaml import Loader
+
+
+refs = files('stormworkflow.refs')
+test_refs = files('tests.data.refs')
+input_v0_0_1 = test_refs.joinpath('input_v0.0.1.yaml')
+input_v0_0_2 = test_refs.joinpath('input_v0.0.2.yaml')
+input_v0_0_3 = test_refs.joinpath('input_v0.0.3.yaml')
+input_latest = refs.joinpath('input.yaml')
+
+
+def read_conf(infile):
+    with open(infile, 'r') as yfile:
+        conf = yaml.load(yfile, Loader=Loader)
+    return conf
+
+
+@pytest.fixture
+def conf_v0_0_1():
+    return read_conf(input_v0_0_1)
+
+
+@pytest.fixture
+def conf_v0_0_2():
+    return read_conf(input_v0_0_2)
+
+
+@pytest.fixture
+def conf_v0_0_3():
+    return read_conf(input_v0_0_3)
+
+
+@pytest.fixture
+def conf_latest():
+    return read_conf(input_latest)

tests/data/refs/input_v0.0.3.yaml

@@ -0,0 +1,48 @@
+---
+input_version: 0.0.3
+
+storm: "florence"
+year: 2018
+suffix: ""
+subset_mesh: 1
+hr_prelandfall: -1
+past_forecast: 1
+hydrology: 0
+use_wwm: 0
+pahm_model: "gahm"
+num_perturb: 2
+sample_rule: "korobov"
+perturb_vars:
+  - "cross_track"
+  - "along_track"
+#  - "radius_of_maximum_winds"
+  - "radius_of_maximum_winds_persistent"
+  - "max_sustained_wind_speed"
+rmw_fill_method: "persistent"
+
+spinup_exec: "pschism_PAHM_TVD-VL"
+hotstart_exec: "pschism_PAHM_TVD-VL"
+
+hpc_solver_nnodes: 3
+hpc_solver_ntasks: 108
+hpc_account: ""
+hpc_partition: ""
+
+RUN_OUT: ""
+L_NWM_DATASET: ""
+L_TPXO_DATASET: ""
+L_LEADTIMES_DATASET: ""
+L_TRACK_DIR: ""
+L_DEM_HI: ""
+L_DEM_LO: ""
+L_MESH_HI: ""
+L_MESH_LO: ""
+L_SHP_DIR: ""
+
+TMPDIR: "/tmp"
+PATH_APPEND: ""
+
+L_SOLVE_MODULES:
+  - "intel/2022.1.2"
+  - "impi/2022.1.2"
+  - "netcdf"

tests/test_input_version.py

@@ -1,40 +1,10 @@
 from copy import deepcopy
-from importlib.resources import files
 
 import pytest
-import yaml
-from packaging.version import Version
-from yaml import Loader, Dumper
 
 from stormworkflow.main import handle_input_version, CUR_INPUT_VER
 
 
-refs = files('tests.data.refs')
-input_v0_0_1 = refs.joinpath('input_v0.0.1.yaml')
-input_v0_0_2 = refs.joinpath('input_v0.0.2.yaml')
-
-
-def read_conf(infile):
-    with open(infile, 'r') as yfile:
-        conf = yaml.load(yfile, Loader=Loader)
-    return conf
-
-
-@pytest.fixture
-def conf_v0_0_1():
-    return read_conf(input_v0_0_1)
-
-
-@pytest.fixture
-def conf_v0_0_2():
-    return read_conf(input_v0_0_2)
-
-
-@pytest.fixture
-def conf_latest(conf_v0_0_2):
-    return conf_v0_0_2
-
-
 def test_no_version_specified(conf_latest):
     conf_latest.pop('input_version')
     with pytest.warns(UserWarning):
@@ -62,6 +32,16 @@ def test_invalid_version_specified(conf_latest):
     assert "invalid version" in str(e.value).lower()
 
 
-def test_v0_0_1_to_v0_0_2(conf_v0_0_1, conf_v0_0_2):
+def test_v0_0_1_to_latest(conf_v0_0_1, conf_latest):
     handle_input_version(conf_v0_0_1)
-    assert conf_v0_0_2 == conf_v0_0_1
+    assert conf_latest == conf_v0_0_1
+
+
+def test_v0_0_2_to_latest(conf_v0_0_2, conf_latest):
+    handle_input_version(conf_v0_0_2)
+    assert conf_latest == conf_v0_0_2
+
+
+def test_v0_0_3_to_latest(conf_v0_0_3, conf_latest):
+    handle_input_version(conf_v0_0_3)
+    assert conf_latest == conf_v0_0_3