Surface air temperature as provided from IABP

   1 """
   2 import sys,os,os.path,glob,gzip
   3 from pylab import *
   4 import scipy.io as sio
   5 from mpl_toolkits.basemap import Basemap
   6 from gmt_tools import *
   7 import grids
   8 import Nio
   9     
  10 print "Reading coordinate lookup table"
  11 elatlon=loadtxt("elatlon",usecols=(1,2)).astype(float32)
  12 
  13 # Download data
  14 # !wget ftp://iabp.apl.washington.edu/pub/IABP/AirT/OI2004.fld.gz
  15 year="2004"
  16 
  17 filename="/scratch/local1/u242023/tmp/OI"+year+".fld.gz"
  18 
  19 
  20 print "Reading data from file ",filename," please wait ... (could take some minutes)"
  21 d=loadtxt(gzip.open(filename))
  22 print "Selecting indices"
  23 
  24 pixel_index_all=d[:,0].astype(int32)
  25 day_all=(d[:,1]/5).astype(int32)
  26 T_air_all=(d[:,2]/100.0).astype(float32)
  27 days=unique(day_all)
  28 
  29 
  30 for day_i in days:
  31 
  32     ind_day=(day_all==day_i) # select time index of one day
  33     pixel_index=pixel_index_all[ind_day] # get pixels indices of one day
  34     lat=elatlon[pixel_index,0] # Select lats and lons from lookup-table
  35     lon=elatlon[pixel_index,1]
  36     T_air=T_air_all[ind_day] # select data of one day
  37 
  38     # GRID data on polar sterographic grid
  39     cs=50.0 # 50 km
  40     region='Arctic'
  41     G=grid_par(region,cs,{})
  42     G['S']=' -S'+str(cs*10)+' '
  43     grdfile="/tmp/xyt_tmp.grd"
  44     outfile="out.ps"
  45     ctabfile="ctab"
  46 
  47     x,y=mapll(lat,lon,1)# polar stereographic projection
  48     
  49     # use GMT NN gridding
  50     nearneighbor(x.astype(float32),y.astype(float32),T_air.astype(float32),grdfile,G)
  51     fid=sio.netcdf_file(grdfile,'r')
  52     T_grid=fid.variables['z'][:,:].copy().astype(float32)
  53 
  54     # write to netcdf file
  55     filename_nc=year+'_'+str(day_i)+'.nc'
  56     F=Nio.open_file(filename_nc,'w') #write access
  57     F.create_dimension('my_y',240)
  58     F.create_dimension('my_x',240)
  59     F.create_dimension('my_time',None)
  60     F.create_variable('Air T','f',('my_y','my_x')) 
  61     F.variables['Air T'][:,:]=T_grid
  62     F.close()
  63 
  64 # Plot data in grdfile
  65 G=makecpt(ctabfile,'jet',-40,30,1,G)
  66 G=anot_par(G)
  67 G['Bll']=' -Ba30g30/a5g5nesw '
  68 G['coast']=' -Dc -W1/0/0/0  '
  69     
  70 cmd('grdimage '+grdfile+' '+opts(G,['Rx','Jx','C','Bx'])+' -Y3.5 -K > '+outfile)
  71 cmd('pscoast '+opts(G,['Rll','Bll','coast'])+'  -O  >> '+outfile)
  72 cmd('gv '+outfile+'&')
  73 

   1 """
   2 Example code to read in the surface air temperature as provided from IABP
   3 Lars Kaleschke, 4.2.2014
   4 
   5 """
   6 import sys,os,os.path,glob,gzip
   7 from pylab import *
   8 from mpl_toolkits.basemap import Basemap
   9 
  10 print "Reading coordinate lookup table"
  11 elatlon=loadtxt("elatlon",usecols=(1,2)).astype(float32)
  12 
  13 # Download data
  14 # !wget ftp://iabp.apl.washington.edu/pub/IABP/AirT/OI2004.fld.gz
  15 
  16 filename="OI2004.fld.gz"
  17 print "Reading data from file ",filename," please wait ... (could take some minutes)"
  18 d=loadtxt(gzip.open(filename))
  19 print "Selecting indices"
  20 
  21 pixel_index_all=d[:,0].astype(int32)
  22 day_all=(d[:,1]/5).astype(int32)
  23 T_air_all=(d[:,2]/100.0).astype(float32)
  24 days=unique(day_all)
  25 
  26 # Select one particular time step (could be used within a loop)
  27 day_i=5
  28 
  29 ind_day=(day_all==day_i) # select time index of one day
  30 pixel_index=pixel_index_all[ind_day] # get pixels indices of one day
  31 lat=elatlon[pixel_index,0] # Select lats and lons from lookup-table
  32 lon=elatlon[pixel_index,1]
  33 
  34 T_air=T_air_all[ind_day] # select data of one day
  35 
  36 print "Creating Basemap instance"
  37 m = Basemap(projection='ortho',lon_0=0.0,lat_0=90.0,resolution='l')
  38 
  39 print "Converting lat,lon -> x,y"
  40 x, y = m(lon,lat)
  41 
  42 print "Selecting unique coordinates to generate new grid"
  43 X=unique(x)
  44 Y=unique(y)[::-1]
  45 
  46 print "Gridding data to new grid"
  47 T_grid=griddata(x,y,T_air,X,Y)
  48 
  49 x_grid,y_grid=meshgrid(X,Y)
  50 ps=linspace(-60,40,101)
  51 
  52 print "Plotting data"
  53 m.contourf(x_grid,y_grid,T_grid,ps,cmap=cm.jet)
  54 m.drawcoastlines()
  55 m.drawmeridians(arange(0, 360, 30))
  56 m.drawparallels(arange(-90, 90, 30))
  57 colorbar()
  58 print "Show plot on screen"
  59 show()
  60 
  61 print "Finished"