5. pySALEplotのサンプル

pySALEPlotの解説&実習 20161109pySALEPlot.pdf の実践編に掲載されてあるサンプル集

DenTmp.py

密度と温度のプロット

import pySALEPlot as psp
import matplotlib.pyplot as plt
from numpy import arange,sqrt,ma

# Need this for the colorbars we will make on the mirrored plot
from mpl_toolkits.axes_grid1 import make_axes_locatable

# This example plotting script designed to plot 
# material and temperature in the Chicxulub example

# If viridis colormap is available, use it here
try:
    plt.set_cmap('viridis')
except:
    plt.set_cmap('YlGnBu_r')

# distances between tracers
def get_distances(s,line):
    x=s.xmark[line]
    y=s.ymark[line]
    return sqrt((x[:-1]-x[1:])**2+(y[:-1]-y[1:])**2)

# Define the maximum separation allowed when plotting lines
maxsep=3.

def make_colorbars(ax, p, f, units):
    # Create axes either side of the plot to place the colorbars
    divider = make_axes_locatable(ax)
    cx1 = divider.append_axes("right", size="5%", pad=0.4)
    cx2 = divider.append_axes("left", size="5%", pad=0.8)
    cb1 = fig.colorbar(p[0], cax=cx1)
    cb1.set_label(psp.longFieldName(f[0])+units[0])
    cb2 = fig.colorbar(p[1], cax=cx2)
    cb2.set_label(psp.longFieldName(f[1])+units[1])
    # Need to set the labels on the left for this colorbar
    cx2.yaxis.tick_left()
    cx2.yaxis.set_label_position('left')

# Make an output directory
dirname='DenTmp'
psp.mkdir_p(dirname)

# Open the datafile
model=psp.opendatfile('jdata.dat')

# Set the distance units to km
model.setScale('km')

# Set up a pylab figure
fig=plt.figure(figsize=(8,4))

ax=fig.add_subplot(111,aspect='equal')

# Loop over timesteps
for i in arange(0,model.nsteps,2):

    # Set the axis labels
    ax.set_xlabel('r [km]')
    ax.set_ylabel('z [km]')

    # Set the axis limits
    ax.set_xlim([-100,100])
    ax.set_ylim([-50,45])

    # Read the time step 'i' from the datafile:
    # read two or more fields by making a list of their abbreviations
    step=model.readStep(['Tmp', 'Den'],i)

    # -- the classic iSALEPlot mirrored setup. Plot the second field
    # -- using negative x values
    p1=ax.pcolormesh(model.x,model.y,step.Den,
            cmap='Oranges') #,vmin=1,vmax=model.nmat+1)
    p2=ax.pcolormesh(-model.x,model.y,step.data[0],
            vmin=0,vmax=1500)

    # Material boundaries
    [ax.contour(model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]
    [ax.contour(-model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]

    # Tracer lines
    for u in range(1,model.tracer_numu):
        tru=model.tru[u]

        # Plot the tracers in horizontal lines, every 5 lines
        for l in arange(0,len(tru.xlines),5):

            # Get the distances between pairs of tracers in xlines
            dist=get_distances(step,tru.xlines[l])
            # Mask the xmark values if separation too big... means the line won't be connected here
            ax.plot(ma.masked_array(step.xmark[tru.xlines[l]][:-1],mask=dist > maxsep*tru.d[0]),
                    step.ymark[tru.xlines[l]][:-1],
                    c='#808080',marker='None',linestyle='-',linewidth=0.5)

    # Add colorbars; only need to do this once
    if i == 0:
        make_colorbars(ax, [p1, p2], [step.plottype[1], step.plottype[0]], [r'[kg/m$^3$]', '[K]'])
    ax.set_title('{: 5.2f} s'.format(step.time))

    # Save the figure
    fig.savefig('{}/DenTmp{:05d}.png'.format(dirname,i),dpi=300)

    # Remove the field, ready for the next timestep to be plotted
    ax.cla()

profile_example1.py

プロットするX座標の変更

import pySALEPlot as psp
import matplotlib.pyplot as plt

# This example plotting script designed to plot 
# a vertical pressure profile of the Chicxulub example

# Open the datafile
model=psp.opendatfile('jdata.dat')

# Set the distance units to mm
model.setScale('km')

# Set up a pyplot figure
fig=plt.figure()
ax=fig.add_subplot(111)

# Set the axis labels
ax.set_xlabel('Pressure [GPa]')
ax.set_ylabel('Depth [km]')

# Set the axis limits
ax.set_ylim(-200,0)
#ax.set_ylim(0,100)

# Read the time steps from the datafile
step0=model.readStep('Pre',0)
step1=model.readStep('Pre',1)

# Plot the pressure field
ax.plot(step1.Pre[0,:]*1e-9,model.yc[0,:],'k:',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[0,:][0]))
ax.plot(step1.Pre[50,:]*1e-9,model.yc[0,:],'m--',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[50,:][0]))
ax.plot(step1.Pre[100,:]*1e-9,model.yc[0,:],'b-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[100,:][0]))
ax.plot(step1.Pre[200,:]*1e-9,model.yc[0,:],'r-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[200,:][0]))

# Add a legend
ax.legend(loc=0)

# Save the figure
fig.savefig('profile1.png')

profile_example2.py

トレーサー粒子の最大圧力(TrP)をプロット

import pySALEPlot as psp
import matplotlib.pyplot as plt

# This example plotting script designed to plot 
# a vertical pressure profile of the Chicxulub example

# Open the datafile
model=psp.opendatfile('jdata.dat')

# Set the distance units to mm
model.setScale('km')

# Set up a pyplot figure
fig=plt.figure()
ax=fig.add_subplot(111)

# Set the axis labels
ax.set_xlabel('r [km]')
ax.set_ylabel('z [km]')

# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)

# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)

# Plot the pressure field
for u in range(model.tracer_numu): #loop over tracer clouds
    tstart = model.tru[u].start
    tend = model.tru[u].end
    scat = ax.scatter(step.xmark[tstart:tend],step.ymark[tstart:tend],
               c=step1.TrP[angle_range]*1e-9,vmin=0,vmax=30,
               cmap='viridis',s=4,linewidths=0) 
cb=fig.colorbar(scat)
cb.set_label('Tracer Prak Pressure [GPa]')

# Save the figure
fig.savefig('profile2.png')

profile_example3.py

ある角度でのトレーサー粒子の最大圧力(TrP)をプロット

import pySALEPlot as psp
import matplotlib.pyplot as plt
import numpy as np

# This example plotting script designed to plot 
# a vertical pressure profile of the Chicxulub example

# Open the datafile
model=psp.opendatfile('jdata.dat')

# Set the distance units to mm
model.setScale('km')

# Set up a pyplot figure
fig=plt.figure()
ax=fig.add_subplot(111)

# Set the axis labels
ax.set_xlabel('r [km]')
ax.set_ylabel('z [km]')

# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)

# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)

# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark)) 
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))
# Plot the pressure field 
scat = ax.scatter(step0.xmark[angle_range],step0.ymark[angle_range],
                  c=step1.TrP[angle_range]*1e-9,vmin=0,vmax=30,
                  cmap='viridis',s=4,linewidths=0)
cb=fig.colorbar(scat)
cb.set_label('Tracer Prak Pressure [GPa]')

# Save the figure
fig.savefig('profile3.png')

profile_example4.py

ある角度でのトレーサー粒子の最大圧力(TrP)と距離をプロット

import pySALEPlot as psp
import matplotlib.pyplot as plt

# This example plotting script designed to plot 
# a vertical pressure profile of the Chicxulub example

# Open the datafile
model=psp.opendatfile('jdata.dat')

# Set the distance units to mm
model.setScale('km')

# Set up a pyplot figure
fig=plt.figure()
ax=fig.add_subplot(111)

# Set the axis labels
ax.set_xlabel('Distance from impact point [km]')
ax.set_ylabel('Peak Pressure [GPa]')

# Set the axis limits
ax.set_xlim(0,30)
ax.set_ylim(0,60)

# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)

# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark)) 
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))

r = np.sqrt(step0.xmark[angle_range]**2+step0.ymark[angle_range]**2)
trp = step1.TrP[angle_range]*1e-9
ax.plot(np.sort(r),trp[np.argsort(r)],'k-')

# Save the figure
fig.savefig('profile4.png')
Last modified:2016/11/21 08:12:24
Keyword(s):
References:[SideMenu]