Python ClawPlotData.getgauge方法代码示例

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
def process_gauge(gaugeNo, output):
    plotdata = ClawPlotData()
    plotdata.outdir = '_output'   # set to the proper output directory
    g = plotdata.getgauge(gaugeNo)
    # g.t is the array of times
    # g.q is the array of values recorded at the gauges (g.q[m,n] is the m`th variable at time `t[n])
    t = g.t
    h = g.q[0,:]
    u = np.where(h>0.001,g.q[1,:]/h,0.)
    v = np.where(h>0.001,g.q[2,:]/h,0.)
    h = h - h[:10].sum()/10.
    hu2 = h*u**2
    hv2 = h*v**2
    header = "Time (s), water depth (m), u-velocity (m/s), v-velocity (m/s), hu^2 (m^3/s^2), hv^2 (m^3/s^2)"
    if 'data_at_gauges_txt' not in os.listdir('./'):
        os.mkdir('./data_at_gauges_txt')
    result = np.vstack((t, h, u, v, hu2, hv2))
    result = result.transpose()
    np.savetxt('./data_at_gauges_txt/'+output, result, delimiter=',', header=header)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
def make_output_for_dakota():
    from clawpack.visclaw.data import ClawPlotData
    plotdata = ClawPlotData()
    plotdata.outdir = os.path.join(os.getcwd() , '_output')   # set to the proper output directory
    gaugeno = 34                  # gauge number to examine
    g = plotdata.getgauge(gaugeno)
    gauge_max = g.q[3,:].max()
    print "Maximum elevation observed at gauge %s: %6.2f meters" \
        % (gaugeno, gauge_max)
    fname = 'results.out'
    f = open(fname,'w')
    f.write("%6.2f\n" % -gauge_max)
    f.close()
    print "Created ",fname

    if 1:
        figure()
        plot(g.t, g.q[3,:])
        title("Gauge %s" % gaugeno)
        ylabel('Elevation (meters)')
        xlabel('time (seconds)')
        fname = 'gauge.png'
        savefig(fname)
        print 'Created ',fname

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
from pylab import *
from scipy import interpolate

from clawpack.visclaw.data import ClawPlotData

plotdata = ClawPlotData()

plotdata.outdir = '_output_1-3sec_alltime'

#tfinal = 4.9 * 3600.
tfinal = 6.4 * 3600.  # for alltime
dt = 1.  # time increment for output files
tout = arange(0., tfinal, dt)   

g = plotdata.getgauge(3333)
p = interpolate.interp1d(g.t, g.q[3,:]) # interpolate surface
g3333_eta = p(tout)

g = plotdata.getgauge(7761)
p = interpolate.interp1d(g.t, g.q[3,:]) # interpolate surface
g7761_eta = p(tout)

g = plotdata.getgauge(1125)
u = g.q[1,:]/g.q[0,:]
v = g.q[2,:]/g.q[0,:]
s = sqrt(u**2 + v**2)
p = interpolate.interp1d(g.t, s) # interpolate speed
g1125_speed = p(tout)

g = plotdata.getgauge(1126)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
    fname = "BM1_leveque_4.txt"

figure(figsize=(8, 12))
clf()

# ---  Gauge 1 ---

d = loadtxt("s1u.txt")
t1u = d[:, 0]
s1u = d[:, 1]

d = loadtxt("s1v.txt")
t1v = d[:, 0]
s1v = d[:, 1]

g = plotdata.getgauge(1)
u1 = g.q[1, :] / g.q[0, :]
v1 = g.q[2, :] / g.q[0, :]
t1 = g.t

subplot(4, 1, 1)
plot(t1u + toffset, s1u, "b", label="Experiment")
plot(t1, u1, "r", label="GeoClaw")
ylabel("u (m/s)")
legend(loc="upper right")
title(plotdata.outdir)

subplot(4, 1, 2)
plot(t1v + toffset, s1v, "b", label="Experiment")
plot(t1, v1, "r", label="GeoClaw")
ylabel("v (m/s)")

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
outdir = '../Runs/HAI1125-26/_output'

tsec = g_obs[1]
thour = tsec / 3600.
eta = g_obs[3]

figure(3)
clf()
plot(thour,eta,'k',linewidth=1)
plot(thour,eta,'k.',linewidth=1,label='Observed')

# computed results:
plotdata = ClawPlotData()
plotdata.outdir = outdir
print "Looking for GeoClaw results in ",plotdata.outdir
g7760 = plotdata.getgauge(7760)

# shift by 10 minutes:
thour = (g7760.t + 600.) / 3600.
plot(thour, g7760.q[3,:],'r',linewidth=2,label='GeoClaw')
xlim(7.5,13)
ylim(-2,1.5)
legend(loc='lower right')
xticks(range(8,14),fontsize=15)
yticks(fontsize=15)
title('Surface elevation at Gauge 7760',fontsize=15)

if 1:
    fname = "TG7760.jpg"
    savefig(fname)
    print "Created ",fname

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
except:
    import matplotlib

    matplotlib.use("Agg")  # set to image backend

import os
from pylab import *
from clawpack.visclaw.data import ClawPlotData

rundir = "../Runs/HAI1107"
outdir = os.path.join(rundir, "_output")
obsdir = os.path.abspath("../Observations")

pd = ClawPlotData()
pd.outdir = outdir
g1 = pd.getgauge(1)
g1107 = pd.getgauge(1107)
g12340 = pd.getgauge(12340)


figure(1)
clf()
t = g1.t / 3600.0
plot(t, g1.q[3, :], "b", linewidth=1, label="Gauge S1")
plot(t, g1107.q[3, :], "r", linewidth=2, label="HAI1107")
# plot(t, g12340.q[3,:], 'g', linewidth=2, label='TG 2340')
xlabel("Hours post-quake")
ylabel("meters")
ylim(-1, 1)
xticks(fontsize=15)
yticks(fontsize=15)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
def compare_gauges(outdir1, outdir2, gaugenos='all', q_components='all',
                    tol=0., verbose=True, plot=False):

    """
    Compare gauge output in two output directories.

    :Input:
     - *outdir1, outdir2* -- output directories
     - *gaugenos* -- list of gauge numbers to compare, or 'all' in which case
       outdir1/gauges.data will be used to determine gauge numbers.
     - *q_components* -- list of components of q to compare.
     - *tol* -- tolerance for checking equality
     - *verbose* -- print out dt and dq for each comparison?
     - *plot* -- if True, will produce a plot for each gauge, with a
       subfigure for each component of q.


    Returns True if everything matches to tolerance *tol*, else False.
    """

    from clawpack.visclaw.data import ClawPlotData
    from matplotlib import pyplot as plt
    
    if gaugenos == 'all':
        # attempt to read from gauges.data:
        try:
            setgauges1 = read_setgauges(outdir1)
            setgauges2 = read_setgauges(outdir2)
        except:
            print '*** could not read gauges.data from one of the outdirs'
            return
        gaugenos = setgauges1.gauge_numbers
        if setgauges2.gauge_numbers != gaugenos:
            print '*** warning -- outdirs have different sets of gauges'

        if len(gaugenos)==0:
            print "*** No gauges found in gauges.data"
            return

    plotdata1 = ClawPlotData()
    plotdata1.outdir = outdir1
    plotdata2 = ClawPlotData()
    plotdata2.outdir = outdir2

    matches = True
    for gaugeno in gaugenos:
        g1 = plotdata1.getgauge(gaugeno,verbose=verbose)
        t1 = g1.t
        q1 = g1.q

        g2 = plotdata2.getgauge(gaugeno,verbose=verbose)
        t2 = g2.t
        q2 = g2.q

        dt = abs(t1-t2).max()
        if verbose:
            print  "Max difference in t[:] at gauge %s is %g" % (gaugeno,dt)
        matches = matches and (dt <= tol)

        if q_components == 'all':
            q_components = range(q1.shape[0])

        for m in q_components:
            dq = abs(q1[m,:]-q2[m,:]).max()
            if verbose:
                print  "Max difference in q[%s] at gauge %s is %g" % (m,gaugeno,dq)
            matches = matches and (dq <= tol)

        if plot:
            plt.figure(gaugeno)
            plt.clf()
            mq = len(q_components)
            for k,m in enumerate(q_components):
                plt.subplot(mq,1,k+1)
                plt.plot(g1.t,g1.q[m,:],'b',label='outdir1')
                plt.plot(g2.t,g2.q[m,:],'r',label='outdir2')
                plt.legend()
                plt.title('q[%s] at gauge number %s' % (m,gaugeno))
        
    return matches       

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
# Folder from out or out2 differ in that the Riemann solver used uses the Lagrangian transformation
# on all of water or only in the interface respectively

# CHOOSE out or outmap in the outdir to choose from non-mapped or mapped version of code
outstr = '_outlim'

#plotdata.outdir = '../../_output'   # set to the proper output directory
#g = plotdata.getgauge(gaugeno)
#p = zeros(g.t.size)
#p = (gammawat - 1.0)*(g.q[3,:] - 0.5*(g.q[1,:]*g.q[1,:] + g.q[2,:]*g.q[2,:])/g.q[0,:]) - gammawat*pinfwat
#p = 0.001*p # Convert to KPa
#tt = g.t*1000000 # Convert to microsec
#plot(tt, p, '-g', label="Level New", linewidth=3)

plotdata.outdir = outstr +'_conv_40x20_lvl6_refrat2-2-2-2-2'   # set to the proper output directory
g = plotdata.getgauge(gaugeno)
p = zeros(g.t.size)
p = (gammawat - 1.0)*(g.q[3,:] - 0.5*(g.q[1,:]*g.q[1,:] + g.q[2,:]*g.q[2,:])/g.q[0,:]) - gammawat*pinfwat
p = 0.001*p # Convert to KPa
tt = g.t*1000000 # Convert to microsec
plot(tt, p, '-k', label="Level 6", linewidth=1)

plotdata.outdir = outstr +'_conv_40x20_lvl5_refrat2-2-2-2'   # set to the proper output directory
g = plotdata.getgauge(gaugeno)
p = zeros(g.t.size)
p = (gammawat - 1.0)*(g.q[3,:] - 0.5*(g.q[1,:]*g.q[1,:] + g.q[2,:]*g.q[2,:])/g.q[0,:]) - gammawat*pinfwat
p = 0.001*p # Convert to KPa
tt = g.t*1000000 # Convert to microsec
plot(tt, p, '-b', label="Level 5", linewidth=1)

# Load and plot ngauge from 40x20 grid AMR run with level 3 AMR and refrinement ratio [2,2]

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
    def check_gauges(self, save=False, 
                           regression_data_path="regression_data.txt",
                           tolerance=1e-14):
        r"""Basic test to assert gauge equality

        Test all gauges found in gauges.data.  Do full comparison of all
        times, levels, components of q.

        :Input:
         - *save* (bool) - If *True* will save the output from this test to 
           the file *regresion_data.txt*.  Default is *False*.
         - *regression_data_path* (path) - Path to the regression test data.
           Defaults to 'regression_data.txt'.
         - *tolerance* (float) - Tolerance used in comparison, defaults to
           *1e-14*.
        """
        from clawpack.visclaw import gaugetools
        from clawpack.visclaw.data import ClawPlotData

        # Get gauge data
        plotdata = ClawPlotData()
        plotdata.outdir = self.temp_path

        setgauges = gaugetools.read_setgauges(plotdata.outdir)
        gauge_numbers = setgauges.gauge_numbers
        
        # read in all gauge data and sort by gauge numbers so we
        # can properly compare.  Note gauges may be written to fort.gauge
        # in random order when OpenMP used.

        for gaugeno in gauge_numbers:
            g = plotdata.getgauge(gaugeno)
            m = len(g.level)
            number = numpy.array(m*[g.number])
            level = numpy.array(g.level)
            data_gauge = numpy.vstack((number,level,g.t,g.q)).T
            try:    
                data = numpy.vstack((data, data_gauge))
            except:
                data = data_gauge  # first time thru loop

        # save the gauge number sorted by gaugeno in case user wants to
        # compare if assertion fails.
        sorted_gauge_file = 'test_gauge.txt'
        sorted_gauge_path = os.path.join(self.temp_path, sorted_gauge_file)
        numpy.savetxt(sorted_gauge_path, data)

        # Get (and save) regression comparison data
        regression_data_file = os.path.join(self.test_path,
                                            regression_data_path)
        if save:
            numpy.savetxt(regression_data_file, data)
            print "Saved new regression data to %s" % regression_data_file
        regression_data = numpy.loadtxt(regression_data_file)

        # if assertion fails, indicate to user what files to compare:
        output_dir = os.path.join(os.getcwd(),
                                     "%s_output" % self.__class__.__name__)
        sorted_gauge_path = os.path.join(output_dir, sorted_gauge_file)
        error_msg = "Full gauge match failed.  Compare these files: " + \
               "\n  %s\n  %s" % (regression_data_file, sorted_gauge_path) + \
               "\nColumns are gaugeno, level, t, q[0:num_eqn]"

        assert numpy.allclose(data, regression_data, tolerance), error_msg

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
b1a = reshape(b1, (9000, 4))
b4 = loadtxt(datadir + "/B4.txt", skiprows=1)
b4a = reshape(b4, (9000, 4))
b6 = loadtxt(datadir + "/B6.txt", skiprows=1)
b6a = reshape(b6, (9000, 4))
b9 = loadtxt(datadir + "/B9.txt", skiprows=1)
b9a = reshape(b9, (9000, 4))


plotdata = ClawPlotData()
plotdata.outdir = "_output_3"

figure(50, figsize=(8, 12))
clf()
for gnum, wg in zip([1, 2, 3, 4], [wg1, wg2, wg3, wg4]):
    g = plotdata.getgauge(gnum)
    subplot(4, 1, gnum)
    plot(t, wg, "b", label="Measured")
    plot(g.t, g.q[3, :], "r", label="GeoClaw")
    xlim(0, 40)
    title("Gauge %s" % gnum)
    ylabel("surface (m)")
legend(loc="upper left")


if 0:
    gauges = {}
    gauges = [
        0,
        1,
        2,

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
g_obs = numpy.loadtxt(os.path.join(TGdir, 'TG_1612340_detided.txt'))

tsec = g_obs[:,0]
thour = tsec / 3600.
eta = g_obs[:,1]

plt.figure(3)
plt.clf()
plt.plot(thour,eta,'k',linewidth=1)
plt.plot(thour,eta,'k.',linewidth=1,label='Observed')

# computed results:
plotdata = ClawPlotData()
plotdata.outdir = outdir
print "Looking for GeoClaw results in ",plotdata.outdir
g = plotdata.getgauge(12340)

# shift by 10 minutes:
thour = (g.t + 600.) / 3600.
plt.plot(thour, g.q[3,:],'r',linewidth=2,label='GeoClaw')
plt.xlim(7,13)
plt.ylim(-1,1)
plt.legend(loc='lower right')
plt.xticks(range(8,14),fontsize=15)
plt.yticks(fontsize=15)
plt.xlabel('Hours post-quake')
plt.ylabel('meters')
plt.title('Surface elevation at TG 1612340',fontsize=15)

if save_plot:
    fname = "../Figures/TG_1612340_compare.png"

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
plotdata = ClawPlotData()

plotdata.outdir = '_output'
#toffset = 7. # hours
toffset = 6.6 # hours fits data at 7760 better

figure(3333, figsize=(12,5))
clf()

# ---  Specified control point ---

CP = loadtxt('../se.dat')
CP_t = CP[:,0] / 60.  # hours since EQ
CP_eta = CP[:,1] + 0.13 + sea_level # adjust to match geoclaw run

g3333 = plotdata.getgauge(3333)
t = g3333.t / 3600. + 7.


plot(CP_t,CP_eta,'k-o',label='Specified')
plot(t, g3333.q[3,:], 'r',label='GeoClaw')
ylabel('meters')
legend(loc='upper right')
xlim(7.0,12.5)
ylim(-1,1.5)
title('Surface elevation at control point')
show()


# ---  Tide Gauge 7760 ---

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
#plotdata.outdir = '_output_1_3sec'
plotdata.outdir = '_output'
toffset = 0.2

port_data = loadtxt(comparison_data_dir + 'port_data.txt')
tAB = port_data[:,0] - toffset
zAB = port_data[:,1] # tsunami only
tTug = port_data[:,3] - toffset
zTug = port_data[:,4] # tsunami only

figure(1, figsize=(12,5))
clf()


g = plotdata.getgauge(1)
t = g.t / 3600.

plot(tAB,zAB,'k-o',label='Observations')
plot(t, g.q[3,:], 'r',lw=2,label='GeoClaw')
ylabel('meters')
legend(loc='upper right')
xlim(12.5, 18)
ylim(-1,1)
#title('Surface elevation at A Beacon')
show()

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

figure(2, figsize=(12,5))
clf()

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
g_obs = numpy.loadtxt(os.path.join(TGdir, 'TG_1615680_detided.txt'))

tsec = g_obs[:,0]
thour = tsec / 3600.
eta = g_obs[:,1]

plt.figure(3)
plt.clf()
plt.plot(thour,eta,'k',linewidth=1)
plt.plot(thour,eta,'k.',linewidth=1,label='Observed')

# computed results:
plotdata = ClawPlotData()
plotdata.outdir = outdir
print "Looking for GeoClaw results in ",plotdata.outdir
g = plotdata.getgauge(5680)

# shift by 10 minutes:
thour = (g.t + 600.) / 3600.
plt.plot(thour, g.q[3,:],'r',linewidth=2,label='GeoClaw')
plt.xlim(7.5,13)
plt.ylim(-3,3)
plt.legend(loc='lower right')
plt.xticks(range(8,14),fontsize=15)
plt.yticks(fontsize=15)
plt.xlabel('Hours post-quake')
plt.ylabel('meters')
plt.title('Surface elevation at TG 1615680',fontsize=15)

if save_plot:
    fname = "../Figures/TG_1615680_compare.png"

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import getgauge [as 别名]
def make_figs(gaugeno,rundir,outdir,veldir):

    # measurement values:
    #mdir = os.path.abspath(os.path.join(obsdir, veldir))
    mdir = os.path.join(obsdir, veldir)
    print "Looking for observations in ",mdir
    
    # Use detided values:
    fname = 'detided_harmonic.txt'
    t_m,u_m,v_m = loadtxt(os.path.join(mdir,fname), unpack=True)
    print "Read detided u,v from ",fname


    # computed results:
    plotdata = ClawPlotData()
    plotdata.outdir = os.path.join(rundir,outdir)
    print "Looking for GeoClaw results in ",plotdata.outdir

    # from PGV.plot_gauge:
    g = plotdata.getgauge(gaugeno)
    gh = g.q[0,:]
    ghu = g.q[1,:]
    ghv = g.q[2,:]
    t = g.t + 600.  # Add 10 minutes for better agreement
    u = where(gh>0.01, ghu/gh * 100., 0.)  # convert to cm/sec
    v = where(gh>0.01, ghv/gh * 100., 0.)  # convert to cm/sec


    #t,speed,u,v,eta = PGV.plot_gauge(gaugeno,plotdata)
    t = t/3600.  # convert to hours

    i_ts = find((t_m >= t.min()) & (t_m <= t.min() + 5.))

    figure(101)
    clf()
    #plot(u_m,v_m,'k')
    plot(u,v,'r',linewidth=3,label='GeoClaw')
    plot(u_m[i_ts],v_m[i_ts],'ko',label='Observed')

    smax = max(abs(u).max(),abs(v).max(),abs(u_m).max(),abs(v_m).max()) * 1.05
    plot([-smax,smax],[0,0],'k')
    plot([0,0],[-smax,smax],'k')
    axis('scaled')
    axis([-smax,smax,-smax,smax])
    legend(loc=('lower right'))
    title('Velocities at gauge %s' % gaugeno)
    #xlabel('u in cm/sec')
    #ylabel('v in cm/sec')
    xlabel('u (east-west velocity) in cm/sec')
    ylabel('v (north-south velocity) in cm/sec')
    fname = "../Figures/figure%sa.png" % gaugeno
    savefig(fname)
    print "Created ",fname

    figure(102)  #,figsize=(8,10))
    clf()
    subplot(2,1,1)
    plot(t_m,u_m,'ko-',linewidth=2)
    plot(t,u,'r',linewidth=3)
    #legend(['Observed','GeoClaw'],'upper left')
    title('u (east-west velocity) at gauge %s' % gaugeno)
    xlim(t.min(),t.max())
    ylim(-smax,smax)
    ylabel('cm / sec')
    #ylim(-80,40)

    subplot(2,1,2)
    plot(t_m,v_m,'ko-',linewidth=2,label="Observed")
    #plot(t_m[i_ts],v_m[i_ts],'bo-',linewidth=2)
    plot(t,v,'r',linewidth=3,label="GeoClaw")
    #legend(loc=('lower left'))
    title('v (north-south velocity) at gauge %s' % gaugeno)
    
    xlim(t.min(),t.max())
    ylim(-smax,smax)
    ylabel('cm/sec')
    xlabel('Hours post-quake')
 
    fname = "../Figures/figure%sb.png" % gaugeno
    savefig(fname)
    print "Created ",fname