Skip to content

Commit

Permalink
Merge pull request #626 from rjleveque/chile2010_setplot
Browse files Browse the repository at this point in the history
update examples/tsunami/chile2010/setplot.py
  • Loading branch information
rjleveque authored Aug 25, 2024
2 parents 88cdbc7 + afc7ff9 commit 3ac6822
Showing 1 changed file with 33 additions and 34 deletions.
67 changes: 33 additions & 34 deletions examples/tsunami/chile2010/setplot.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,18 +57,14 @@ def addgauges(current_data):

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes('pcolor')
plotaxes.title = 'Surface'
plotaxes.scaled = True

def fixup(current_data):
import pylab
addgauges(current_data)
t = current_data.t
t = t / 3600. # hours
pylab.title('Surface at %4.2f hours' % t, fontsize=20)
pylab.xticks(fontsize=15)
pylab.yticks(fontsize=15)
plotaxes.afteraxes = fixup
plotaxes.title = 'Surface at time h:m:s'
plotaxes.aspect_latitude = -30.
plotaxes.xticks_fontsize = 10
plotaxes.yticks_fontsize = 10
plotaxes.xlabel = 'longitude'
plotaxes.ylabel = 'latitude'

plotaxes.afteraxes = addgauges

# Water
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
Expand All @@ -78,6 +74,8 @@ def fixup(current_data):
plotitem.pcolor_cmin = -0.2
plotitem.pcolor_cmax = 0.2
plotitem.add_colorbar = True
plotitem.colorbar_shrink = 0.8
plotitem.colorbar_extend = 'both'
plotitem.amr_celledges_show = [0,0,0]
plotitem.patchedges_show = 1

Expand Down Expand Up @@ -114,18 +112,37 @@ def fixup(current_data):

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = 'auto'
plotaxes.ylimits = 'auto'
plotaxes.time_scale = 1/3600. # convert to hours
plotaxes.xlimits = [0, 9]
plotaxes.ylimits = [-0.3, 0.3]
plotaxes.title = 'Surface'
plotaxes.title_fontsize = 15
plotaxes.time_label = 'time (hours)'
plotaxes.ylabel = 'surface elevation (m)'
plotaxes.grid = True

def add_obs(current_data):
from pylab import plot, legend
gaugeno = current_data.gaugeno
if gaugeno == 32412:
try:
tgauge = TG32412[:,0] / 3600. # convert to hours
plot(tgauge, TG32412[:,1], 'r', linewidth=1.0)
legend(['GeoClaw','Obs'],loc='lower right')
except: pass

plotaxes.afteraxes = add_obs

# Plot surface as blue curve:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = 3
plotitem.plot_var = 3 # eta
plotitem.plotstyle = 'b-'
plotitem.kwargs = {'linewidth':1.8}


# Plot topo as green curve:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.show = False
plotitem.show = False # not being used

def gaugetopo(current_data):
q = current_data.q
Expand All @@ -137,24 +154,6 @@ def gaugetopo(current_data):
plotitem.plot_var = gaugetopo
plotitem.plotstyle = 'g-'

def add_zeroline(current_data):
from pylab import plot, legend, xticks, floor, axis, xlabel
t = current_data.t
gaugeno = current_data.gaugeno

if gaugeno == 32412:
try:
plot(TG32412[:,0], TG32412[:,1], 'r')
legend(['GeoClaw','Obs'],loc='lower right')
except: pass
axis((0,t.max(),-0.3,0.3))

plot(t, 0*t, 'k')
n = int(floor(t.max()/3600.) + 2)
xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)])
xlabel('time (hours)')

plotaxes.afteraxes = add_zeroline


#-----------------------------------------
Expand Down

0 comments on commit 3ac6822

Please sign in to comment.