1 #module general/gmt_tools.py
   2 # LK 14/08/2008
   3 
   4 from polar_projection import *
   5 from grids import def_regions
   6 import pipes,os,tempfile
   7 #Python 2.4
   8 #import PyNGL.Nio as Nio
   9 #Python 2.6
  10 #import Nio
  11 import scipy.io as io
  12 from numpy import array
  13 #run ../general/gmt_tools.py
  14 
  15 def cmd(cmdline):
  16     print cmdline
  17     os.system(cmdline)
  18     return
  19 
  20 def grid_par(region,cs,G):
  21     regions=def_regions()
  22     #print regions
  23     sgn,cds= regions[region][0],regions[region][1:5]
  24     x0,y0,x2,y2=cds[0],cds[1],cds[2],cds[3]
  25     xs,ys=(x0-x2)/cs,(y0-y2)/cs
  26     G['Rx']=' -R'+str(x2)+'/'+str(x0)+'/'+str(y2)+'/'+str(y0)+' '
  27     xwidth,ywidth=xs*cs,ys*cs
  28     lat2,lon2=mapxy(x0,y0,sgn)
  29     lat0,lon0=mapxy(x2,y2,sgn)
  30     map_x_size, map_y_size = 1.65, 2.45     #A4
  31     map_x_size, map_y_size = 1.5, 2.22     #A4 smaller
  32     
  33     map_ratio=map_x_size/map_y_size
  34     img_ratio=xwidth/ywidth
  35     width=xwidth
  36     paper=map_x_size
  37     #paper=map_y_size #Landscape
  38     latlon_corners=str(lon0)+'/'+str(lat0)+'/'+str(lon2)+'/'+str(lat2)
  39     scalell=long(width*10000.0/paper)
  40     if sgn==1:
  41         grid = '-Js-45/90/70/1:'
  42     if sgn==-1:
  43         grid = '-Js0/-90/-70/1:'
  44     latlon_par=' -R'+latlon_corners+'r '+grid+str(scalell)+' '
  45     G['Rll']=latlon_par
  46     scale=1/width*10.0*paper
  47     G['Jx']=' -Jx'+str(scale)+'c '
  48     G['I']=' -I'+str(cs)+' '
  49     G['S']=' -S'+str(cs*2)+' '
  50     G['N']=' -N4 '
  51     G['Bll']=' -Ba10mg10m/a5mg5mNEsW ' # fine
  52     G['Bll']=' -Ba10g10/a5g5neSW ' # coarse
  53     G['Bx']=' -Bnesw '
  54     G['coast']=' -Df -W0.5/0/255/0 '
  55 
  56     return G
  57 
  58 def anot_par(G): # Obsolete, now in grid_par
  59     G['Bll']=' -Ba10mg10m/a5mg5mNEsW ' # fine
  60     G['Bll']=' -Ba10g10/a5g5neSW ' # coarse
  61     G['Bx']=' -Bnesw '
  62     G['coast']=' -Dc -W0.5/0/255/0 '
  63     return G
  64 
  65 
  66 def makecpt(ctabfile,ctab,min,max,step,G,**kw):
  67    
  68     s=''
  69     if kw.has_key('I'):
  70         s+=' -I '
  71     cmd('makecpt -C'+ctab+s+' -D -Z -T'+str(min)+'/'+str(max)+'/'+str(step)+' >'+ctabfile)
  72     G['C']=' -C'+ctabfile
  73     return G
  74 
  75 
  76 
  77 def opts(G,opts):
  78     s=' '
  79     for o in opts:
  80         s=s+G[o]+' '
  81     return s
  82 
  83 def nearneighbor(xab,yab,I,filename,G):
  84     cmd='nearneighbor  '+G['Rx']+'  -F '+G['N']+G['I']+G['S']+' -bis -G'+filename
  85     print cmd
  86     pipe=pipes.Template()
  87     pipe.append(cmd,'-.')
  88     pipe.open(filename+'.pipe', 'w').write(array([xab.flatten(),yab.flatten(),I.flatten()]).swapaxes(0,1).tostring())
  89     return
  90 
  91 def xyz2grd(xab,yab,I,filename,G):
  92     cmd='xyz2grd  '+G['Rx']+' -V -F  '+G['I']+' -bis -G'+filename
  93     pipe=pipes.Template()
  94     pipe.append(cmd,'-.')
  95     pipe.open(filename+'.pipe', 'w').write(array([xab.flatten(),yab.flatten(),I.flatten()]).swapaxes(0,1).tostring())
  96     return
  97 
  98 def readnc(filename):
  99     g=io.netcdf_file(filename,'r')
 100     return array(g.variables['z'][:,:]).copy()
 101 
 102 #def readnc(filename):
 103 #    g=Nio.open_file(filename,'r')
 104 #    return g.variables['z'][:,:]
 105 
 106 def write_table(filename,x,y,z):
 107     l=open(filename,'w')
 108     for i in range(z.shape[0]):
 109         l.write(str(x[i])+' '+str(y[i])+' '+str(z[i])+'\n') # e.g. Lat lon TB
 110     l.close()
 111     return
 112 
 113 def write_table_flag(filename,x,y,z,f):
 114     l=open(filename,'w')
 115     for i in range(z.shape[0]):
 116         if f[i]==0:
 117             l.write(str(x[i])+' '+str(y[i])+' '+str(z[i])+'\n') # e.g. Lat lon TB
 118     l.close()
 119     return
 120 
 121 def x0y0x2y2(region):
 122     regions=def_regions()
 123     sgn,cds= regions[region][0],regions[region][1:5]
 124     return cds[0],cds[1],cds[2],cds[3]
 125 
 126 
 127 def xsys(region,cs):
 128     x0,y0,x2,y2=x0y0x2y2(region)
 129     xs,ys=int((x0-x2)/cs),int((y0-y2)/cs)
 130     return xs,ys
 131 
 132 def XYgrid(region,cs):
 133     x0,y0,x2,y2=x0y0x2y2(region)
 134     y=linspace(y2+cs/2,y0-cs/2,(y0-y2)/cs)
 135     x=linspace(x2+cs/2,x0-cs/2,(x0-x2)/cs)
 136     return meshgrid(x.astype(float32),y.astype(float32))
 137 
 138 def test_region(region):
 139     cs=25.0
 140     outfile='test.ps'
 141     print 'Drawing map for region '+region+' in '+outfile
 142     G=grid_par(region,cs,{})
 143 #    G=anot_par(G)
 144     cmd('gmtset ELLIPSOID WGS-84 OBLIQUE_ANOTATION 0 PAPER_MEDIA A4+ PAGE_ORIENTATION portrait PSIMAGE_FORMAT hex PLOT_DEGREE_FORMAT  ddd:mm:ss')
 145     G['Bx']=' -Ba1000g100/a1000g100NEsw '
 146     cmd('gmtset GRID_CROSS_SIZE 0.01i')
 147     cmd('psbasemap '+opts(G,['Rx','Bx','Jx'])+' -K  >'+outfile)
 148     cmd('gmtset GRID_CROSS_SIZE 0')
 149     cmd('pscoast '+opts(G,['Rll','Bll','coast'])+' -K -O >> '+outfile)
 150     titlestr=region+' '+opts(G,['Rx','Rll'])
 151     cmd('echo "0.00 -0.07 8 0.0 0 0 '+titlestr+'" | pstext -N -R0/1/0/1 -JX21/29 -O >> '+outfile)
 152     cmd('gv '+outfile)
 153 
 154 
 155 def print_regions():
 156     for r in def_regions().keys():
 157         print r,def_regions()[r]
 158 
 159 
 160 if __name__=='__main__':
 161     print_regions()