[cig-commits] r5335 - mc/3D/CitcomS/trunk/visual
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Nov 21 12:18:03 PST 2006
Author: tan2
Date: 2006-11-21 12:18:02 -0800 (Tue, 21 Nov 2006)
New Revision: 5335
Added:
mc/3D/CitcomS/trunk/visual/plot_annulus.py
mc/3D/CitcomS/trunk/visual/plot_layer.py
Modified:
mc/3D/CitcomS/trunk/visual/zslice.py
Log:
Added scripts for 2D temperature plots, using GMT commands
Added: mc/3D/CitcomS/trunk/visual/plot_annulus.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/plot_annulus.py 2006-11-20 23:49:02 UTC (rev 5334)
+++ mc/3D/CitcomS/trunk/visual/plot_annulus.py 2006-11-21 20:18:02 UTC (rev 5335)
@@ -0,0 +1,254 @@
+#!/usr/bin/env python
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#<LicenseText>
+#
+# plot_annulus.py by Eh Tan.
+# Copyright (C) 2002-2005, California Institute of Technology.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#</LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+"""Plot an annulus cross-section from full spherical CitcomS output
+
+usage: plot_annulus.py modelname step
+ modelname: prefix of CitcomS datafile
+ step: time step to plot
+
+input:
+ modelname.capXX.step
+
+output:
+ modelname.step.z###.tgrd
+ modelname.xsect.step.ps
+"""
+
+
+def get_radius(modelname, step):
+ ### get nodez ###
+ capfile = get_capfile(modelname, 0, step)
+ data = open(capfile)
+ header = data.readline()
+ nodex, nodey, nodez = header.split('x')
+ nodez = int(nodez)
+
+ ### read z coordinate ###
+ radius = range(nodez)
+ for i in range(nodez):
+ radius[i] = float(data.readline().split()[2])
+
+ data.close()
+
+ return radius
+
+
+def get_capfile(modelname, cap, step):
+ return '%s.cap%02d.%d' % (modelname, cap, step)
+
+
+def great_circle_proj():
+ from math import pi, sin, cos, tan, atan
+ while 1:
+ text = '''Choose method to specify great circle path:
+ (1) a center and an azimuth
+ (2) a center and another point on the path
+ (3) a center and a rotation pole position\n'''
+ option = int(raw_input(text))
+
+ if 1 <= option <= 3:
+ lon = float(raw_input('Center longitude: '))
+ lat = float(raw_input('Center latitude: '))
+ center = '-C%f/%f -L-180/180' % (lon, lat)
+ if option == 1:
+ az = float(raw_input('Azimuth (in degrees clockwise from north): '))
+ proj = '%s -A%f' % (center, az)
+ break
+
+ if option == 2:
+ elon = float(raw_input('2nd point longitude: '))
+ elat = float(raw_input('2nd point latitude: '))
+ r2d = 180.0 / pi
+ ## transfer to azimuth mode
+ b = (90 - elat) / r2d
+ a = (90 - lat) / r2d
+ delta = (elon - lon) / r2d
+ if abs(lat) == 90:
+ ## on the pole
+ print 'pole cannot be the center.'
+ continue
+ elif (elon - lon) % 180 == 0:
+ ## on the same meridian
+ az = 0
+ elif lat == 0 and elat == 0:
+ ## on the equator
+ az = 90
+ else:
+ az = atan((sin(a)/tan(b) - cos(a)*cos(delta)) / sin(delta))
+ az = 90 - r2d * az
+
+ proj = '%s -A%f' % (center, az)
+ break
+
+ if option == 3:
+ lon = float(raw_input('Pole longitude: '))
+ lat = float(raw_input('Pole latitude: '))
+ proj = '%s -T%f/%f' % (center, lon, lat)
+ break
+
+ else:
+ print 'Incorrect mode!\n'
+ continue
+
+ return proj
+
+
+#######################################################################
+# Main
+#######################################################################
+
+import os, sys
+import zslice
+
+if len(sys.argv) != 3:
+ print __doc__
+ sys.exit(0)
+
+modelname = sys.argv[1]
+step = int(sys.argv[2])
+
+
+### get radius ###
+radius = get_radius(modelname, step)
+nodez = len(radius)
+#print radius
+
+
+### slice the capfile for all layers ###
+for cap in range(12):
+ capfile = get_capfile(modelname, cap, step)
+ exist = 1
+ for layer in range(nodez):
+ zfile = zslice.zslicefile(capfile, layer)
+ exist = exist and os.path.isfile(zfile)
+
+ if not exist:
+ #print 'Creating zslice files'
+ zslice.zslice(capfile)
+
+
+### create great circle path ###
+gcproj = great_circle_proj()
+gcfile = 'circle.xyp'
+az_resolution = 0.5
+command = 'project %(gcproj)s -G%(az_resolution)f > %(gcfile)s' % vars()
+os.system(command)
+
+
+### create cross section along the great circle ###
+
+## range of layers to plot
+botlayer = 0
+toplayer = nodez - 1
+
+bounds = '0/360/-90/90'
+xsectfile = '%s.xsect.%d.xyz' % (modelname, step)
+out = open(xsectfile,'w')
+
+for layer in range(botlayer, toplayer+1):
+ ## gather the filenames of all zfiles
+ all_zfiles = ''
+ for cap in range(12):
+ capfile = get_capfile(modelname, cap, step)
+ zfile = zslice.zslicefile(capfile, layer)
+ all_zfiles = all_zfiles + ' ' + zfile
+
+ ## create a grdfile for each layer
+ grdfile = '%s.%d.z%03d.tgrd' % (modelname, step, layer)
+ if not os.path.isfile(grdfile):
+ command = '''
+cut -d' ' -f1,2,6 %(all_zfiles)s | \
+ surface -I%(az_resolution)s -G%(grdfile)s -R%(bounds)s -N1
+''' % vars()
+ os.system(command)
+
+ ## sampling the grdfile along the great circle
+ xyptfp = os.popen('grdtrack %(gcfile)s -G%(grdfile)s -Lg' % vars())
+
+ ## write the sampled results (azimuth, r, temperature) to a xect file
+ for line in xyptfp.readlines():
+ xypt = line.split()
+ out.write('%s\t%f\t%s\n' % (xypt[2], radius[layer], xypt[3]) )
+ xyptfp.close()
+
+out.close()
+
+
+### Plotting ###
+#print 'Plotting'
+psfile = '%s.xsect.%d.ps' % (modelname, step)
+mapwidth = 12.0
+proj = 'H0/%f' % mapwidth
+yshift = mapwidth * 1.5
+cptfile = 'zz.cpt'
+
+## colorbar length and location
+cbarh = mapwidth / 2
+cbxshift = mapwidth + .5
+cbyshift = cbarh / 2
+
+## plot the temperature field at mid-depth and a great circle
+grdfile = '%s.%d.z%03d.tgrd' % (modelname, step, int(nodez/2))
+command = '''
+makecpt -Cpolar -T0/1/.1 > %(cptfile)s
+
+grdimage %(grdfile)s -C%(cptfile)s -Bg360 -R%(bounds)s -J%(proj)s -X1.5 -Y%(yshift)f -P -K > %(psfile)s
+
+pscoast -R%(bounds)s -J%(proj)s -W -Dc -K -O >> %(psfile)s
+
+psxy %(gcfile)s -R%(bounds)s -J%(proj)s -W6. -O -K >> %(psfile)s
+
+gmtset ANOT_FONT_SIZE 9
+psscale -C%(cptfile)s -D%(cbxshift)f/%(cbyshift)f/%(cbarh)f/0.25 -K -O >> %(psfile)s
+''' % vars()
+os.system(command)
+
+
+## create a polar coordinate plot of the xsection
+##
+## TODO: there is always a gap on the left side of the annulus. How to fix it?
+grdfile2 = 'xsection.grd'
+bounds2 = '-180/180/%f/%f' % (radius[botlayer], radius[toplayer])
+r_resolution = (radius[toplayer] - radius[botlayer]) / 100
+resolution = '%f/%f' % (az_resolution, r_resolution)
+yshift = mapwidth * 1.2
+
+command = '''
+surface %(xsectfile)s -G%(grdfile2)s -I%(resolution)s -R%(bounds2)s
+
+grdimage %(grdfile2)s -C%(cptfile)s -JP%(mapwidth)fa -B30ns -R%(bounds2)s -X0.2 -Y-%(yshift)f -P -O >> %(psfile)s
+
+rm -f label.txt %(cptfile)s %(gcfile)s %(xsectfile)s %(grdfile2)s
+''' % vars()
+os.system(command)
+
+
+# version
+# $Id$
+
+# End of file
Property changes on: mc/3D/CitcomS/trunk/visual/plot_annulus.py
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mime-type
+ text/script
Name: svn:keywords
+ Id
Added: mc/3D/CitcomS/trunk/visual/plot_layer.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/plot_layer.py 2006-11-20 23:49:02 UTC (rev 5334)
+++ mc/3D/CitcomS/trunk/visual/plot_layer.py 2006-11-21 20:18:02 UTC (rev 5335)
@@ -0,0 +1,173 @@
+#!/usr/bin/env python
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+#<LicenseText>
+#
+# plot_layer.py by Eh Tan.
+# Copyright (C) 2002-2005, California Institute of Technology.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#</LicenseText>
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+"""Plotting CitcomS temperature at the specified layer
+
+usage: plot_layer.py modelname caps step layer
+ modelname: prefix of the combined capfile(s)
+ caps: 1 (regional) or 12 (full) caps
+ step: time step to plot
+ layer: which layer to plot (0 to nodez-1)
+
+input file:
+ modelname.capXX.step - CitcomS combined cap file
+ infile
+
+output file:
+ modelname.capXX.step.zYYY - extracted layer YYY from the cap file(s)
+ modelname.capXX.step.zYYY.ps - Postscript image of layer YYY
+
+"""
+
+
+def read_params(infile):
+ ## open input file and read header
+ nodex, nodey, nodez = open(infile).readline().split('x')
+ nodex = int(nodex)
+ nodey = int(nodey)
+ nodez = int(nodez)
+ #print nodex, nodey, nodez
+ return nodex, nodey, nodez
+
+
+def find_minmax(zfile, nodex, nodey):
+ fp = open(zfile)
+ x = range(nodex*nodey)
+ y = range(nodex*nodey)
+ n = 0
+ for line in fp.readlines():
+ x[n] = float(line.split()[0])
+ y[n] = float(line.split()[1])
+ n = n + 1
+ xmin = min(x)
+ xmax = max(x)
+ ymin = min(y)
+ ymax = max(y)
+
+ fp.close()
+
+ #print xmin, xmax, ymin, ymax
+ return xmin, xmax, ymin, ymax
+
+
+#######################################################################
+# Main
+#######################################################################
+import sys, os
+
+if len(sys.argv) != 5:
+ print __doc__
+ sys.exit(1)
+
+modelname = sys.argv[1]
+caps = int(sys.argv[2])
+step = int(sys.argv[3])
+layer = int(sys.argv[4])
+
+## read model parameters
+infile = '%s.cap%02d.%d' % (modelname, 0, step)
+nodex, nodey, nodez = read_params(infile)
+
+
+## slice the capfile(s), results saved in zfile
+import zslice
+all_zfiles = ''
+for cap in range(caps):
+ capfile = '%s.cap%02d.%d' % (modelname, cap, step)
+ zfile = zslice.zslicefile(capfile, layer)
+ all_zfiles = all_zfiles + ' ' + zfile
+ if not os.path.exists(zfile):
+ zslice.zslice(capfile, layer)
+
+
+#######################################################################
+# define GMT parameters
+#######################################################################
+
+## width of the plot
+mapwidth = 8.0
+
+if caps == 1:
+ ## find min/max of the coordinate
+ xmin, xmax, ymin, ymax = find_minmax(zfile, nodex, nodey)
+ bounds = '%f/%f/%f/%f' % (xmin,xmax,ymin,ymax)
+ proj = 'M%f' % mapwidth
+ resolution = '%f/%f' % ( (xmax-xmin)/nodex, (ymax-ymin)/nodey )
+else:
+ ## map centered at Pacific
+ bounds = '0/360/-90/90'
+ proj = 'H180/%d' % mapwidth
+ ## map centered at Greenwich
+ #bounds = '-180/180/-90/90'
+ #proj = 'H0/%d' % mapwidth
+ resolution = '0.5'
+
+cptfile = 'zz.cpt'
+grdfile = '%s.%d.z%03d.tgrd' % (modelname, step, layer)
+psfile = '%s.%d.z%03d.ps' % (modelname, step, layer)
+
+
+#print 'Plotting...'
+
+#######################################################################
+## call GMT commands to do the followings:
+## 1. cut the 1st, 2nd, 6th columns (lat, lon, temperature) of zfiles
+## 2. using GMT's "surface" to generate the grdfile
+## 3. generate color palette
+## 4. plot the grdfile
+## 5. plot the coastlines
+## 6. set the font size for the colorbar
+## 7. plot the colorbar
+## 8. remove the grdfile and cptfile
+#######################################################################
+command = '''
+cut -d' ' -f1,2,6 %(all_zfiles)s | \
+ surface -I%(resolution)s -G%(grdfile)s -R%(bounds)s
+
+makecpt -Cpolar -T0/1/.1 > %(cptfile)s
+
+grdimage %(grdfile)s -C%(cptfile)s -Bg90 -R%(bounds)s -J%(proj)s -X1.5 -Y5.0 -P -K > %(psfile)s
+
+pscoast -R%(bounds)s -J%(proj)s -Bg90 -W -Dc -K -O >> %(psfile)s
+
+gmtset ANOT_FONT_SIZE 9
+
+psscale -C%(cptfile)s -D8.5/2.25/4.0/0.25 -O >> %(psfile)s
+
+rm -f %(grdfile)s %(cptfile)s
+''' % vars()
+
+#print command
+os.system(command)
+
+#print 'Done'
+
+
+# version
+# $Id$
+
+# End of file
Property changes on: mc/3D/CitcomS/trunk/visual/plot_layer.py
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mime-type
+ text/script
Name: svn:keywords
+ Id
Modified: mc/3D/CitcomS/trunk/visual/zslice.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/zslice.py 2006-11-20 23:49:02 UTC (rev 5334)
+++ mc/3D/CitcomS/trunk/visual/zslice.py 2006-11-21 20:18:02 UTC (rev 5335)
@@ -46,64 +46,74 @@
prefix.z###
"""
-all = ('zslice')
+all = ('zslice', 'zslicefile')
+
+def zslicefile(prefix, layer):
+ return '%s.z%03d' % (prefix, layer)
+
+
def zslice(prefix, *ilayers):
from math import pi
r2d = 180.0 / pi
- # open input file and read header
- input = open(prefix)
- nodex, nodey, nodez = input.readline().split('x')
+ ## open input file and read header
+ infile = open(prefix)
+ nodex, nodey, nodez = infile.readline().split('x')
nodez = int(nodez)
#print nodez, ilayers
+ ## validate ilayers
layers = check_layers(ilayers, nodez)
nlayer = len(layers)
- # allocate arrays
+ ## allocate arrays
output = range(nlayer)
lines = range(nodez)
- # open output files
+ ## open output files
for i in range(nlayer):
- zfile = '%s.z%03d' % (prefix, layers[i])
+ zfile = zslicefile(prefix, layers[i])
output[i] = open(zfile, 'w')
- while 1:
- for j in range(nodez):
- lines[j] = input.readline()
+ try:
+ while 1:
+ ## read nodez lines
+ for j in range(nodez):
+ lines[j] = infile.readline()
- # file is empty or EOF?
- if not lines[nodez-1]:
- break
+ ## file is empty or EOF?
+ if not lines[nodez-1]:
+ break
- for i in range(nlayer):
- layer = layers[i]
+ for i in range(nlayer):
+ layer = layers[i]
- # spilt the first 3 columns only, the rest will be
- # output as-is
- data = lines[layer].split(' ', 3)
- lat = 90 - float(data[0])*r2d
- lon = float(data[1])*r2d
- output[i].write( '%f %f %s' % (lon, lat, data[3]) )
+ ## spilt the first 3 columns only, the rest will be
+ ## output as-is
+ data = lines[layer].split(' ', 3)
+ lat = 90 - float(data[0])*r2d
+ lon = float(data[1])*r2d
+ output[i].write( '%f %f %s' % (lon, lat, data[3]) )
- input.close()
+ infile.close()
- for i in range(nlayer):
- output[i].close()
+ finally:
+ for i in range(nlayer):
+ output[i].close()
+
return
def check_layers(layers, nodez):
if layers == ():
- # if empty, we will slice every layer
+ ## if empty, we will slice every layer
layers = range(nodez)
else:
- # otherwise, check bounds of layers
+ ## otherwise, check bounds of layers
for layer in layers:
if not (0<= layer < nodez):
raise ValueError, 'layer out of range (0-%d)' % (nodez - 1)
@@ -112,7 +122,7 @@
-# if run as a script
+## if run as a script
if __name__ == '__main__':
import sys
More information about the cig-commits
mailing list