Example 2

    
    from Numeric import *
    from pyclimate.diffoperators import *
    from Scientific.IO.NetCDF import *
    from pyclimate.ncstruct import *
    
    inc = NetCDFFile("NCARreanal200.nc")
    hgt = inc.variables["hgt"]
    lats = inc.variables["lat"]
    lons = inc.variables["lon"]
    hgtdata = hgt[0,:,:]
    grad = HGRADIENT(lats[:], lons[:])
    geosfactor = 2 * 7.29e-5 * sin(lats[:] * 3.14159 / 180.) / 9.81
    geosfactor = geosfactor[:,NewAxis] * ones(hgtdata.shape)
    geosfactor = geosfactor + equal(geosfactor, 0.0) 
    hg = grad.hgradient(hgtdata)
    wind = (-hg[1]/geosfactor, hg[0]/geosfactor) 
    
    dims = hgt.dimensions
    onc = nccopystruct("gwind200.nc", inc, dims, dims, dims)
    gw_u = onc.createVariable("gw_u", Float32, dims)
    gw_u.long_name = "geostrophic wind zonal component"
    gw_u.units = "m/s"
    gw_v = onc.createVariable("gw_v", Float32, dims)
    gw_v.long_name = "geostrophic wind meridional component"
    gw_v.units = "m/s"
    gw_u[0,:,:] = wind[0].astype(Float32)
    gw_v[0,:,:] = wind[1].astype(Float32)
    onc.close()
    

    grads -pb << EOF
    sdfopen NCARreanal200.nc
    set grads off
    d hgt
    d ugrd;vgrd
    run cbarn
    draw title 200 mb geopotencial height surface and real wind
    enable print tmp.gm
    print
    disable print
    quit
    EOF
    
    gxgif -i tmp.gm -o ejemplo2source.gif
    rm -f tmp.gm
    
    grads -pb << EOF
    sdfopen gwind200.nc
    sdfopen NCARreanal200.nc
    set grads off
    d hgt.2
    d gw_u;gw_v
    run cbarn
    draw title 200 mb geostrophic wind generated by ejemplo2.py
    enable print tmp.gm
    print
    disable print
    quit
    EOF
    
    gxgif -i tmp.gm -o ejemplo2.gif
    rm -f tmp.gm