Merge remote-tracking branch 'clawpack/master' into switch-pytest

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/2d/bouss/Makefile.bouss

@@ -167,6 +167,6 @@ COMMON_SOURCES += \
   $(BOUSSLIB)/resetBoussStuff.f \
   $(BOUSSLIB)/simpleBound.f90 \
   $(BOUSSLIB)/umfpack_support.f \
-  $(BOUSSLIB)/rpn2_geoclaw.f \
-  $(BOUSSLIB)/rpt2_geoclaw_sym.f \
-  $(BOUSSLIB)/geoclaw_riemann_utils.f \
+  $(CLAW)/riemann/src/rpn2_geoclaw.f \
+  $(CLAW)/riemann/src/rpt2_geoclaw.f \
+  $(CLAW)/riemann/src/geoclaw_riemann_utils.f \