Python ClawPlotData.outdir方法代码示例

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [as 别名]
def read_data(outdir="_output", adjoint=False):
    
    from numpy import loadtxt
    fname = outdir + '/fort.H'
    B = loadtxt(fname)
    print "Loaded B"
    
    pd = ClawPlotData()
    pd.outdir = outdir

    times = []
    qxt = []
    for frameno in range(5001):
        try:
            frame = pd.getframe(frameno)
        except:
            break
        q = frame.state.q
        t = frame.state.t
        q[0,:] = B + q[0,:]
        qxt.append(q)
        times.append(t)
    
    x = frame.state.patch.x.centers
    x = x
    X,T = np.meshgrid(x,times)
    qxt = np.array(qxt)
    if adjoint:
        qxt = np.flipud(qxt)  # reverse t for adjoint
    return X,T,qxt

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [as 别名]
def make_getgauge(outdir='_output'):
    """
    Create function getgauge that will grab one set of gauge data
    from the fort.gauge file in directory specified by outdir.
    """
    from clawpack.visclaw.data import ClawPlotData
    plotdata = ClawPlotData()
    plotdata.outdir = outdir
    getgauge = plotdata.getgauge
    return getgauge

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [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 outdir [as 别名]
def plotclaw(outdir='.', plotdir='_plots', setplot = 'setplot.py',format='ascii'):
    """
    Create html and/or latex versions of plots.

    INPUT:
        setplot is a module containing a function setplot that will be called
                to set various plotting parameters.
        format specifies the format of the files output from Clawpack
    """

    from clawpack.visclaw.data import ClawPlotData
    from clawpack.visclaw import plotpages

    plotdata = ClawPlotData()
    plotdata.outdir = outdir
    plotdata.plotdir = plotdir
    plotdata.setplot = setplot
    plotdata.format = format

    plotpages.plotclaw_driver(plotdata, verbose=False, format=format)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [as 别名]
def read_data(outdir="_output", adjoint=False):
    pd = ClawPlotData()
    pd.outdir = outdir

    times = []
    qxt = []
    for frameno in range(5001):
        try:
            frame = pd.getframe(frameno)
        except:
            break
        q = frame.state.q
        t = frame.state.t
        qxt.append(q)
        times.append(t)
    
    x = frame.state.patch.x.centers
    x = x
    X,T = np.meshgrid(x,times)
    qxt = np.array(qxt)
    if adjoint:
        qxt = np.flipud(qxt)  # reverse t for adjoint
    return X,T,qxt

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [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 outdir [as 别名]
g_obs=tidegauge.read_tide_gauge('1617760__2011-03-11_to_2011-03-13.csv')
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"

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

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

plotdata = ClawPlotData()

if 1:
    plotdata.outdir = '_output'
    toffset = 92.

if 0:
    plotdata.outdir = '_output_manning_0.025'
    toffset = 92.
if 0:
    plotdata.outdir = '_output_manning015_cfl090'
    toffset = 92.
if 0:
    plotdata.outdir = '_output_manning015_cfl089'
    toffset = 96.

figure(50,figsize=(18,12))
clf()

# ---  Gauge 1 ---

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

d = loadtxt('s1v.txt')

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [as 别名]
import matplotlib
matplotlib.use('Agg')
from pylab import *
#from pyclaw.plotters.data import ClawPlotData
from clawpack.visclaw.data import ClawPlotData

Runs = '../Runs/'
pd = ClawPlotData()
pd.outdir = Runs + 'HAI1107/_output'
g1 = pd.getgauge(1)
g1107 = pd.getgauge(1107)
g12340 = pd.getgauge(12340)


figure(1)
clf()
t = g1.t / 3600.
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')

legend(loc='upper right')
title('Surface elevation at Honolulu gauges')



figure(2)
clf()

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [as 别名]
    matplotlib  # see if it's already been imported (interactive session)
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)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [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)

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [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 outdir [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 outdir [as 别名]
def plotclaw(outdir='.', plotdir='_plots', setplot = 'setplot.py',
             format='ascii', msgfile='', frames=None, verbose=False):
    """
    Create html and/or latex versions of plots.

    INPUT:
        setplot is a module containing a function setplot that will be called
                to set various plotting parameters.
        format specifies the format of the files output from Clawpack
    """

    from clawpack.visclaw.data import ClawPlotData
    from clawpack.visclaw import plotpages

    plotdata = ClawPlotData()
    plotdata.outdir = outdir
    plotdata.plotdir = plotdir
    plotdata.setplot = setplot
    plotdata.format = format
    plotdata.msgfile = msgfile


    frametools.call_setplot(plotdata.setplot, plotdata)


    if plotdata.num_procs is None:
        plotdata.num_procs = int(os.environ.get("OMP_NUM_THREADS", 1))

    # Make sure plotdata.parallel is False in some cases:


    if plotdata.parallel:
        assert type(setplot) in [str, bool, type(None)], \
                "*** Parallel plotting is not supported when ClawPlotData " \
                + "attribute setplot is a function."

    if plotdata.parallel and (plotdata.num_procs > 1):

        # If this is the original call then we need to split up the work and 
        # call this function again

        # First set up plotdir:

        plotdata._parallel_todo = 'initialize'
        plotpages.plotclaw_driver(plotdata, verbose=False, format=format)

        if frames is None:
            if plotdata.num_procs is None:
                plotdata.num_procs = int(os.environ.get("OMP_NUM_THREADS", 1))


            frames = [[] for n in xrange(plotdata.num_procs)]
            framenos = frametools.only_most_recent(plotdata.print_framenos,
                                                   outdir)

            # don't use more procs than frames or infinite loop!!
            num_procs = min(plotdata.num_procs, len(framenos))

            for (n, frame) in enumerate(framenos):
                frames[n%num_procs].append(frame)

            # Create subprocesses to work on each
            plotclaw_cmd = "python %s" % __file__
            process_queue = []
            for n in xrange(num_procs):
                plot_cmd = "%s %s %s %s" % (plotclaw_cmd, 
                                            outdir,
                                            plotdir, 
                                            setplot)
                plot_cmd = plot_cmd + " " + " ".join([str(i) for i in frames[n]])

                process_queue.append(subprocess.Popen(plot_cmd, shell=True))

            poll_interval = 5
            try:
                while len(process_queue) > 0:
                    time.sleep(poll_interval)
                    for process in process_queue:
                        if process.poll() is not None:
                            process_queue.remove(process)
                    if verbose:
                        print "Number of processes currently:",len(process_queue)
            
            # Stop child processes if interrupt was caught or something went 
            # wrong
            except KeyboardInterrupt:
                print "ABORTING: A keyboard interrupt was caught.  All " + \
                      "child processes will be terminated as well."
                for process in process_queue:
                    process.terminate()
                raise

            except:
                print "ERROR: An error occurred while waiting for " + \
                      "plotting processes to complete.  Aborting all " + \
                      "child processes."
                for process in process_queue:
                    process.terminate()
                raise 

#.........这里部分代码省略.........

# 需要导入模块: from clawpack.visclaw.data import ClawPlotData [as 别名]
# 或者: from clawpack.visclaw.data.ClawPlotData import outdir [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