Refactor and extend NetCDF input support for topography and storm surge (met) forcing

examples/storm-surge/ike/test_ike.py

@@ -17,7 +17,7 @@
     ),
     pytest.param(
         {"num_cells": [29 * 4, 24 * 4], 
-         "amr_levels_max": 6, 
+         "amr_levels_max": 4, 
          "num_output_times": 16},
         id="fine",
         marks=pytest.mark.slow,

examples/storm-surge/isaac/Makefile

@@ -24,7 +24,7 @@ PLOTDIR = _plots               # Directory for plots
 # FFLAGS += -DNETCDF -I$(NETCDF4_DIR)/include -L$(NETCDF4_DIR)/lib
 # LFLAGS += -I$(NETCDF4_DIR)/include -L$(NETCDF4_DIR)/lib -lnetcdff
 FFLAGS += -DNETCDF $(NETCDF_FFLAGS)
-LFLAGS += $(NETCDF_LFLAGS)
+LFLAGS += $(FFLAGS) $(NETCDF_LFLAGS)
 
 # ---------------------------------
 # package sources for this program:

examples/storm-surge/isaac/setplot.py

@@ -13,6 +13,7 @@
 import clawpack.amrclaw.data as amrclaw
 import clawpack.geoclaw.data as geodata
 import clawpack.geoclaw.util as geoutil
+import clawpack.geoclaw.surge.storm as stormtools
 import clawpack.geoclaw.surge.plot as surgeplot
 
 try:
@@ -42,7 +43,8 @@ def setplot(plotdata=None):
     friction_data = geodata.FrictionData()
     friction_data.read(os.path.join(plotdata.outdir, 'friction.data'))
 
-    # Load storm track
+    # Load storm data
+    storm = stormtools.Storm(surge_data.storm_file, file_format="data")
     track = surgeplot.track_data(os.path.join(plotdata.outdir, 'fort.track'))
 
     # Set afteraxes function
@@ -173,17 +175,19 @@ def friction_after_axes(cd):
     # ========================================================================
     #  Figures for gauges
     # ========================================================================
+    # Map GeoClaw gauge number to NOAA gauge number and location/name
+    gauge_mapping = {1: [8761724, "Grand Isle, LA"],
+                        2: [8760922, 'Pilots Station East, SW Pass, LA']}
+
+    # Storm Timepoints
+    landfall_time = np.datetime64(datetime.datetime(2012, 8, 29, 0))
+    begin_date = datetime.datetime(2012, 8, 27)
+    end_date = datetime.datetime(2012, 8, 31)
+
     def plot_observed(current_data):
         """Fetch and plot gauge data for gauges used."""
 
-        # Map GeoClaw gauge number to NOAA gauge number and location/name
-        gauge_mapping = {1: [8761724, "Grand Isle, LA"],
-                         2: [8760922, 'Pilots Station East, SW Pass, LA']}
-
         station_id, station_name = gauge_mapping[current_data.gaugesoln.id]
-        landfall_time = np.datetime64(datetime.datetime(2012, 8, 29, 0))
-        begin_date = datetime.datetime(2012, 8, 27)
-        end_date = datetime.datetime(2012, 8, 31)
 
         # Fetch data if needed
         date_time, water_level, tide = geoutil.fetch_noaa_tide_data(station_id,
@@ -227,6 +231,49 @@ def plot_observed(current_data):
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
     plotitem.plot_var = surgeplot.gauge_dry_regions
     plotitem.kwargs = {"color":'lightcoral', "linewidth":5}
+
+
+    # Plot wind and pressure at each gauge
+    def storm_gauge_afteraxes(cd, storm):
+        # Map GeoClaw gauge number to NOAA gauge number and location/name
+        station_id, station_name = gauge_mapping[cd.gaugesoln.id]
+
+        # Convert to seconds relative to landfall
+        # t = (cd.gaugesoln.t - landfall_time) / np.timedelta64(1, 's')
+        t = cd.gaugesoln.t
+        t /= (24 * 60**2)
+
+        w_ax = plt.gca()
+        color ="tab:blue"
+        w_ax.plot(t, surgeplot.gauge_wind(cd), label="wind", color=color)
+        w_ax.set_ylabel("Wind (m/s)", color=color)
+        w_ax.set_ylim(wind_limits)
+        w_ax.tick_params(axis='y', labelcolor=color)
+
+        P_ax = w_ax.twinx()
+        color = "tab:red"
+        P_ax.set_ylabel("Pressure (mbar)", color=color)
+        P_ax.plot(t, surgeplot.gauge_pressure(cd) * 1e-2, label="pressure", color=color)
+        P_ax.set_ylim(pressure_limits)
+        P_ax.tick_params(axis='y', labelcolor=color)
+
+        w_ax.set_title(f"Storm Fileds at {station_name}")
+
+        plt.gcf().tight_layout()
+
+    plotfigure = plotdata.new_plotfigure(name='Gauge Wind', figno=301,
+                                         type='each_gauge')
+    plotfigure.show = True
+    plotfigure.clf_each_gauge = True
+    plotfigure.kwargs = {"figsize": [6.4 * 1.2,
+                                     4.8 * 1.0]}
+
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.time_scale = 1 / (24 * 60**2)
+    plotaxes.grid = True
+    plotaxes.xlimits = [-2, 1.5]
+    plotaxes.time_label = "Days relative to 2012-12-26 00:00 UTC"
+    plotaxes.afteraxes = lambda cd: storm_gauge_afteraxes(cd, storm)
 
     #  Gauge Location Plot
     def gauge_location_afteraxes(cd):

examples/storm-surge/isaac/setrun.py

@@ -322,6 +322,9 @@ def setrun(claw_pkg='geoclaw'):
     rundata.gaugedata.gauges.append([2, -89.95826667, 29.26460167,
                                      rundata.clawdata.t0,
                                      rundata.clawdata.tfinal])
+
+    # Record the wind and pressure fields at the gauges
+    rundata.gaugedata.aux_out_fields = [4, 5, 6]
 
     # ------------------------------------------------------------------
     # GeoClaw specific parameters:
@@ -409,34 +412,6 @@ def setgeo(rundata):
     data.storm_specification_type = 'holland80'
     # data.storm_specification_type = 'data'
     data.storm_file = (Path() / 'isaac.storm').resolve()
-
-    isaac = Storm()
-    atcf_path = (Path(__file__).parent / "bal092012.dat").resolve()
-    isaac.read(path=atcf_path, file_format="ATCF")
-    # Calculate landfall time - Need to specify as the file above does not
-    # include this info (~2345 UTC - 6:45 p.m. CDT - on August 28)
-    isaac.time_offset = np.datetime64("2012-08-29T00:00")
-
-    if data.storm_specification_type == 'holland80':
-        isaac.write(data.storm_file, file_format='geoclaw')
-
-    elif data.storm_specification_type == "data":
-        # ASCII
-        isaac.file_format = 'NWS12'
-        isaac.file_paths = []
-        isaac.file_paths.append((Path() / "isaac.PRE").resolve())
-        isaac.file_paths.append((Path() / "isaac.WIN").resolve())
-
-        # NetCDF
-        # isaac.file_format = "NWS13"
-        # isaac.file_paths.append((Path() / "isaac.nc").resolve())
-        # # For NetCDF we need to set this to the epoch or things do not work yet
-        # isaac.time_offset = np.datetime64("1970-01-01")
-
-        isaac.write(data.storm_file, file_format='data')
-
-    else:
-        raise ValueError("Invalid storm specification type.")
 
     # =======================
     #  Set Variable Friction
@@ -463,6 +438,49 @@ def setgeo(rundata):
     # ----------------------
 
 
+def write_storm_file(rundata):
+    """Write the storm file configured in rundata.surge_data.
+
+    Separated from setgeo so that importing/calling setrun() does not write
+    files as a side effect (important for test runners that call setrun() and
+    then substitute their own storm file).
+    """
+    data = rundata.surge_data
+
+    isaac = Storm()
+    atcf_path = (Path(__file__).parent / "bal092012.dat").resolve()
+    isaac.read(path=atcf_path, file_format="ATCF")
+    # Calculate landfall time - Need to specify as the file above does not
+    # include this info (~2345 UTC - 6:45 p.m. CDT - on August 28)
+    isaac.time_offset = np.datetime64("2012-08-29T00:00")
+
+    if data.storm_specification_type == 'holland80':
+        isaac.write(data.storm_file, file_format='geoclaw')
+
+    elif data.storm_specification_type == "data":
+        # ASCII
+        isaac.file_format = 'NWS12'
+        isaac.file_paths = [
+            (Path(__file__).parent / "isaac.PRE").resolve(),
+            (Path(__file__).parent / "isaac.WIN").resolve(),
+        ]
+
+        # NetCDF
+        # isaac.file_format = "NWS13"
+        # isaac.file_paths.append((Path() / "isaac.nc").resolve())
+        # # For NetCDF we need to set this to the epoch or things do not work yet
+        # isaac.time_offset = np.datetime64("1970-01-01")
+
+        # Stretch (>1) or compress (<1) the storm time axis relative to the
+        # first time step.  1.0 = no change; 2.0 = storm moves twice as slow.
+        isaac.storm_time_scale = 1.0
+
+        isaac.write(data.storm_file, file_format='data')
+
+    else:
+        raise ValueError("Invalid storm specification type.")
+
+
 if __name__ == '__main__':
     # Set up run-time parameters and write all data files.
     import sys
@@ -471,4 +489,5 @@ def setgeo(rundata):
     else:
         rundata = setrun()
 
+    write_storm_file(rundata)
     rundata.write()