diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..bcf9649 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,36 @@ +name: tests + +on: + push: + branches: + - main + paths: + - '**.py' + - '.github/workflows/tests.yml' + - 'pyproject.toml' + pull_request: + branches: + - main + +jobs: + test: + name: test + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ ubuntu-latest ] + python-version: [ '3.9', '3.10', '3.11' ] + steps: + - name: clone repository + uses: actions/checkout@v4 + - name: conda virtual environment + uses: mamba-org/setup-micromamba@v1 + with: + init-shell: bash + environment-file: environment.yml + - name: install the package + run: pip install ".[dev]" + shell: micromamba-shell {0} + - name: run tests + run: pytest + shell: micromamba-shell {0} diff --git a/README.md b/README.md index 8ec0ca3..dcb8066 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,16 @@ TO BE COMPLETED... ## References -- Daneshvar, F., et al. Tech Report (TODO) +- Daneshvar, F., Mani, S., Pringle, W. J., Moghimi, S., Myers, E., (2024). +_Probabilistic Prediction of Tropical Cyclone (TC) Driven Flood Depth_ +_and Extent with a User-Friendly Python Module._ +NOAA Technical Memorandum NOS CS 59. https://repository.library.noaa.gov/view/noaa/66123 + +- Daneshvar, F., Mani, S., Pringle, W. J., Moghimi, S., Myers, E., ( 2023). +_Evaluating the effect of inland hydrology and hurricane model on +the probabilistic prediction of tropical cyclone (TC) driven coastal flooding._ +American Geophysical Union (AGU) Fall Meeting, December 2023, San Francisco, CA. + - Pringle, W. J., Mani, S., Sargsyan, K., Moghimi, S., Zhang, Y. J., Khazaei, B., Myers, E. (January 2023). _Surrogate-Assisted Bayesian Uncertainty Quantification for diff --git a/environment.yml b/environment.yml index 28cf939..ae94725 100644 --- a/environment.yml +++ b/environment.yml @@ -14,7 +14,7 @@ dependencies: - fiona - gdal - geoalchemy2 - - geopandas>=0.13 + - geopandas - geos - hdf5 - importlib_metadata<8 # Fix issue with esmpy Author import @@ -30,7 +30,7 @@ dependencies: - pyarrow - pygeos - pyproj - - python<3.11 + - python<3.12 - pytz - shapely>=2 - rasterio @@ -42,4 +42,4 @@ dependencies: - tqdm - udunits2 - utm - - xarray==2023.7.0 + - xarray diff --git a/pyproject.toml b/pyproject.toml index 8a8a1e7..446219d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,7 +21,7 @@ description = "A set of scripts to generate probabilistic storm surge results!" license = {file = "LICENSE"} -requires-python = ">= 3.8, < 3.11" +requires-python = ">= 3.9, < 3.12" dependencies = [ "cartopy", @@ -30,22 +30,22 @@ dependencies = [ "cfgrib", "cfunits", "chaospy>=4.2.7", - "coupledmodeldriver>=1.6.6", + "coupledmodeldriver>=1.7.1.post1", "colored-traceback", "cmocean", "dask", "dask-jobqueue", - "ensembleperturbation>=1.1.2", + "ensembleperturbation>=1.3.4", # perturb feature "fiona", "geoalchemy2", - "geopandas==0.14", # to address AttributeError. Should be fixed later in EnsemblePerturbation -# "geopandas", + "geopandas", "matplotlib", "mpi4py", "netCDF4", "numpy", "numba", "ocsmesh==1.5.3", + "packaging", "pandas", "pyarrow", "pygeos", @@ -54,7 +54,7 @@ dependencies = [ "pytz", "pyyaml", "shapely>=2", - "stormevents==2.2.4", + "stormevents>=2.3.4", # rmax option "rasterio", "requests", "rtree", @@ -63,7 +63,12 @@ dependencies = [ "typing-extensions", "tqdm", "utm", - "xarray==2023.7.0", + "xarray", +] + +[project.optional-dependencies] +dev = [ + "pytest" ] [tool.setuptools_scm] @@ -94,4 +99,4 @@ download_data = "stormworkflow.prep.download_data:cli" setup_ensemble = "stormworkflow.prep.setup_ensemble:cli" combine_ensemble = "stormworkflow.post.combine_ensemble:cli" analyze_ensemble = "stormworkflow.post.analyze_ensemble:cli" -storm_roc_curve = "stormworkflow.post.ROC_single_run:cli" +storm_roc_curve = "stormworkflow.post.storm_roc_curve:cli" diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 5cc4a85..0000000 --- a/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -chaospy>=4.2.7 -stormevents==2.2.3 -pyschism>=0.1.15 -coupledmodeldriver>=1.6.6 -ensembleperturbation>=1.1.2 diff --git a/stormworkflow/main.py b/stormworkflow/main.py index e649f6c..4b4c029 100644 --- a/stormworkflow/main.py +++ b/stormworkflow/main.py @@ -2,20 +2,113 @@ import logging import os import shlex +import warnings from importlib.resources import files from argparse import ArgumentParser from pathlib import Path -import stormworkflow import yaml +from packaging.version import Version try: from yaml import CLoader as Loader, CDumper as Dumper except ImportError: from yaml import Loader, Dumper +import stormworkflow _logger = logging.getLogger(__file__) +CUR_INPUT_VER = Version('0.0.5') +VER_UPDATE_FUNCS = [] + + +def _input_version(prev, curr): + def decorator(handler): + def wrapper(inout_conf): + ver = Version(inout_conf['input_version']) + + # Only update config if specified version matches the + # assumed one + if ver != Version(prev): + return ver + + # TODO: Need return values? + handler(inout_conf) + + return Version(curr) + global VER_UPDATE_FUNCS + VER_UPDATE_FUNCS.append(wrapper) + return wrapper + return decorator + + +@_input_version('0.0.1', '0.0.2') +def _handle_input_v0_0_1_to_v0_0_2(inout_conf): + + _logger.info( + "Adding perturbation variables for persistent RMW perturbation" + ) + inout_conf['perturb_vars'] = [ + 'cross_track', + 'along_track', + 'radius_of_maximum_winds_persistent', + 'max_sustained_wind_speed', + ] + + +@_input_version('0.0.2', '0.0.3') +def _handle_input_v0_0_2_to_v0_0_3(inout_conf): + + _logger.info( + "Adding RMW fill method default to persistent" + ) + inout_conf['rmw_fill_method'] = 'persistent' + + +@_input_version('0.0.3', '0.0.4') +def _handle_input_v0_0_3_to_v0_0_4(inout_conf): + + _logger.info( + "Path to observations" + ) + inout_conf['NHC_OBS'] = '' + + +@_input_version('0.0.4', '0.0.5') +def _handle_input_v0_0_4_to_v0_0_5(inout_conf): + + _logger.info("Adding perturbation features") + inout_conf['perturb_features'] = [ + 'isotach_adjustment', + ] + + +def handle_input_version(inout_conf): + + if 'input_version' not in inout_conf: + ver = CUR_INPUT_VER + warnings.warn( + f"`input_version` is NOT specified in `input.yaml`; assuming {ver}" + ) + inout_conf['input_version'] = str(ver) + return + + ver = Version(inout_conf['input_version']) + + if ver > CUR_INPUT_VER: + raise ValueError( + f"Input version not supported! Max version supported is {CUR_INPUT_VER}" + ) + + for fn in VER_UPDATE_FUNCS: + 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}" + ) + def main(): parser = ArgumentParser() @@ -28,12 +121,17 @@ def main(): infile = args.configuration if infile is None: - _logger.warn('No input configuration provided, using reference file!') + warnings.warn( + 'No input configuration provided, using reference file!' + ) infile = refs.joinpath('input.yaml') with open(infile, 'r') as yfile: conf = yaml.load(yfile, Loader=Loader) + handle_input_version(conf) + # TODO: Write out the updated config as a yaml file + wf = scripts.joinpath('workflow.sh') run_env = os.environ.copy() diff --git a/stormworkflow/post/analyze_ensemble.py b/stormworkflow/post/analyze_ensemble.py index 44a78a3..f12ea8b 100644 --- a/stormworkflow/post/analyze_ensemble.py +++ b/stormworkflow/post/analyze_ensemble.py @@ -408,19 +408,19 @@ def cli(): parser.add_argument('-t', '--tracks-dir', type=Path) parser.add_argument('-s', '--sequential', action='store_true') - main(parser.parse_args()) - - -if __name__ == '__main__': cluster = SLURMCluster(cores=16, processes=1, memory="500GB", - account="compute", + account="nos-surge", #PW:"compute" ; Hercules:"nos-surge" or "nosofs" walltime="08:00:00", - header_skip=['--mem'], - interface="eth0") - cluster.scale(6) - client = Client(cluster) + header_skip=['--mem']) #, +# interface="eth0") # only for PW + cluster.scale(6) + client = Client(cluster) print(client) + main(parser.parse_args()) + + +if __name__ == '__main__': cli() diff --git a/stormworkflow/post/ROC_single_run.py b/stormworkflow/post/storm_roc_curve.py similarity index 87% rename from stormworkflow/post/ROC_single_run.py rename to stormworkflow/post/storm_roc_curve.py index 26dbd2c..fe7bde5 100644 --- a/stormworkflow/post/ROC_single_run.py +++ b/stormworkflow/post/storm_roc_curve.py @@ -10,6 +10,8 @@ from pathlib import Path from cartopy.feature import NaturalEarthFeature +import geodatasets + os.environ['USE_PYGEOS'] = '0' import geopandas as gpd @@ -117,12 +119,17 @@ def calculate_POD_FAR(hit, miss, false_alarm, correct_neg): def main(args): - storm_name = args.storm_name.capitalize() - storm_year = args.storm_year + storm = args.storm.capitalize() + year = args.year leadtime = args.leadtime - prob_nc_path = Path(args.prob_nc_path) obs_df_path = Path(args.obs_df_path) - save_dir = args.save_dir + ensemble_dir = Path(args.ensemble_dir) + + output_directory = ensemble_dir / 'analyze/linear_k1_p1_n0.025' + prob_nc_path = output_directory / 'probabilities.nc' + + if leadtime == -1: + leadtime = 48 # *.nc file coordinates thresholds_ft = [3, 6, 9] # in ft @@ -150,7 +157,7 @@ def main(args): # Load obs file, extract storm obs points and coordinates df_obs = pd.read_csv(obs_df_path) - Event_name = f'{storm_name}_{storm_year}' + Event_name = f'{storm}_{year}' df_obs_storm = df_obs[df_obs.Event == Event_name] obs_coordinates = stack_station_coordinates( df_obs_storm.Longitude.values, df_obs_storm.Latitude.values @@ -159,10 +166,12 @@ def main(args): # Load probabilities.nc file ds_prob = xr.open_dataset(prob_nc_path) - gdf_countries = gpd.GeoSeries( - NaturalEarthFeature(category='physical', scale='10m', name='land',).geometries(), - crs=4326, - ) + gdf_countries = gpd.read_file(geodatasets.get_path('naturalearth land')) + + # gdf_countries = gpd.GeoSeries( + # NaturalEarthFeature(category='physical', scale='10m', name='land',).geometries(), + # crs=4326, + # ) # Loop through thresholds and sources and find corresponding values from probabilities.nc threshold_count = -1 @@ -186,10 +195,10 @@ def main(args): df_obs_storm, f'{source}_prob', gdf_countries, - f'Probability of {source} exceeding {thresholds_ft[threshold_count]} ft \n {storm_name}, {storm_year}, {leadtime}-hr leadtime', + f'Probability of {source} exceeding {thresholds_ft[threshold_count]} ft \n {storm}, {year}, {leadtime}-hr leadtime', os.path.join( - save_dir, - f'prob_{source}_above_{thresholds_ft[threshold_count]}ft_{storm_name}_{storm_year}_{leadtime}-hr.png', + output_directory, + f'prob_{source}_above_{thresholds_ft[threshold_count]}ft_{storm}_{year}_{leadtime}-hr.png', ), ) @@ -212,7 +221,7 @@ def main(args): ds_ROC = xr.Dataset( coords=dict( threshold=thresholds_ft, - storm=[storm_name], + storm=[storm], leadtime=[leadtime], source=sources, prob=probabilities, @@ -232,9 +241,7 @@ def main(args): FAR=(['threshold', 'storm', 'leadtime', 'source', 'prob'], FAR_arr), ), ) - ds_ROC.to_netcdf( - os.path.join(save_dir, f'{storm_name}_{storm_year}_{leadtime}hr_leadtime_POD_FAR.nc') - ) + ds_ROC.to_netcdf(os.path.join(output_directory, f'{storm}_{year}_{leadtime}hr_POD_FAR.nc')) # plot ROC curves marker_list = ['s', 'x'] @@ -262,12 +269,11 @@ def main(args): plt.xlabel('False Alarm Rate') plt.ylabel('Probability of Detection') - plt.title( - f'{storm_name}_{storm_year}, {leadtime}-hr leadtime, {threshold} ft threshold' - ) + plt.title(f'{storm}_{year}, {leadtime}-hr leadtime, {threshold} ft threshold') plt.savefig( os.path.join( - save_dir, f'ROC_{storm_name}_{leadtime}hr_leadtime_{threshold}_ft.png' + output_directory, + f'ROC_{storm}_{year}_{leadtime}hr_leadtime_{threshold}_ft.png', ) ) plt.close() @@ -276,20 +282,11 @@ def main(args): def cli(): parser = argparse.ArgumentParser() - parser.add_argument('--storm_name', help='name of the storm', type=str) - - parser.add_argument('--storm_year', help='year of the storm', type=int) - + parser.add_argument('--storm', help='name of the storm', type=str) + parser.add_argument('--year', help='year of the storm', type=int) parser.add_argument('--leadtime', help='OFCL track leadtime hr', type=int) - - parser.add_argument('--prob_nc_path', help='path to probabilities.nc', type=str) - - parser.add_argument('--obs_df_path', help='Path to observations dataframe', type=str) - - # optional - parser.add_argument( - '--save_dir', help='directory for saving analysis', default=os.getcwd(), type=str - ) + parser.add_argument('--obs_df_path', help='path to NHC obs data', type=str) + parser.add_argument('--ensemble-dir', help='path to ensemble.dir', type=str) main(parser.parse_args()) diff --git a/stormworkflow/prep/hurricane_data.py b/stormworkflow/prep/hurricane_data.py index 5e3382a..faf729e 100644 --- a/stormworkflow/prep/hurricane_data.py +++ b/stormworkflow/prep/hurricane_data.py @@ -17,12 +17,14 @@ import pandas as pd import geopandas as gpd +import geodatasets from searvey.coops import COOPS_TidalDatum from searvey.coops import COOPS_TimeZone from searvey.coops import COOPS_Units 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, @@ -37,7 +39,7 @@ format='%(asctime)s,%(msecs)d %(levelname)-8s [%(filename)s:%(lineno)d] %(message)s', datefmt='%Y-%m-%d:%H:%M:%S') - +NE_LOW_ADMIN = 'https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip' def trackstart_from_file( leadtime_file: Optional[pathlib.Path], @@ -143,12 +145,15 @@ def main(args): hr_before_landfall = args.hours_before_landfall lead_times = args.lead_times track_dir = args.preprocessed_tracks_dir - countries_shpfile = args.countries_polygon + rmw_fill = RMWFillMethod[args.rmw_fill.lower()] if hr_before_landfall < 0: hr_before_landfall = 48 - ne_low = gpd.read_file(countries_shpfile) + # Caching for next steps! + geodatasets.get_path('naturalearth land') + + ne_low = gpd.read_file(NE_LOW_ADMIN) shp_US = ne_low[ne_low.NAME_EN.isin(['United States of America', 'Puerto Rico'])].unary_union logger.info("Fetching hurricane info...") @@ -180,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'] @@ -190,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, ) @@ -416,10 +446,10 @@ def cli(): ) parser.add_argument( - "--countries-polygon", - type=pathlib.Path, - help="Shapefile containing country polygons", - ) + '--rmw-fill', + type=str, + help="Method to use to fill missing RMW data for OFCL track", + ) args = parser.parse_args() diff --git a/stormworkflow/prep/setup_ensemble.py b/stormworkflow/prep/setup_ensemble.py index 71d2218..4e78cb7 100644 --- a/stormworkflow/prep/setup_ensemble.py +++ b/stormworkflow/prep/setup_ensemble.py @@ -22,7 +22,7 @@ from coupledmodeldriver.generate import SCHISMRunConfiguration from coupledmodeldriver.generate.schism.script import SchismEnsembleGenerationJob from coupledmodeldriver.generate import generate_schism_configuration -from ensembleperturbation.perturbation.atcf import perturb_tracks +from ensembleperturbation.perturbation.atcf import perturb_tracks, PerturberFeatures from pylib_essentials.schism_file import ( read_schism_hgrid_cached, schism_bpfile, @@ -107,6 +107,24 @@ def _fix_nwm_issues(ensemble_dir, hires_shapefile): _relocate_source_sink(pth, hires_shapefile) +def _fix_hotstart_issue(ensemble_dir): + hotstart_dirs = ensemble_dir.glob('runs/*') + for pth in hotstart_dirs: + nm_list = f90nml.read(pth / 'param.nml') + nm_list['opt']['dramp'] = 0.0 + nm_list['opt']['drampbc'] = 0.0 + nm_list['opt']['dramp_ss'] = 0.0 + nm_list['opt']['drampwind'] = 0.0 + nm_list.write(pth / 'param.nml', force=True) + +def _fix_veg_parameter_issue(ensemble_dir): + # See https://github.com/schism-dev/pyschism/issues/126 + param_nmls = ensemble_dir.glob('**/param.nml') + for pth in param_nmls: + nm_list = f90nml.read(pth) + nm_list['core']['nbins_veg_vert'] = 2 + nm_list.write(pth, force=True) + def main(args): track_path = args.track_file @@ -119,6 +137,9 @@ def main(args): use_wwm = args.use_wwm with_hydrology = args.with_hydrology pahm_model = args.pahm_model + setup_features = PerturberFeatures.NONE + for feat in args.perturb_features: + setup_features |= PerturberFeatures[feat.upper()] workdir = out_dir mesh_file = mesh_dir / 'mesh_w_bdry.grd' @@ -165,6 +186,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, @@ -178,12 +202,7 @@ def main(args): perturbations=args.num_perturbations, directory=workdir / 'track_files', storm=workdir / 'track_to_perturb.dat', - variables=[ - 'cross_track', - 'along_track', - 'radius_of_maximum_winds', # TODO: add option for persistent - 'max_sustained_wind_speed', - ], + variables=args.variables, sample_from_distribution=args.sample_from_distribution, sample_rule=args.sample_rule, quadrature=args.quadrature, @@ -192,6 +211,7 @@ def main(args): overwrite=True, file_deck=file_deck, advisories=[advisory], + features=setup_features, ) if perturb_start != model_start_time: @@ -267,6 +287,8 @@ def main(args): } ) + _fix_hotstart_issue(workdir) + _fix_veg_parameter_issue(workdir) # For newer SCHISM version if with_hydrology: _fix_nwm_issues(workdir, hires_reg) if use_wwm: @@ -327,6 +349,8 @@ def parse_arguments(): argument_parser.add_argument('--use-wwm', action='store_true') argument_parser.add_argument('--with-hydrology', action='store_true') argument_parser.add_argument('--pahm-model', choices=['gahm', 'symmetric'], default='gahm') + argument_parser.add_argument('--perturb-features', nargs='+', type=str, default=[PerturberFeatures.ISOTACH_ADJUSTMENT.name]) + argument_parser.add_argument('--variables', nargs='+', type=str) argument_parser.add_argument('name', help='name of the storm', type=str) diff --git a/stormworkflow/refs/input.yaml b/stormworkflow/refs/input.yaml index 1e8141e..3c629c1 100644 --- a/stormworkflow/refs/input.yaml +++ b/stormworkflow/refs/input.yaml @@ -1,5 +1,5 @@ --- -input_version: 0.0.1 +input_version: 0.0.5 storm: "florence" year: 2018 @@ -12,6 +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" +perturb_features: + - "isotach_adjustment" + spinup_exec: "pschism_PAHM_TVD-VL" hotstart_exec: "pschism_PAHM_TVD-VL" @@ -30,6 +40,7 @@ L_DEM_LO: "" L_MESH_HI: "" L_MESH_LO: "" L_SHP_DIR: "" +NHC_OBS: "" TMPDIR: "/tmp" PATH_APPEND: "" diff --git a/stormworkflow/scripts/workflow.sh b/stormworkflow/scripts/workflow.sh index e9fab6c..2e22618 100755 --- a/stormworkflow/scripts/workflow.sh +++ b/stormworkflow/scripts/workflow.sh @@ -48,12 +48,15 @@ function init { done logfile=$run_dir/versions.info + version $logfile stormworkflow version $logfile stormevents version $logfile ensembleperturbation + version $logfile coupledmodeldriver + version $logfile pyschism version $logfile ocsmesh echo "SCHISM: see solver.version each outputs dir" >> $logfile - cp $input_file $run_dir/input.yaml + cp $input_file $run_dir/input_asis.yaml echo $run_dir } @@ -74,7 +77,7 @@ hurricane_data \ --hours-before-landfall "$hr_prelandfall" \ --lead-times "$L_LEADTIMES_DATASET" \ --preprocessed-tracks-dir "$L_TRACK_DIR" \ - --countries-polygon "$L_SHP_DIR/ne_110m_cultural/ne_110m_admin_0_countries.shp" \ + --rmw-fill "$rmw_fill_method" \ $storm $year 2>&1 | tee "${run_dir}/output/head_hurricane_data.out" @@ -132,8 +135,10 @@ PREP_KWDS+=" --sample-rule $sample_rule" PREP_KWDS+=" --date-range-file $run_dir/setup/dates.csv" PREP_KWDS+=" --nwm-file $L_NWM_DATASET" PREP_KWDS+=" --tpxo-dir $L_TPXO_DATASET" +PREP_KWDS+=" --variables $perturb_vars" if [ $use_wwm == 1 ]; then PREP_KWDS+=" --use-wwm"; fi if [ $hydrology == 1 ]; then PREP_KWDS+=" --with-hydrology"; fi +PREP_KWDS+=" --perturb-features $perturb_features" PREP_KWDS+=" --pahm-model $pahm_model" export PREP_KWDS # NOTE: We need to wait because run jobs depend on perturbation dirs! @@ -142,7 +147,7 @@ setup_id=$(sbatch \ --wait \ --job-name=prep_$tag \ --parsable \ - --export=ALL,PREP_KWDS,STORM=$storm,YEAR=$year,IMG="$L_IMG_DIR/prep.sif" \ + --export=ALL,PREP_KWDS,STORM=$storm,YEAR=$year \ $run_dir/slurm/prep.sbatch \ ) @@ -150,7 +155,6 @@ setup_id=$(sbatch \ echo "Launching runs" SCHISM_SHARED_ENV="" SCHISM_SHARED_ENV+="ALL" -SCHISM_SHARED_ENV+=",IMG=$L_IMG_DIR/solve.sif" SCHISM_SHARED_ENV+=",MODULES=$L_SOLVE_MODULES" spinup_id=$(sbatch \ --nodes $hpc_solver_nnodes --ntasks $hpc_solver_ntasks \ @@ -181,5 +185,5 @@ sbatch \ --output "${run_dir}/output/slurm-%j.post.out" \ --job-name=post_$tag \ -d afterok${joblist} \ - --export=ALL,IMG="$L_IMG_DIR/prep.sif",ENSEMBLE_DIR="$run_dir/setup/ensemble.dir/" \ + --export=ALL,ENSEMBLE_DIR="$run_dir/setup/ensemble.dir/" \ $run_dir/slurm/post.sbatch diff --git a/stormworkflow/slurm/mesh.sbatch b/stormworkflow/slurm/mesh.sbatch index 757b9cd..503381a 100644 --- a/stormworkflow/slurm/mesh.sbatch +++ b/stormworkflow/slurm/mesh.sbatch @@ -1,7 +1,7 @@ #!/bin/bash #SBATCH --parsable #SBATCH --exclusive -#SBATCH --time=00:30:00 +#SBATCH --time=01:00:00 #SBATCH --nodes=1 set -ex diff --git a/stormworkflow/slurm/post.sbatch b/stormworkflow/slurm/post.sbatch index 4ef59f8..33ef357 100644 --- a/stormworkflow/slurm/post.sbatch +++ b/stormworkflow/slurm/post.sbatch @@ -2,6 +2,7 @@ #SBATCH --parsable #SBATCH --time=05:00:00 #SBATCH --nodes=1 +#SBATCH --exclusive set -ex @@ -12,3 +13,10 @@ combine_ensemble \ analyze_ensemble \ --ensemble-dir $ENSEMBLE_DIR \ --tracks-dir $ENSEMBLE_DIR/track_files + +storm_roc_curve \ + --storm ${storm} \ + --year ${year} \ + --leadtime ${hr_prelandfall} \ + --obs_df_path ${NHC_OBS} \ + --ensemble-dir $ENSEMBLE_DIR diff --git a/stormworkflow/slurm/prep.sbatch b/stormworkflow/slurm/prep.sbatch index 8beb298..d5dd455 100644 --- a/stormworkflow/slurm/prep.sbatch +++ b/stormworkflow/slurm/prep.sbatch @@ -1,7 +1,7 @@ #!/bin/bash #SBATCH --parsable #SBATCH --exclusive -#SBATCH --time=00:30:00 +#SBATCH --time=01:00:00 #SBATCH --nodes=1 set -ex diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..e62a5af --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,48 @@ +from importlib.resources import files + +import pytest +import yaml +from yaml import Loader + + +refs = files('stormworkflow.refs') +test_refs = files('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_v0_0_4 = test_refs.joinpath('input_v0.0.4.yaml') +input_v0_0_5 = test_refs.joinpath('input_v0.0.5.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_v0_0_4(): + return read_conf(input_v0_0_4) + +@pytest.fixture +def conf_v0_0_5(): + return read_conf(input_v0_0_5) + +@pytest.fixture +def conf_latest(): + return read_conf(input_latest) diff --git a/tests/data/refs/input_v0.0.1.yaml b/tests/data/refs/input_v0.0.1.yaml new file mode 100644 index 0000000..1e8141e --- /dev/null +++ b/tests/data/refs/input_v0.0.1.yaml @@ -0,0 +1,40 @@ +--- +input_version: 0.0.1 + +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" +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" diff --git a/tests/data/refs/input_v0.0.2.yaml b/tests/data/refs/input_v0.0.2.yaml new file mode 100644 index 0000000..f438167 --- /dev/null +++ b/tests/data/refs/input_v0.0.2.yaml @@ -0,0 +1,46 @@ +--- +input_version: 0.0.2 + +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" +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 +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" diff --git a/tests/data/refs/input_v0.0.3.yaml b/tests/data/refs/input_v0.0.3.yaml new file mode 100644 index 0000000..02df434 --- /dev/null +++ b/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" diff --git a/tests/data/refs/input_v0.0.4.yaml b/tests/data/refs/input_v0.0.4.yaml new file mode 100644 index 0000000..30aefa2 --- /dev/null +++ b/tests/data/refs/input_v0.0.4.yaml @@ -0,0 +1,49 @@ +--- +input_version: 0.0.4 + +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: "" +NHC_OBS: "" + +TMPDIR: "/tmp" +PATH_APPEND: "" + +L_SOLVE_MODULES: + - "intel/2022.1.2" + - "impi/2022.1.2" + - "netcdf" diff --git a/tests/data/refs/input_v0.0.5.yaml b/tests/data/refs/input_v0.0.5.yaml new file mode 100644 index 0000000..3c629c1 --- /dev/null +++ b/tests/data/refs/input_v0.0.5.yaml @@ -0,0 +1,51 @@ +--- +input_version: 0.0.5 + +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" +perturb_features: + - "isotach_adjustment" + +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: "" +NHC_OBS: "" + +TMPDIR: "/tmp" +PATH_APPEND: "" + +L_SOLVE_MODULES: + - "intel/2022.1.2" + - "impi/2022.1.2" + - "netcdf" diff --git a/tests/test_input_version.py b/tests/test_input_version.py new file mode 100644 index 0000000..9d22429 --- /dev/null +++ b/tests/test_input_version.py @@ -0,0 +1,57 @@ +from copy import deepcopy + +import pytest + +from stormworkflow.main import handle_input_version, CUR_INPUT_VER + + +def test_no_version_specified(conf_latest): + conf_latest.pop('input_version') + with pytest.warns(UserWarning): + handle_input_version(conf_latest) + + assert conf_latest['input_version'] == str(CUR_INPUT_VER) + + +def test_invalid_version_specified(conf_latest): + + invalid_1 = deepcopy(conf_latest) + invalid_1['input_version'] = ( + f'{CUR_INPUT_VER.major}.{CUR_INPUT_VER.minor}.{CUR_INPUT_VER.micro + 1}' + ) + with pytest.raises(ValueError) as e: + handle_input_version(invalid_1) + + assert "max" in str(e.value).lower() + + + invalid_2 = deepcopy(conf_latest) + invalid_2['input_version'] = 'a.b.c' + with pytest.raises(ValueError) as e: + handle_input_version(invalid_2) + assert "invalid version" in str(e.value).lower() + + +def test_v0_0_1_to_latest(conf_v0_0_1, conf_latest): + handle_input_version(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 + + +def test_v0_0_4_to_latest(conf_v0_0_4, conf_latest): + handle_input_version(conf_v0_0_4) + assert conf_latest == conf_v0_0_4 + + +def test_v0_0_5_to_latest(conf_v0_0_5, conf_latest): + handle_input_version(conf_v0_0_5) + assert conf_latest == conf_v0_0_5