PyEarthScience: Python examples for Earth Scientists

In [4]:

    

from   mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
import matplotlib.pyplot as plt
from   netCDF4 import *
import numpy as np


    

In [5]:

    

%matplotlib inline


    

In [6]:

    

#-- open netcdf file
nc = Dataset('/Users/k204045/NCL/general/data/new_data/rectilinear_grid_2D.nc', mode='r')

#-- read variable
u10 = nc.variables['u10'][0,:,:]
v10 = nc.variables['v10'][0,:,:]
lat = nc.variables['lat'][::-1]
lon = nc.variables['lon'][:]

u, lonsout = addcyclic(u10, lon)
v, lonsout = addcyclic(v10, lon)
print "lon[0]: ", lonsout[0], "lon[-1]: ", lonsout[-1]
print "lat[0]: ", lat[0], "lat[-1]: ", lat[-1]
print lonsout[:]
print lat[:]


    

    




lon[0]:  -180.0 lon[-1]:  180.0
lat[0]:  -88.572168514 lat[-1]:  88.572168514
[-180.    -178.125 -176.25  -174.375 -172.5   -170.625 -168.75  -166.875
 -165.    -163.125 -161.25  -159.375 -157.5   -155.625 -153.75  -151.875
 -150.    -148.125 -146.25  -144.375 -142.5   -140.625 -138.75  -136.875
 -135.    -133.125 -131.25  -129.375 -127.5   -125.625 -123.75  -121.875
 -120.    -118.125 -116.25  -114.375 -112.5   -110.625 -108.75  -106.875
 -105.    -103.125 -101.25   -99.375  -97.5    -95.625  -93.75   -91.875
  -90.     -88.125  -86.25   -84.375  -82.5    -80.625  -78.75   -76.875
  -75.     -73.125  -71.25   -69.375  -67.5    -65.625  -63.75   -61.875
  -60.     -58.125  -56.25   -54.375  -52.5    -50.625  -48.75   -46.875
  -45.     -43.125  -41.25   -39.375  -37.5    -35.625  -33.75   -31.875
  -30.     -28.125  -26.25   -24.375  -22.5    -20.625  -18.75   -16.875
  -15.     -13.125  -11.25    -9.375   -7.5     -5.625   -3.75    -1.875
    0.       1.875    3.75     5.625    7.5      9.375   11.25    13.125
   15.      16.875   18.75    20.625   22.5     24.375   26.25    28.125
   30.      31.875   33.75    35.625   37.5     39.375   41.25    43.125
   45.      46.875   48.75    50.625   52.5     54.375   56.25    58.125
   60.      61.875   63.75    65.625   67.5     69.375   71.25    73.125
   75.      76.875   78.75    80.625   82.5     84.375   86.25    88.125
   90.      91.875   93.75    95.625   97.5     99.375  101.25   103.125
  105.     106.875  108.75   110.625  112.5    114.375  116.25   118.125
  120.     121.875  123.75   125.625  127.5    129.375  131.25   133.125
  135.     136.875  138.75   140.625  142.5    144.375  146.25   148.125
  150.     151.875  153.75   155.625  157.5    159.375  161.25   163.125
  165.     166.875  168.75   170.625  172.5    174.375  176.25   178.125
  180.   ]
[-88.57216851 -86.72253095 -84.86197029 -82.99894164 -81.13497684
 -79.27055903 -77.40588808 -75.54106145 -73.67613231 -71.81113211
 -69.94608065 -68.08099099 -66.21587211 -64.35073041 -62.48557052
 -60.62039593 -58.75520927 -56.8900126  -55.02480754 -53.15959537
 -51.29437714 -49.4291537  -47.56392575 -45.69869388 -43.83345858
 -41.96822027 -40.1029793  -38.23773599 -36.37249059 -34.50724334
 -32.64199444 -30.77674406 -28.91149237 -27.0462395  -25.18098558
 -23.31573073 -21.45047504 -19.58521861 -17.71996153 -15.85470387
 -13.98944571 -12.12418712 -10.25892817  -8.39366891  -6.5284094
  -4.66314971  -2.79788988  -0.93262997   0.93262997   2.79788988
   4.66314971   6.5284094    8.39366891  10.25892817  12.12418712
  13.98944571  15.85470387  17.71996153  19.58521861  21.45047504
  23.31573073  25.18098558  27.0462395   28.91149237  30.77674406
  32.64199444  34.50724334  36.37249059  38.23773599  40.1029793
  41.96822027  43.83345858  45.69869388  47.56392575  49.4291537
  51.29437714  53.15959537  55.02480754  56.8900126   58.75520927
  60.62039593  62.48557052  64.35073041  66.21587211  68.08099099
  69.94608065  71.81113211  73.67613231  75.54106145  77.40588808
  79.27055903  81.13497684  82.99894164  84.86197029  86.72253095
  88.57216851]

In [7]:

    

#-- create figure and axes instances
fig = plt.figure(figsize=(8,8))
ax  = fig.add_axes([0.1,0.4,0.7,0.7])

#-- create map
map = Basemap(projection='cyl',llcrnrlat= -90.,urcrnrlat= 90.,\
              resolution='c',  llcrnrlon=-180.,urcrnrlon=180.)

#-- draw coastlines, state and country boundaries, edge of map
map.drawcoastlines()
map.drawstates()
map.drawcountries()

#-- create and draw meridians and parallels grid lines
map.drawparallels(np.arange( -90., 90.,30.),labels=[1,0,0,0],fontsize=10)
map.drawmeridians(np.arange(-180.,180.,30.),labels=[0,0,0,1],fontsize=10)

#-- convert latitude/longitude values to plot x/y values
x, y = map(*np.meshgrid(lon,lat))

#-- transform vector and coordinate data
veclon = u10.shape[1]/2                    #-- only every 2nd vector
veclat = u10.shape[0]/2                    #-- only every 2nd vector
uproj,vproj,xx,yy = map.transform_vector(u,v,lonsout,lat,veclon,veclat,returnxy=True,masked=True)

#-- create vector plot on map
vecplot = map.quiver(xx,yy,uproj,vproj,scale=600)
qk = plt.quiverkey(vecplot, 0.2, -0.2, 20, '20 m/s', labelpos='W')  #-- position and reference label

#-- add plot title
plt.title('Wind velocity')

plt.show()
#plt.savefig('plot_vector_matplotlib.png')