Skip to content

Commit

Permalink
Update Katrina setplot
Browse files Browse the repository at this point in the history
  • Loading branch information
@mandli
mandli committed Sep 17, 2024
1 parent 3f24587 commit cd2f018
Showing 1 changed file with 91 additions and 94 deletions.
185 changes: 91 additions & 94 deletions katrina_2005/setplot.py
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@

"""
"""
Set up the plot figures, axes, and items to be done for each frame.
This module is imported by the plotting routines and then the
function setplot is called to set the plot parameters.
"""

import os
Expand All @@ -19,7 +19,7 @@
import clawpack.amrclaw.data
import clawpack.geoclaw.data

import clawpack.geoclaw.surge.plot as surge
import clawpack.geoclaw.surge.plot as surgeplot

try:
from setplotfg import setplotfg
Expand All @@ -28,7 +28,7 @@

def setplot(plotdata):
r"""Setplot function for surge plotting"""


plotdata.clearfigures() # clear any old figures,axes,items data
plotdata.format = 'binary'
Expand All @@ -44,11 +44,11 @@ def setplot(plotdata):
friction_data.read(os.path.join(plotdata.outdir,'friction.data'))

# Load storm track32
track = surge.track_data(os.path.join(plotdata.outdir,'fort.track'))
track = surgeplot.track_data(os.path.join(plotdata.outdir,'fort.track'))

# Set afteraxes function
def surge_afteraxes(cd):
return surge.surge_afteraxes(cd, track, plot_direction=False)
return surgeplot.surge_afteraxes(cd, track, plot_direction=False)

# Limits for plots
dx = 0.5
Expand Down Expand Up @@ -104,9 +104,9 @@ def surge_afteraxes(cd):
plotaxes.ylimits = ylimits
plotaxes.afteraxes = surge_afteraxes

surge.add_surface_elevation(plotaxes, bounds=surface_limits)
surgeplot.add_surface_elevation(plotaxes, bounds=surface_limits)
plotaxes.plotitem_dict['surface'].amr_patchedges_show = [0] * 10
surge.add_land(plotaxes)
surgeplot.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10

# ======================================================================
Expand All @@ -123,124 +123,120 @@ def surge_afteraxes(cd):
plotaxes.ylimits = ylimits
plotaxes.afteraxes = surge_afteraxes

surge.add_speed(plotaxes, bounds=speed_limits)
surgeplot.add_speed(plotaxes, bounds=speed_limits)
plotaxes.plotitem_dict['speed'].amr_patchedges_show = [0] * 10
surge.add_land(plotaxes)
surgeplot.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10

# ======================================================================
# Wind Field
# ======================================================================
plotfigure = plotdata.new_plotfigure(name='Wind Speed - %s' % name)
plotfigure.show = surge_data.wind_forcing

plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Wind Field"
plotaxes.scaled = True
plotaxes.xlimits = xlimits
plotaxes.ylimits = ylimits
plotaxes.afteraxes = surge_afteraxes
# Friction field
plotfigure = plotdata.new_plotfigure(name='Friction')
plotfigure.show = friction_data.variable_friction

surge.add_wind(plotaxes, bounds=wind_limits)
plotaxes.plotitem_dict['wind'].amr_patchedges_show = [0] * 10
surge.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Manning's N Coefficients"
plotaxes.scaled = True
plotaxes.xlimits = regions[0]["limits"][0]
plotaxes.ylimits = regions[0]["limits"][1]
plotaxes.afteraxes = surge_afteraxes

surgeplot.add_friction(plotaxes, bounds=friction_bounds)

# ========================================================================
# Hurricane forcing
# ========================================================================
# Friction field
plotfigure = plotdata.new_plotfigure(name='Friction - %s' % name)
plotfigure.show = friction_data.variable_friction
# Wind Field
plotfigure = plotdata.new_plotfigure(name='Wind Speed')
plotfigure.show = surge_data.wind_forcing

plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Manning's N Coefficients"
plotaxes.scaled = True
plotaxes.xlimits = xlimits
plotaxes.ylimits = ylimits
plotaxes.afteraxes = surge_afteraxes
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Wind Field"
plotaxes.scaled = True
plotaxes.xlimits = regions[0]["limits"][0]
plotaxes.ylimits = regions[0]["limits"][1]
plotaxes.afteraxes = surge_afteraxes

surge.add_friction(plotaxes, bounds=friction_bounds)
surgeplot.add_wind(plotaxes, bounds=wind_limits)
plotaxes.plotitem_dict['wind'].amr_patchedges_show = [0] * 10
surgeplot.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10

# Pressure field
plotfigure = plotdata.new_plotfigure(name='Pressure - %s' % name)
plotfigure.show = surge_data.pressure_forcing

plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Pressure Field"
plotaxes.scaled = True
plotaxes.xlimits = xlimits
plotaxes.ylimits = ylimits
plotaxes.afteraxes = surge_afteraxes

surge.add_pressure(plotaxes, bounds=pressure_limits)
plotaxes.plotitem_dict['pressure'].amr_patchedges_show = [0] * 10
surge.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10
# Pressure field
plotfigure = plotdata.new_plotfigure(name='Pressure')
plotfigure.show = surge_data.pressure_forcing

plotaxes = plotfigure.new_plotaxes()
plotaxes.title = "Pressure Field"
plotaxes.scaled = True
plotaxes.xlimits = regions[0]["limits"][0]
plotaxes.ylimits = regions[0]["limits"][1]
plotaxes.afteraxes = surge_afteraxes

surgeplot.add_pressure(plotaxes, bounds=pressure_limits)
plotaxes.plotitem_dict['pressure'].amr_patchedges_show = [0] * 10
surgeplot.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10


# ========================================================================
# Figures for gauges
# ========================================================================
plotfigure = plotdata.new_plotfigure(name='Surface & Topo', figno=300, \
type='each_gauge')
plotfigure.show = True
plotfigure.clf_each_gauge = True

stations = [('8761724', 'Grand Isle, LA'),
('8760922', 'Pilots Station East, SW Pass, LA'),
('8735180', 'Dauphin Island, AL')]

landfall_time = np.datetime64('2005-08-29T11:10')
begin_date = datetime.datetime(2005, 8, 26)
end_date = datetime.datetime(2005, 8, 31)

def get_actual_water_levels(station_id):
# Fetch water levels and tide predictions for given station
date_time, water_level, tide = fetch_noaa_tide_data(station_id,
begin_date, end_date)

# Calculate times relative to landfall
secs_rel_landfall = (date_time - landfall_time) / np.timedelta64(1, 's')
days_rel_landfall = (date_time - landfall_time) / np.timedelta64(1,'s') / 86400

# Subtract tide predictions from measured water levels
water_level -= tide

return secs_rel_landfall, water_level
# Find the mean values
# Data imported every 6 minutes (i.e. 360 seconds)
num_data_pts = (end_date - begin_date).total_seconds() / 360
mean_value = np.sum(water_level) / num_data_pts
# water_level -= mean_value
print(f"{station_id}: {mean_value}")

return days_rel_landfall, water_level

def gauge_afteraxes(cd):
station_id, station_name = stations[cd.gaugeno - 1]
secs_rel_landfall, actual_level = get_actual_water_levels(station_id)
def plot_observed(cd):
station_id, station_name = stations[cd.gaugeno-1]
days_rel_landfall, actual_level = get_actual_water_levels(station_id)

axes = plt.gca()
axes.plot(secs_rel_landfall, actual_level, 'g')
ax = plt.gca()
ax.plot(days_rel_landfall, actual_level, 'g')

# Fix up plot - in particular fix time labels
axes.set_title(station_name)
axes.set_xlabel('Days relative to landfall')
axes.set_ylabel('Surface (m)')
axes.set_xlim(np.array([-3, 1]) * 86400)
axes.set_ylim([-0.5, 2.5])
axes.set_xticks(np.linspace(-3, 1, 5) * 86400)
axes.set_xticklabels([r"$-3$", r"$-2$", r"$-1$", r"$0$", r"$1$"])
axes.grid(True)

# Plot wind speed using second scale
wind_speed = np.linalg.norm(cd.gaugesoln.q[8:10, :], axis=0)
axes2 = axes.twinx()
axes2.plot(cd.gaugesoln.t, wind_speed, 'r.', markersize=0.3)
axes2.set_ylabel('Wind Speed (m/s)')
axes2.set_ylim([-2.5, 52.5])
plotfigure = plotdata.new_plotfigure(name='Surface & Topo', figno=300, \
type='each_gauge')
plotfigure.show = True
plotfigure.clf_each_gauge = True

stations = [('8761724', 'Grand Isle, LA'),
('8760922', 'Pilots Station East, SW Pass, LA'),
('8735180', 'Dauphin Island, AL')]

landfall_time = np.datetime64('2005-08-29T11:10')
begin_date = datetime.datetime(2005, 8, 26)
end_date = datetime.datetime(2005, 8, 31)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.afteraxes = gauge_afteraxes
plotaxes.time_scale = 1 / (24 * 60**2)
plotaxes.xlimits = [-3, 1]
plotaxes.ylimits = [-0.5, 2.5]
plotaxes.title = "Surface"
plotaxes.ylabel = "Surface (m)"
plotaxes.time_label = "Days relative to landfall"
plotaxes.afteraxes = plot_observed
plotaxes.grid = True

# Plot surface as blue curve:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = 3
plotitem.plotstyle = 'b-'

plotitem.plot_var = surgeplot.gauge_surface
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = surgeplot.gauge_dry_regions
plotitem.kwargs = {"color":'lightcoral', "linewidth":5}

#
# Gauge Location Plot
Expand All @@ -262,14 +258,14 @@ def gauge_location_afteraxes(cd):
plotaxes.ylimits = [28.0, 31.0]
plotaxes.afteraxes = gauge_location_afteraxes

surge.add_surface_elevation(plotaxes, bounds=surface_limits)
surgeplot.add_surface_elevation(plotaxes, bounds=surface_limits)
plotaxes.plotitem_dict['surface'].amr_patchedges_show = [0] * 10
surge.add_land(plotaxes)
surgeplot.add_land(plotaxes)
plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10


#-----------------------------------------

# Parameters used only when creating html and/or latex hardcopy
# e.g., via pyclaw.plotters.frametools.printframes:

Expand All @@ -284,6 +280,7 @@ def gauge_location_afteraxes(cd):
plotdata.latex_figsperline = 2 # layout of plots
plotdata.latex_framesperline = 1 # layout of plots
plotdata.latex_makepdf = False # also run pdflatex?
plotdata.parallel = True # parallel plotting

return plotdata

0 comments on commit cd2f018

Please sign in to comment.