Surface air temperature as provided from IABP

   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"