= Surface air temperature as provided from IABP =
== Calculate Lat, Lon for gridded data ==

{{{#!python

import sys,os,os.path,glob,gzip
from pylab import *
import scipy.io as sio
from mpl_toolkits.basemap import Basemap
from gmt_tools import *
from polar_projection import *
import grids

fn='/scratch/local1/loew/climate_data_2014/project_B/data/data_yearly/2004_mean.nc'

fid=sio.netcdf_file(fn,'r')
T=fid.variables['Air_T'][0,:,:]

region='Arctic'
cs=50.0
X,Y=XYgrid(region,cs)
lat,lon=mapxy(X,Y,1) # calculate lat, lon
###############################


m = Basemap(projection='ortho',lon_0=0.0,lat_0=90.0,resolution='l')
x, y = m(lon,lat)
ps=linspace(-60,40,101)
m.contourf(x,y,T,ps,cmap=cm.jet)
m.drawcoastlines()
m.drawmeridians(arange(0, 360, 30))
m.drawparallels(arange(-90, 90, 30))
colorbar()
show()
}}}


== Re-gridding using GMT nearneighbor ==
{{attachment:out.png}}

{{{#!python
import sys,os,os.path,glob,gzip
from pylab import *
import scipy.io as sio
from mpl_toolkits.basemap import Basemap
from gmt_tools import *
import grids
import Nio
    
print "Reading coordinate lookup table"
elatlon=loadtxt("elatlon",usecols=(1,2)).astype(float32)

# Download data
# !wget ftp://iabp.apl.washington.edu/pub/IABP/AirT/OI2004.fld.gz
year="2004"

filename="/scratch/local1/u242023/tmp/OI"+year+".fld.gz"


print "Reading data from file ",filename," please wait ... (could take some minutes)"
d=loadtxt(gzip.open(filename))
print "Selecting indices"

pixel_index_all=d[:,0].astype(int32)
day_all=(d[:,1]/5).astype(int32)
T_air_all=(d[:,2]/100.0).astype(float32)
days=unique(day_all)


for day_i in days:

    ind_day=(day_all==day_i) # select time index of one day
    pixel_index=pixel_index_all[ind_day] # get pixels indices of one day
    lat=elatlon[pixel_index,0] # Select lats and lons from lookup-table
    lon=elatlon[pixel_index,1]
    T_air=T_air_all[ind_day] # select data of one day

    # GRID data on polar sterographic grid
    cs=50.0 # 50 km
    region='Arctic'
    G=grid_par(region,cs,{})
    G['S']=' -S'+str(cs*10)+' '
    grdfile="/tmp/xyt_tmp.grd"
    outfile="out.ps"
    ctabfile="ctab"

    x,y=mapll(lat,lon,1)# polar stereographic projection
    
    # use GMT NN gridding
    nearneighbor(x.astype(float32),y.astype(float32),T_air.astype(float32),grdfile,G)
    fid=sio.netcdf_file(grdfile,'r')
    T_grid=fid.variables['z'][:,:].copy().astype(float32)

    # write to netcdf file
    filename_nc=year+'_'+str(day_i)+'.nc'
    F=Nio.open_file(filename_nc,'w') #write access
    F.create_dimension('my_y',240)
    F.create_dimension('my_x',240)
    F.create_dimension('my_time',None)
    F.create_variable('Air T','f',('my_y','my_x')) 
    F.variables['Air T'][:,:]=T_grid
    F.close()

# Plot data in grdfile
G=makecpt(ctabfile,'jet',-40,30,1,G)
G=anot_par(G)
G['Bll']=' -Ba30g30/a5g5nesw '
G['coast']=' -Dc -W1/0/0/0  '
    
cmd('grdimage '+grdfile+' '+opts(G,['Rx','Jx','C','Bx'])+' -Y3.5 -K > '+outfile)
cmd('pscoast '+opts(G,['Rll','Bll','coast'])+'  -O  >> '+outfile)
cmd('gv '+outfile+'&')
}}}


== Re-gridding using Python griddata (slow) ==

{{attachment:iabp.png}}

{{{#!python
"""
Example code to read in the surface air temperature as provided from IABP
Lars Kaleschke, 4.2.2014

"""
import sys,os,os.path,glob,gzip
from pylab import *
from mpl_toolkits.basemap import Basemap

print "Reading coordinate lookup table"
elatlon=loadtxt("elatlon",usecols=(1,2)).astype(float32)

# Download data
# !wget ftp://iabp.apl.washington.edu/pub/IABP/AirT/OI2004.fld.gz

filename="OI2004.fld.gz"
print "Reading data from file ",filename," please wait ... (could take some minutes)"
d=loadtxt(gzip.open(filename))
print "Selecting indices"

pixel_index_all=d[:,0].astype(int32)
day_all=(d[:,1]/5).astype(int32)
T_air_all=(d[:,2]/100.0).astype(float32)
days=unique(day_all)

# Select one particular time step (could be used within a loop)
day_i=5

ind_day=(day_all==day_i) # select time index of one day
pixel_index=pixel_index_all[ind_day] # get pixels indices of one day
lat=elatlon[pixel_index,0] # Select lats and lons from lookup-table
lon=elatlon[pixel_index,1]

T_air=T_air_all[ind_day] # select data of one day

print "Creating Basemap instance"
m = Basemap(projection='ortho',lon_0=0.0,lat_0=90.0,resolution='l')

print "Converting lat,lon -> x,y"
x, y = m(lon,lat)

print "Selecting unique coordinates to generate new grid"
X=unique(x)
Y=unique(y)[::-1]

print "Gridding data to new grid"
T_grid=griddata(x,y,T_air,X,Y)

x_grid,y_grid=meshgrid(X,Y)
ps=linspace(-60,40,101)

print "Plotting data"
m.contourf(x_grid,y_grid,T_grid,ps,cmap=cm.jet)
m.drawcoastlines()
m.drawmeridians(arange(0, 360, 30))
m.drawparallels(arange(-90, 90, 30))
colorbar()
print "Show plot on screen"
show()

print "Finished"
}}}