Merge pull request clawpack#612 from mandli/update-surge-gauge-plot

Fix up gauge plotting for storm surge

examples/storm-surge/ike/setplot.py

@@ -155,38 +155,23 @@ def friction_after_axes(cd):
     plotfigure.show = True
     plotfigure.clf_each_gauge = True
 
-    # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes()
-    plotaxes.xlimits = [-2, 1]
-    # plotaxes.xlabel = "Days from landfall"
-    # plotaxes.ylabel = "Surface (m)"
-    plotaxes.ylimits = [-1, 5]
-    plotaxes.title = 'Surface'
-
-    def gauge_afteraxes(cd):
-
-        axes = plt.gca()
-        surgeplot.plot_landfall_gauge(cd.gaugesoln, axes)
-
-        # Fix up plot - in particular fix time labels
-        axes.set_title('Station %s' % cd.gaugeno)
-        axes.set_xlabel('Days relative to landfall')
-        axes.set_ylabel('Surface (m)')
-        axes.set_xlim([-2, 1])
-        axes.set_ylim([-1, 5])
-        axes.set_xticks([-2, -1, 0, 1])
-        axes.set_xticklabels([r"$-2$", r"$-1$", r"$0$", r"$1$"])
-        axes.grid(True)
-    plotaxes.afteraxes = gauge_afteraxes
-
-    # Plot surface as blue curve:
+    plotaxes.time_scale = 1 / (24 * 60**2)
+    plotaxes.grid = True
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = "Surface"
+    plotaxes.ylabel = "Surface (m)"
+    plotaxes.time_label = "Days relative to landfall"
+
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
-    # plotitem.plot_var = 3
-    # plotitem.plotstyle = 'b-'
-
-    #
+    plotitem.plot_var = surgeplot.gauge_surface
+    # Plot red area if gauge is dry
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.plot_var = surgeplot.gauge_dry_regions
+    plotitem.kwargs = {"color":'lightcoral', "linewidth":5}
+
     #  Gauge Location Plot
-    #
     def gauge_location_afteraxes(cd):
         plt.subplots_adjust(left=0.12, bottom=0.06, right=0.97, top=0.97)
         surge_afteraxes(cd)

examples/storm-surge/isaac/setplot.py

@@ -1,20 +1,18 @@
-
-from __future__ import absolute_import
-from __future__ import print_function
+#!/usr/bin/env python
 
 import os
+import warnings
+import datetime
 
-import numpy
+import numpy as np
 import matplotlib.pyplot as plt
-import datetime
 
 import clawpack.visclaw.colormaps as colormap
 import clawpack.visclaw.gaugetools as gaugetools
 import clawpack.clawutil.data as clawutil
 import clawpack.amrclaw.data as amrclaw
 import clawpack.geoclaw.data as geodata
-
-
+import clawpack.geoclaw.util as geoutil
 import clawpack.geoclaw.surge.plot as surgeplot
 
 try:
@@ -88,9 +86,9 @@ def friction_after_axes(cd):
     regions = {"Gulf": {"xlimits": (clawdata.lower[0], clawdata.upper[0]),
                         "ylimits": (clawdata.lower[1], clawdata.upper[1]),
                         "figsize": (6.4, 4.8)},
-                "Louisiana": {"xlimits": (-92, -83),
-                               "ylimits": (27.5, 30.5),
-                               "figsize": (8, 2.7)}}
+               "Louisiana": {"xlimits": (-92, -83),
+                             "ylimits": (27.5, 30.5),
+                             "figsize": (8, 2.7)}}
 
     for (name, region_dict) in regions.items():
 
@@ -175,40 +173,82 @@ def friction_after_axes(cd):
     # ========================================================================
     #  Figures for gauges
     # ========================================================================
+    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,
+                                                                    begin_date, 
+                                                                    end_date)
+        if water_level is None:
+            print("*** Could not fetch gauge {}.".format(station_id))
+        else:
+            # Convert to seconds relative to landfall
+            t = (date_time - landfall_time) / np.timedelta64(1, 's')
+            t /= (24 * 60**2)
+
+            # Detide
+            water_level -= tide
+
+            # Plot data
+            ax = plt.gca()
+            ax.plot(t, water_level, color='lightgray', marker='x')
+            ax.set_title(station_name)
+            ax.legend(['Computed', "Observed"])
+
+
     plotfigure = plotdata.new_plotfigure(name='Gauge Surfaces', figno=300,
                                          type='each_gauge')
     plotfigure.show = True
     plotfigure.clf_each_gauge = True
 
-    # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes()
-    plotaxes.xlimits = [-2, 1]
-    # plotaxes.xlabel = "Days from landfall"
-    # plotaxes.ylabel = "Surface (m)"
-    plotaxes.ylimits = [-1, 5]
-    plotaxes.title = 'Surface'
-
-    def gauge_afteraxes(cd):
-
-        axes = plt.gca()
-        landfall = 0.
-        surgeplot.plot_landfall_gauge(cd.gaugesoln, axes, landfall=landfall)
-
-        # Fix up plot - in particular fix time labels
-        axes.set_title('Station %s' % cd.gaugeno)
-        axes.set_xlabel('Days relative to landfall')
-        axes.set_ylabel('Surface (m)')
-        axes.set_xlim([-2, 1])
-        axes.set_ylim([-1, 5])
-        axes.set_xticks([-2, -1, 0, 1])
-        axes.set_xticklabels([r"$-2$", r"$-1$", r"$0$", r"$1$"])
-        axes.grid(True)
-    plotaxes.afteraxes = gauge_afteraxes
-
-    # Plot surface as blue curve:
+    plotaxes.time_scale = 1 / (24 * 60**2)
+    plotaxes.grid = True
+    plotaxes.xlimits = [-2, 1.5]
+    plotaxes.ylimits = [-0.25, 1]
+    plotaxes.title = "Surface"
+    plotaxes.ylabel = "Surface (m)"
+    plotaxes.time_label = "Days relative to landfall"
+    plotaxes.afteraxes = plot_observed
+
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.plot_var = surgeplot.gauge_surface
+    # Plot red area if gauge is dry
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
-    plotitem.plot_var = 3
-    plotitem.plotstyle = 'b-'
+    plotitem.plot_var = surgeplot.gauge_dry_regions
+    plotitem.kwargs = {"color":'lightcoral', "linewidth":5}
+
+    #  Gauge Location Plot
+    def gauge_location_afteraxes(cd):
+        plt.subplots_adjust(left=0.12, bottom=0.06, right=0.97, top=0.97)
+        surge_afteraxes(cd)
+        gaugetools.plot_gauge_locations(cd.plotdata, gaugenos='all',
+                                        format_string='ko', add_labels=True)
+
+    plotfigure = plotdata.new_plotfigure(name="Gauge Locations")
+    plotfigure.show = True
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.title = 'Gauge Locations'
+    plotaxes.scaled = True
+    plotaxes.xlimits = regions['Louisiana']['xlimits']
+    plotaxes.ylimits = regions['Louisiana']['ylimits']
+    plotaxes.afteraxes = gauge_location_afteraxes
+    surgeplot.add_surface_elevation(plotaxes, bounds=surface_limits)
+    surgeplot.add_land(plotaxes, bounds=[0.0, 20.0])
+    plotaxes.plotitem_dict['surface'].amr_patchedges_show = [0] * 10
+    plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10
 
     # -----------------------------------------
     # Parameters used only when creating html and/or latex hardcopy
@@ -217,7 +257,7 @@ def gauge_afteraxes(cd):
     plotdata.printfigs = True                # print figures
     plotdata.print_format = 'png'            # file format
     plotdata.print_framenos = 'all'          # list of frames to print
-    plotdata.print_gaugenos = [1, 2, 3, 4]   # list of gauges to print
+    plotdata.print_gaugenos = 'all'          # list of gauges to print
     plotdata.print_fignos = 'all'            # list of figures to print
     plotdata.html = True                     # create html files of plots?
     plotdata.latex = True                    # create latex file of plots?

examples/storm-surge/isaac/setrun.py

@@ -313,17 +313,16 @@ def setrun(claw_pkg='geoclaw'):
     regions = rundata.regiondata.regions
     # to specify regions of refinement append lines of the form
     #  [minlevel,maxlevel,t1,t2,x1,x2,y1,y2]
-    # Gauges from Ike AWR paper (2011 Dawson et al)
-    rundata.gaugedata.gauges.append([1, -95.04, 29.07,
-                                     rundata.clawdata.t0,
-                                     rundata.clawdata.tfinal])
-    rundata.gaugedata.gauges.append([2, -94.71, 29.28,
-                                     rundata.clawdata.t0,
-                                     rundata.clawdata.tfinal])
-    rundata.gaugedata.gauges.append([3, -94.39, 29.49,
+
+    # Pilots Station East, S.W. Pass, LA - 28°55.9429'N, 89°24.4445'W
+    # https://tidesandcurrents.noaa.gov/stationhome.html?id=8760922
+    rundata.gaugedata.gauges.append([1, -89.40740833, 28.93238167,
                                      rundata.clawdata.t0,
                                      rundata.clawdata.tfinal])
-    rundata.gaugedata.gauges.append([4, -94.13, 29.58,
+
+    # Grand Isle, LA - 29°15.8761'N 89°57.4960'W
+    # https://tidesandcurrents.noaa.gov/stationhome.html?id=8761724
+    rundata.gaugedata.gauges.append([2, -89.95826667, 29.26460167,
                                      rundata.clawdata.t0,
                                      rundata.clawdata.tfinal])
 

src/python/geoclaw/surge/plot.py

@@ -83,20 +83,25 @@ def gauge_locations(current_data, gaugenos='all'):
                                     add_labels=True, xoffset=0.02,
                                     yoffset=0.02)
 
+def gauge_dry_regions(cd, dry_tolerance=1e-16):
+    """Masked array of zeros where gauge is dry."""
+    return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) > dry_tolerance, 
+                              np.zeros(cd.gaugesoln.q[0,:].shape))
 
-def gaugetopo(current_data):
-    q = current_data.q
-    h = q[0, :]
-    eta = q[3, :]
-    topo = eta - h
-    return topo
+def gauge_surface(cd, dry_tolerance=1e-16):
+    """Sea surface at gauge masked when dry."""
+    return np.ma.masked_where(np.abs(cd.gaugesoln.q[0, :]) < dry_tolerance, 
+                              cd.gaugesoln.q[3, :])
 
 
 def plot_landfall_gauge(gauge, axes, landfall=0.0, style='b', kwargs={}):
     """Plot gauge data on the axes provided
 
     This will transform the plot so that it is relative to the landfall value
     provided.
+
+    This can be done using `plotaxes.time_scale` instead so this function will
+    be deprecated and removed in a future release.
     """
     axes = plt.gca()
 

src/python/geoclaw/util.py

@@ -361,7 +361,8 @@ def get_cache_path(product):
                                end_date.strftime(cache_date_fmt))
         filename = '{}_{}_{}'.format(time_zone, datum, units)
         abs_cache_dir = os.path.abspath(cache_dir)
-        return os.path.join(abs_cache_dir, product, station, dates, filename)
+        return os.path.join(abs_cache_dir, product, str(station), dates, 
+                            filename)
 
     def save_to_cache(cache_path, data):
         # make parent directories if they do not exist