Feature/update upstream (#38)

* change nodata_value to -99999 in topotools

Previous value -9999 did not have enough digits since new etopo 2022
data has this as an actual topography value in the Mariana Trench
(when written with `Z_format='%.0f'`).

* Fix bug related to tracking pressure at gauges

* simplify make_fgout_animation.py use of update_artists

These do not need to be passed into update, unpacked, and repacked,
since the objects created in the script will be used into update.
If blit==False then they are not needed at all.

* use image backend Agg in make_fgout_animation.py so animation size agrees with specified figure size

* set blit=False and auto-detect number of fgout frames

* clean up make_fgout_animation.py

* added new chile2010_fgmax-fgout/make_fgout_animation_with_transect.py

* Add dZ_format parameter to DTopography.write function

The default is now '%.3f', millimeter resolution, making smaller dtopo
files than previously.

* handle 0 radius

* fix type on set_pressure

* fix duplicate r

* make comparison fp-safe

* handle underflow

* avoid underflow in wind setting too

* point to riemann/src for Riemann solvers

* Remove Riemann solvers from src/2d/bouss, Makefile.common points to riemann/src

Note that rpn2_geoclaw.f and geoclaw_riemann_utils.f were update in riemann/src
to handle 5 waves and rpt2_geoclaw_sym.f was discarded in favor of
standard rpt2.f, which works fine with 5 waves.

* Initial CI GitHub action script

* Update testing.yml

* Add checking out of clawpack and geoclaw

* Add lint and testing to CI

* Disable linting for the time being

* Re-enable linting but only for geoclaw

* Exclude old_topotools.py

Do not want to touch the old topotools file for reference.

* Fix up gauge plotting for storm surge

Includes a couple of minor bugs related to gauge plotting.  Major change
involves how we now plot the gauge data.

* Add dry gauge plotting

* create gauge filenames that allow more than 5 digits in gauge number

If fewer than 5 digits, still add zero padding e.g. gauge00012.txt
but no longer throws an error for large gauge numbers, e.g. gauge1234567.txt

* cleaner way to zero pad only if fewer than 5 digits using I0.5 format

* Cleanup and reimplement reading of ATCF and writing of GeoClaw storms

* Minor string comparison bug fix

* Extract name of storm from file name

* Minor tweaks and fixes

* Fix missing comma

* AutoPEP8 surge code

* add fgout output_style parameter and support for array of output_times

* remove trailing whitespaces in fgout_module.f990

* Refactor fgout_module.f90 so it works for either GeoClaw or D-Claw

and support in fgout_tools.py for new dclaw attribute dclaw to set in setrun.py
to indicate D-Claw, in which case 7 components of q are output instead of 4.
Support for eventually indicating fewer components to output.

* Use geoclaw module rho

* Remove module level parameters that were not needed or were specific

These mostly pertained to the CLE code.

* Initial implementation of rotation control

* Minor bugfixes and rearranging

* Fix bugs in the non-spherical coordinates for storms"

* fix fgout_tools.FGoutGrid.read_fgout_grids_data for time array

Now that `output_style==2` is supported for fgout grids (an array of times,
see clawpack#617), the function
fgout_tools.FGoutGrid.read_fgout_grids_data needs to be fixed to properly
read in the new format of `fgout_grids.data`, also note that `nout`
now comes before `tstart` and `tend` for `output_style==1`.

---------

Co-authored-by: Randy LeVeque <rjl@uw.edu>
Co-authored-by: Kyle Mandli <kyle.Mandli@gmail.com>

.github/workflows/testing.yml

@@ -40,12 +40,7 @@ jobs:
     - name: Setup geoclaw
       run: |
         cd geoclaw
-        if [[ "${{ github.event_name }}" == "pull_request" ]]; then
-          git fetch origin pull/${{ github.event.pull_request.number }}/merge:PR
-          git checkout PR
-        else
-          git checkout ${{ github.ref_name }}
-        fi
+        git checkout ${{ github.ref }}
     - name: Lint with flake8
       run: |
         cd geoclaw

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 \