diff --git a/examples/storm-surge/ike/setplot.py b/examples/storm-surge/ike/setplot.py index 5784f9ce5..865093be4 100644 --- a/examples/storm-surge/ike/setplot.py +++ b/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) diff --git a/examples/storm-surge/isaac/setplot.py b/examples/storm-surge/isaac/setplot.py index 5fb7aebc3..8fe715b15 100644 --- a/examples/storm-surge/isaac/setplot.py +++ b/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? diff --git a/examples/storm-surge/isaac/setrun.py b/examples/storm-surge/isaac/setrun.py index d6d63eb93..4e8299286 100644 --- a/examples/storm-surge/isaac/setrun.py +++ b/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]) diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py index 254d74f7b..372ca1073 100644 --- a/src/python/geoclaw/surge/plot.py +++ b/src/python/geoclaw/surge/plot.py @@ -83,13 +83,15 @@ 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={}): @@ -97,6 +99,9 @@ def plot_landfall_gauge(gauge, axes, landfall=0.0, style='b', kwargs={}): 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() diff --git a/src/python/geoclaw/util.py b/src/python/geoclaw/util.py index 816043f05..b358c9546 100644 --- a/src/python/geoclaw/util.py +++ b/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