Python data.ClawPlotData类代码示例

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

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

def setplot(plotdata=None):
#--------------------------
    
    """ 
    Specify what is to be plotted at each frame.
    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
    Output: a modified version of plotdata.
    
    """ 

    if plotdata is None:
        from clawpack.visclaw.data import ClawPlotData
        plotdata = ClawPlotData()

    plotdata.clearfigures()

    # Figures corresponding to Figure 9.5 of LeVeque, "Finite Volume
    # Methods for Hyperbolic Problems," 2002 (though more of them)

    # Tuples of (variable name, variable number)
    figdata = [('Pressure', 0),
               ('Velocity', 1)]

    # Afteraxes function: draw a vertical dashed line at the interface
    # between different media
    def draw_interface(current_data):
        import pylab
        pylab.plot([0., 0.], [-1000., 1000.], 'k--')

    for varname, varid in figdata:
        plotfigure = plotdata.new_plotfigure(name=varname, figno=varid)

        plotaxes = plotfigure.new_plotaxes()
        plotaxes.xlimits = [-5., 5.]
        plotaxes.ylimits = [-0.5, 1.5]    # Good for both vars because of near-unit impedance
        plotaxes.title = varname
        plotaxes.afteraxes = draw_interface

        plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
        plotitem.plot_var = varid
        plotitem.color = 'b'

    plotdata.printfigs = True          # Whether to output figures
    plotdata.print_format = 'png'      # What type of output format
    plotdata.print_framenos = 'all'    # Which frames to output
    plotdata.print_fignos = 'all'      # Which figures to print
    plotdata.html = True               # Whether to create HTML files
    plotdata.latex = False             # Whether to make LaTeX output

    return plotdata

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)

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)

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

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

    rundata.clawdata.num_cells[0] = mx
    rundata.clawdata.num_output_times = mframe	# output_times=1 won't work for high grid resolution
    rundata.clawdata.tfinal = 0.100000e+00
    rundata.clawdata.lower[0] = xlower
    rundata.clawdata.upper[0] = xupper 
    rundata.clawdata.order = 1
    rundata.clawdata.bc_lower[0] = 'periodic'#'user'   # at xlower
    rundata.clawdata.bc_upper[0] = 'periodic'#'extrap'   # at xupper
    rundata.write()
    runclaw(xclawcmd='xclaw',outdir=outdir)	# xclaw.exe file produced after make .exe 

    #Get the material parameters
    aux = np.loadtxt(outdir+'/fort.a0000',skiprows=5)	# don't delate skiprows or set it equal 6 

    plotdata = ClawPlotData()
    plotdata.outdir=outdir

    #Read in the solution
    dat = plotdata.getframe(mframe)
    u = dat.q[0,:]

#    print u

    stress = np.exp(u*aux) - 1	# not sure why here don't need aux[0,:]
    stress = u

    #Compute parameter for error calculation
    reshape_para = mx_exact/mx

    #Compute error

try:
    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)

def setplot(plotdata=None):
#--------------------------
    
    """ 
    Specify what is to be plotted at each frame.
    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
    Output: a modified version of plotdata.
    
    """ 

    from clawpack.visclaw import colormaps

    if plotdata is None:
        from clawpack.visclaw.data import ClawPlotData
        plotdata = ClawPlotData()


    plotdata.clearfigures()  # clear any old figures,axes,items data
    

    # Figure for density - pcolor
    plotfigure = plotdata.new_plotfigure(name='Density', figno=0)

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.xlimits = [0,1]
    plotaxes.ylimits = [0,1]
    plotaxes.title = 'Density'
    plotaxes.scaled = True
    plotaxes.afteraxes = addgauges

    # Set up for item on these axes:
    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
    plotitem.plot_var = 0
    #plotitem.pcolor_cmap = colormaps.yellow_red_blue
    plotitem.pcolor_cmin = 0.
    plotitem.pcolor_cmax = 2.
    plotitem.add_colorbar = True
    plotitem.amr_patchedges_show = [0]
    plotitem.amr_celledges_show = [0]


    # Figure for density - Schlieren
    plotfigure = plotdata.new_plotfigure(name='Schlieren', figno=1)

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.xlimits = [0,1]
    plotaxes.ylimits = [0,1]
    plotaxes.title = 'Density'
    plotaxes.scaled = True      # so aspect ratio is 1

    # Set up for item on these axes:
    plotitem = plotaxes.new_plotitem(plot_type='2d_schlieren')
    plotitem.schlieren_cmin = 0.0
    plotitem.schlieren_cmax = 1.0
    plotitem.plot_var = 0
    plotitem.add_colorbar = False


    # Figure for grid cells
    plotfigure = plotdata.new_plotfigure(name='cells', figno=2)

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.xlimits = [0,1]
    plotaxes.ylimits = [0,1]
    plotaxes.title = 'Grid patches'
    plotaxes.scaled = True

    # Set up for item on these axes:
    plotitem = plotaxes.new_plotitem(plot_type='2d_patch')
    plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee']
    plotitem.amr_celledges_show = [1,0]
    plotitem.amr_patchedges_show = [1]


    #-----------------------------------------
    # Figures for gauges
    #-----------------------------------------
    plotfigure = plotdata.new_plotfigure(name='q', figno=300, \
                    type='each_gauge')
    plotfigure.clf_each_gauge = True

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.xlimits = [0,1]
    plotaxes.ylimits = [0,1]
    plotaxes.title = 'Density'

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


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

    plotdata.printfigs = True                # print figures
#.........这里部分代码省略.........

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       

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 pylab import *
from clawpack.visclaw.data import ClawPlotData
from scipy import interpolate

# Aux Parameters
gammawat = 7.15
pinfwat = 300000000.0
rhow = 1000.0
gaugeno = 1 #1 and 2

plotdata = ClawPlotData()

# 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

def setplot(plotdata=None):
#--------------------------
    
    """ 
    Specify what is to be plotted at each frame.
    Input:  plotdata, an instance of pyclaw.plotters.data.ClawPlotData.
    Output: a modified version of plotdata.
    
    """ 

    import os

    import numpy as np
    import matplotlib.pyplot as plt

    from clawpack.visclaw import geoplot, gaugetools, colormaps

    import clawpack.clawutil.data as clawutil
    import clawpack.amrclaw.data as amrclaw
    import clawpack.geoclaw.data

    import clawpack.geoclaw.multilayer.plot as ml_plot

    if plotdata is None:
        from clawpack.visclaw.data import ClawPlotData
        plotdata = ClawPlotData()

    from numpy import linspace



    plotdata.clearfigures()  # clear any old figures,axes,items data
    plotdata.save_frames = False

    # Load data from output
    clawdata = clawutil.ClawInputData(2)
    clawdata.read(os.path.join(plotdata.outdir,'claw.data'))
    amrdata = amrclaw.AmrclawInputData(clawdata)
    amrdata.read(os.path.join(plotdata.outdir,'amr.data'))
    geodata = clawpack.geoclaw.data.GeoClawData()
    geodata.read(os.path.join(plotdata.outdir,'geoclaw.data'))
    multilayer_data = clawpack.geoclaw.data.MultilayerData()
    multilayer_data.read(os.path.join(plotdata.outdir,'multilayer.data'))

    # To plot gauge locations on pcolor or contour plot, use this as
    # an afteraxis function:

    def addgauges(current_data):
        from clawpack.visclaw import gaugetools
        gaugetools.plot_gauge_locations(current_data.plotdata, \
             gaugenos='all', format_string='ko', add_labels=True)

    # ========================================================================
    #  Generic helper functions
    # ========================================================================
    def pcolor_afteraxes(current_data):
        # bathy_ref_lines(current_data)
        gauge_locations(current_data)
        
    def contour_afteraxes(current_data):
        # gauge_locations(current_data)
        # m_to_km_labels()
        plt.hold(True)
        pos = -80.0 * (23e3 / 180) + 500e3 - 5e3
        plt.plot([pos,pos],[-300e3,300e3],'b',[pos-5e3,pos-5e3],[-300e3,300e3],'y')
        plt.hold(False)
        wind_contours(current_data)
        bathy_ref_lines(current_data)
        
    def profile_afteraxes(current_data):
        pass
        
    def gauge_locations(current_data,gaugenos='all'):
        plt.hold(True)
        gaugetools.plot_gauge_locations(current_data.plotdata, \
             gaugenos=gaugenos, format_string='kx', add_labels=True)
        plt.hold(False)
    
    # ========================================================================
    # Axis limits
    # xlimits = [amrdata.xlower,amrdata.xupper]
    xlimits = [-100.0, 100.0]

    # ylimits = [amrdata.ylower,amrdata.yupper]
    ylimits = [-100.0, 100.0]
    
    eta = [multilayer_data.eta[0],multilayer_data.eta[1]]

    top_surface_limits = [eta[0]-10,eta[0]+10]
    internal_surface_limits = [eta[1]-5,eta[1]+5]
    depth_limits = [0.0, 0.4]
    top_speed_limits = [0.0,0.1]
    internal_speed_limits = [0.0,0.03]
    

    # ========================================================================
    #  Surface Elevations
    # ========================================================================
    plotfigure = plotdata.new_plotfigure(name='Surface', figno=0)
    plotfigure.show = True
#.........这里部分代码省略.........

def setplot(plotdata=None):
#--------------------------
    
    """ 
    Specify what is to be plotted at each frame.
    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
    Output: a modified version of plotdata.
    
    """ 

    if plotdata is None:
        from clawpack.visclaw.data import ClawPlotData
        plotdata = ClawPlotData()

    plotdata.clearfigures()  # clear any old figures,axes,items data

    # Figure for pressure and velocity:
    plotfigure = plotdata.new_plotfigure(name='Pressure and Velocity', figno=1)

    # Pressure:
    # ---------

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.axescmd = 'subplot(2,1,1)'   # top figure
    plotaxes.xlimits = 'auto'
    plotaxes.ylimits = [-.5,1.]
    plotaxes.title = 'Pressure'

    # Set up for item on these axes:
    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
    plotitem.plot_var = 0
    plotitem.plotstyle = '-'
    plotitem.color = 'b'


    # Velocity:
    # ---------

    # Set up for axes in this figure:
    plotaxes = plotfigure.new_plotaxes()
    plotaxes.axescmd = 'subplot(2,1,2)'   # bottom figure
    plotaxes.xlimits = 'auto'
    plotaxes.ylimits = [-1.,1.]
    plotaxes.title = 'Velocity'

    # Set up for item on these axes:
    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
    plotitem.plot_var = 1
    plotitem.plotstyle = '-'
    plotitem.color = 'b'

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

    plotdata.printfigs = True                # print figures
    plotdata.print_format = 'png'            # file format
    plotdata.print_framenos = 'all'          # list of frames to print
    plotdata.print_fignos = 'all'            # list of figures to print
    plotdata.html = True                     # create html files of plots?
    plotdata.html_homelink = '../README.html'
    plotdata.latex = True                    # create latex file of plots?
    plotdata.latex_figsperline = 2           # layout of plots
    plotdata.latex_framesperline = 1         # layout of plots
    plotdata.latex_makepdf = False           # also run pdflatex?

    return plotdata