Netcdf and Basemap

Main goals of this lecture:

  1. Open a netcdf file and make a simple plot
  2. Make a cross section plot
  3. Extract a specific point to a timeseries

note : The data are available at: https://drive.google.com/open?id=0B9LDS6UoXaUoUTJ2YVdFeWpXVlk

1. Open a netcdf file and make a simple plot


In [1]:
import numpy as np
from scipy.io import netcdf
from matplotlib import pyplot as plt
%matplotlib inline

In [4]:
f = netcdf.netcdf_file('./data/r2_reanalysis.nc', 'r')

The file can be iterated to show all of the variables


In [5]:
for var in f.variables:
    print var


time_bnds
prate
tmin
level
tmax
v
lon
u
t
time
lat
u10
v10
level_2
level_3

Now we can extract the variables that we want to use. Using the [:] indicates that a numpy array we be extract, otherwise an python object will be extracted


In [6]:
u = f.variables['u']
v = f.variables['v']
t = f.variables['t']

longitude = f.variables['lon'][:]
latitude = f.variables['lat'][:]
level = f.variables['level'][:]
time = f.variables['time']

print u.shape, longitude.shape, latitude.shape, level.shape, time.shape
print time.units


(182, 17, 73, 144) (144,) (73,) (17,) (182,)
hours since 1800-1-1 00:00:00

In [7]:
print latitude
print longitude
print level


[ 90.   87.5  85.   82.5  80.   77.5  75.   72.5  70.   67.5  65.   62.5
  60.   57.5  55.   52.5  50.   47.5  45.   42.5  40.   37.5  35.   32.5
  30.   27.5  25.   22.5  20.   17.5  15.   12.5  10.    7.5   5.    2.5
   0.   -2.5  -5.   -7.5 -10.  -12.5 -15.  -17.5 -20.  -22.5 -25.  -27.5
 -30.  -32.5 -35.  -37.5 -40.  -42.5 -45.  -47.5 -50.  -52.5 -55.  -57.5
 -60.  -62.5 -65.  -67.5 -70.  -72.5 -75.  -77.5 -80.  -82.5 -85.  -87.5
 -90. ]
[   0.     2.5    5.     7.5   10.    12.5   15.    17.5   20.    22.5
   25.    27.5   30.    32.5   35.    37.5   40.    42.5   45.    47.5
   50.    52.5   55.    57.5   60.    62.5   65.    67.5   70.    72.5
   75.    77.5   80.    82.5   85.    87.5   90.    92.5   95.    97.5
  100.   102.5  105.   107.5  110.   112.5  115.   117.5  120.   122.5
  125.   127.5  130.   132.5  135.   137.5  140.   142.5  145.   147.5
  150.   152.5  155.   157.5  160.   162.5  165.   167.5  170.   172.5
  175.   177.5  180.   182.5  185.   187.5  190.   192.5  195.   197.5
  200.   202.5  205.   207.5  210.   212.5  215.   217.5  220.   222.5
  225.   227.5  230.   232.5  235.   237.5  240.   242.5  245.   247.5
  250.   252.5  255.   257.5  260.   262.5  265.   267.5  270.   272.5
  275.   277.5  280.   282.5  285.   287.5  290.   292.5  295.   297.5
  300.   302.5  305.   307.5  310.   312.5  315.   317.5  320.   322.5
  325.   327.5  330.   332.5  335.   337.5  340.   342.5  345.   347.5
  350.   352.5  355.   357.5]
[ 1000.   925.   850.   700.   600.   500.   400.   300.   250.   200.
   150.   100.    70.    50.    30.    20.    10.]

netcdf uses a type of compression for the variables that requires to be multiplied by an scale factor (scale_factor) and added an offset value (add_offset). For example, the scale_factor for the zonal wind U is 0.01 and it's offset value is 187.65. If we don't apply this changes the first element of the u array is u[0,0,0,0] = -18573 and if we apply those changes it becomes u[0,0,0,0] = 1.92. The module netCDF4 (http://unidata.github.io/netcdf4-python/) automagically handles those attributes saving up some time.


In [8]:
def unpack(var):
    return var[:] * var.scale_factor + var.add_offset

In [9]:
u = unpack(u)
v = unpack(v)
t = unpack(t)

The np.meshgrid returns a coordinate matrice from coordinate vector. In this example it will associate each of the longitudes in longitude vector to each of the latitudes in latitude and vice versa. Thereby, an 2D array will be created


In [10]:
lons, lats = np.meshgrid(longitude, latitude)
print lons.shape, lats.shape


(73, 144) (73, 144)

To plot that data we can use Basemap http://matplotlib.org/basemap/users/index.html

" The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python (...) Basemap is geared toward the needs of earth scientists, particular oceanographers and meteorologists (...) "


In [11]:
from mpl_toolkits.basemap import Basemap

A simple image can be plotted using the min and max values in the lats and lons variables. We can create a map object m defining the lower corner and uper corner longitude and latitude, llcrnrlon(lat) and urcrnrlon(lat). Then we can use the matplotlib methods to make the temperature and wind plot.


In [12]:
fig = plt.figure(figsize=(16,35))
m = Basemap(llcrnrlon=lons.min(), llcrnrlat=lats.min(), urcrnrlon=lons.max(), urcrnrlat=lats.max())

m.drawcoastlines()
m.drawcountries()

skip = 2

cs = m.contourf(lons, lats, t[0, 0, :, :])
qv = m.quiver(lons[::skip, ::skip], lats[::skip, ::skip], u[0, 0, ::skip, ::skip], v[0, 0, ::skip, ::skip])


There are many different projections to choose from in Basemap:

  • Azimuthal Equidistant Projection
  • Gnomonic Projection
  • Orthographic Projection
  • Geostationary Projection
  • Near-Sided Perspective Projection
  • Mollweide Projection
  • Hammer Projection
  • Robinson Projection
  • Eckert IV Projection
  • Kavrayskiy VII Projection
  • McBryde-Thomas Flat Polar Quartic
  • Sinusoidal Projection
  • Equidistant Cylindrical Projection
  • Cassini Projection
  • Mercator Projection
  • Transverse Mercator Projection
  • Oblique Mercator Projection
  • Polyconic Projection
  • Miller Cylindrical Projection
  • Gall Stereographic Projection
  • Cylindrial Equal-Area Projection
  • Lambert Conformal Projection
  • Lambert Azimuthal Equal Area Projection
  • Stereographic Projection
  • Equidistant Conic Projection
  • Albers Equal Area Projection
  • Polar Stereographic Projection
  • Polar Lambert Azimuthal Projection
  • Polar Azimuthal Equidistant Projection
  • van der Grinten Projection

In [13]:
fig = plt.figure(figsize=(20,20))

ax1 = fig.add_subplot(121)
ax1.set_title('Orthographic Projection')

m = Basemap(projection='ortho',lat_0=-15,lon_0=300, ax=ax1)

m.drawcoastlines()
m.drawcountries()

skip = 2

x, y = m(lons, lats)

cs = m.contourf(x, y, t[0, 0, :, :])
qv = m.quiver(x[::skip, ::skip], y[::skip, ::skip], u[0, 0, ::skip, ::skip], v[0, 0, ::skip, ::skip])

###############################################################################################

ax2 = fig.add_subplot(122)
ax2.set_title('Polar Azimuthal Equidistant Projection')

m = Basemap(projection='spaeqd',boundinglat=-10,lon_0=270, ax=ax2)

m.drawcoastlines()
m.drawcountries()

skip = 2

x, y = m(lons, lats)

cs = m.contourf(x, y, t[0, 0, :, :])
qv = m.quiver(x[::skip, ::skip], y[::skip, ::skip], u[0, 0, ::skip, ::skip], v[0, 0, ::skip, ::skip])


Question: What would happen if i wanted to make a plot between -180 and 180 instead of 0 to 360?

Answer: Since my longitude and my variables are ordered from 0 to 360, only a part of the data would be plotted (the part between 0 and 180). The part between -180 and 0 which equals 180 to 360, in the data between 0 and 360 would not be plotted


In [14]:
fig = plt.figure(figsize=(16,35))
m = Basemap(llcrnrlon=-180, llcrnrlat=lats.min(), urcrnrlon=180, urcrnrlat=lats.max())

m.drawcoastlines()
m.drawcountries()

skip = 2

cs = m.contourf(lons, lats, t[0, 0, :, :])
qv = m.quiver(lons[::skip, ::skip], lats[::skip, ::skip], u[0, 0, ::skip, ::skip], v[0, 0, ::skip, ::skip])


To reorder the variables and the longitude we can create a flip_grid function, that takes the longitude and variable to be reordered as arguments. This function will define a filter fltr that will be True for every longitude higher than 180 and False for values below. The filter will then be used to concatenate the values between -180 to -2.5 with the values between 0 to 177.5 for my longitude vector and my variable


In [15]:
def flip_grid(var, lons):
    fltr = lons >= 180
    # fltr =  [False False False ... True  True  True]
    newlons = np.concatenate(((lons - 360)[fltr], lons[~fltr]), axis=-1)
    # newlons = [-180 -177.5 -175 ... -5 -2.5 ] concatenated with [0 2.5 5 ... 175 177.5]
    # newlons = [-180 -177.5 -175 ... 175 177.5 180]
    if var.ndim == 2:
        newvar = np.concatenate((var[:, fltr], var[:, ~fltr]), axis=-1)
    elif var.ndim == 3:
        newvar = np.concatenate((var[:, :, fltr], var[:, :, ~fltr]), axis=-1)
    elif var.ndim == 4:
        newvar = np.concatenate((var[:, :, :, fltr], var[:, :, :, ~fltr]), axis=-1)        
        
    return newvar, newlons

The flip_grid is then used to reorder u, v and t variables


In [16]:
u, newlon = flip_grid(u, longitude)
v, newlon = flip_grid(v, longitude)
t, newlon = flip_grid(t, longitude)

print newlon


[-180.  -177.5 -175.  -172.5 -170.  -167.5 -165.  -162.5 -160.  -157.5
 -155.  -152.5 -150.  -147.5 -145.  -142.5 -140.  -137.5 -135.  -132.5
 -130.  -127.5 -125.  -122.5 -120.  -117.5 -115.  -112.5 -110.  -107.5
 -105.  -102.5 -100.   -97.5  -95.   -92.5  -90.   -87.5  -85.   -82.5
  -80.   -77.5  -75.   -72.5  -70.   -67.5  -65.   -62.5  -60.   -57.5
  -55.   -52.5  -50.   -47.5  -45.   -42.5  -40.   -37.5  -35.   -32.5
  -30.   -27.5  -25.   -22.5  -20.   -17.5  -15.   -12.5  -10.    -7.5
   -5.    -2.5    0.     2.5    5.     7.5   10.    12.5   15.    17.5
   20.    22.5   25.    27.5   30.    32.5   35.    37.5   40.    42.5
   45.    47.5   50.    52.5   55.    57.5   60.    62.5   65.    67.5
   70.    72.5   75.    77.5   80.    82.5   85.    87.5   90.    92.5
   95.    97.5  100.   102.5  105.   107.5  110.   112.5  115.   117.5
  120.   122.5  125.   127.5  130.   132.5  135.   137.5  140.   142.5
  145.   147.5  150.   152.5  155.   157.5  160.   162.5  165.   167.5
  170.   172.5  175.   177.5]

In [17]:
lons, lats = np.meshgrid(newlon, latitude)

In [18]:
fig = plt.figure(figsize=(16,35))
m = Basemap(llcrnrlon=lons.min(), llcrnrlat=lats.min(), urcrnrlon=lons.max(), urcrnrlat=lats.max())

m.drawcoastlines()
m.drawcountries()

skip = 2

cs = m.contourf(lons, lats, t[0, 0, :, :])
qv = m.quiver(lons[::skip, ::skip], lats[::skip, ::skip], u[0, 0, ::skip, ::skip], v[0, 0, ::skip, ::skip])


2. Make a cross section plot

Next, a function find_nearest is created to find the nearest grid point for a given location, based on the distance between points


In [19]:
def find_nearest(x, y, gridx, gridy):

    distance = (gridx - x)**2 + (gridy - y)**2
    idx = np.where(distance == distance.min())
    # idx = (array([45]), array([0]))
    
    return [idx[0][0], idx[1][0]]

Using, for example, the latitude and longitude of São Paulo (lat = -23.650000 and lon = -46.616667)


In [20]:
idx = find_nearest(-46.616667, -23.650000, lons, lats)

print lons[idx[0], idx[1]], lats[idx[0], idx[1]]


-47.5 -22.5

Selecting the index idx[1], which means the longitude we can make a plot for u[0, :, :, idx[1]] which is the zonal wind for all levels, all latitudes in the nearest longitude to -46.616667


In [21]:
fig = plt.figure(figsize=(20,8))
ax = fig.add_subplot(111)

newlats, levs = np.meshgrid(lats[:, idx[1]], level)

clevs = np.linspace(-50, 50, 8)
cs = ax.contourf(newlats, levs, u[0, :, :, idx[1]], clevs, cmap='PiYG', extend='both')
cbar = plt.colorbar(cs)
ax.invert_yaxis()


3. Extract a specific point to a timeseries

The idea is to extract a time series for the variables near the surface. This means a dataframe of this kind:

2016-01-01 300.700012
2016-01-02 299.500000
2016-01-03 299.279999
2016-01-04 298.349976
...

First we need to take a look in the time variable...


In [22]:
print time[:]


[ 1893408.  1893432.  1893456.  1893480.  1893504.  1893528.  1893552.
  1893576.  1893600.  1893624.  1893648.  1893672.  1893696.  1893720.
  1893744.  1893768.  1893792.  1893816.  1893840.  1893864.  1893888.
  1893912.  1893936.  1893960.  1893984.  1894008.  1894032.  1894056.
  1894080.  1894104.  1894128.  1894152.  1894176.  1894200.  1894224.
  1894248.  1894272.  1894296.  1894320.  1894344.  1894368.  1894392.
  1894416.  1894440.  1894464.  1894488.  1894512.  1894536.  1894560.
  1894584.  1894608.  1894632.  1894656.  1894680.  1894704.  1894728.
  1894752.  1894776.  1894800.  1894824.  1894848.  1894872.  1894896.
  1894920.  1894944.  1894968.  1894992.  1895016.  1895040.  1895064.
  1895088.  1895112.  1895136.  1895160.  1895184.  1895208.  1895232.
  1895256.  1895280.  1895304.  1895328.  1895352.  1895376.  1895400.
  1895424.  1895448.  1895472.  1895496.  1895520.  1895544.  1895568.
  1895592.  1895616.  1895640.  1895664.  1895688.  1895712.  1895736.
  1895760.  1895784.  1895808.  1895832.  1895856.  1895880.  1895904.
  1895928.  1895952.  1895976.  1896000.  1896024.  1896048.  1896072.
  1896096.  1896120.  1896144.  1896168.  1896192.  1896216.  1896240.
  1896264.  1896288.  1896312.  1896336.  1896360.  1896384.  1896408.
  1896432.  1896456.  1896480.  1896504.  1896528.  1896552.  1896576.
  1896600.  1896624.  1896648.  1896672.  1896696.  1896720.  1896744.
  1896768.  1896792.  1896816.  1896840.  1896864.  1896888.  1896912.
  1896936.  1896960.  1896984.  1897008.  1897032.  1897056.  1897080.
  1897104.  1897128.  1897152.  1897176.  1897200.  1897224.  1897248.
  1897272.  1897296.  1897320.  1897344.  1897368.  1897392.  1897416.
  1897440.  1897464.  1897488.  1897512.  1897536.  1897560.  1897584.
  1897608.  1897632.  1897656.  1897680.  1897704.  1897728.  1897752.]

Dates in netcdf files are generally computed as hours or days since an especific date. The ncdump -h for file r2.nc:

    double time(time) ;
            time:standard_name = "time" ;
            time:long_name = "Time" ;
            time:bounds = "time_bnds" ;
            time:units = "hours since 1800-1-1 00:00:00" ;
            time:calendar = "standard" ;
            time:axis = "T" ;

In [23]:
print time.units


hours since 1800-1-1 00:00:00

We may now import datetime and timedelta to make a list with datetime objects instead of those float. In the netCDF4 this convertion is as easy as dates = netCDF4.num2date(time[:], time.units)


In [24]:
from datetime import datetime, timedelta

In [25]:
def get_dates(time):
    splited_date = time.units.split()
    # splited_date = ['hours', 'since', '1800-1-1', '00:00:00']
    start_date_string = ' '.join(splited_date[-2:])
    # start_date_string = 1800-1-1 00:00:00
    # convert string into datetime object
    start_date = datetime.strptime(start_date_string, '%Y-%m-%d %H:%M:%S')
    
    dates = [start_date + timedelta(hours=i) for i in time[:]]
    return dates

In [26]:
dates = get_dates(time)
print dates


[datetime.datetime(2016, 1, 1, 0, 0), datetime.datetime(2016, 1, 2, 0, 0), datetime.datetime(2016, 1, 3, 0, 0), datetime.datetime(2016, 1, 4, 0, 0), datetime.datetime(2016, 1, 5, 0, 0), datetime.datetime(2016, 1, 6, 0, 0), datetime.datetime(2016, 1, 7, 0, 0), datetime.datetime(2016, 1, 8, 0, 0), datetime.datetime(2016, 1, 9, 0, 0), datetime.datetime(2016, 1, 10, 0, 0), datetime.datetime(2016, 1, 11, 0, 0), datetime.datetime(2016, 1, 12, 0, 0), datetime.datetime(2016, 1, 13, 0, 0), datetime.datetime(2016, 1, 14, 0, 0), datetime.datetime(2016, 1, 15, 0, 0), datetime.datetime(2016, 1, 16, 0, 0), datetime.datetime(2016, 1, 17, 0, 0), datetime.datetime(2016, 1, 18, 0, 0), datetime.datetime(2016, 1, 19, 0, 0), datetime.datetime(2016, 1, 20, 0, 0), datetime.datetime(2016, 1, 21, 0, 0), datetime.datetime(2016, 1, 22, 0, 0), datetime.datetime(2016, 1, 23, 0, 0), datetime.datetime(2016, 1, 24, 0, 0), datetime.datetime(2016, 1, 25, 0, 0), datetime.datetime(2016, 1, 26, 0, 0), datetime.datetime(2016, 1, 27, 0, 0), datetime.datetime(2016, 1, 28, 0, 0), datetime.datetime(2016, 1, 29, 0, 0), datetime.datetime(2016, 1, 30, 0, 0), datetime.datetime(2016, 1, 31, 0, 0), datetime.datetime(2016, 2, 1, 0, 0), datetime.datetime(2016, 2, 2, 0, 0), datetime.datetime(2016, 2, 3, 0, 0), datetime.datetime(2016, 2, 4, 0, 0), datetime.datetime(2016, 2, 5, 0, 0), datetime.datetime(2016, 2, 6, 0, 0), datetime.datetime(2016, 2, 7, 0, 0), datetime.datetime(2016, 2, 8, 0, 0), datetime.datetime(2016, 2, 9, 0, 0), datetime.datetime(2016, 2, 10, 0, 0), datetime.datetime(2016, 2, 11, 0, 0), datetime.datetime(2016, 2, 12, 0, 0), datetime.datetime(2016, 2, 13, 0, 0), datetime.datetime(2016, 2, 14, 0, 0), datetime.datetime(2016, 2, 15, 0, 0), datetime.datetime(2016, 2, 16, 0, 0), datetime.datetime(2016, 2, 17, 0, 0), datetime.datetime(2016, 2, 18, 0, 0), datetime.datetime(2016, 2, 19, 0, 0), datetime.datetime(2016, 2, 20, 0, 0), datetime.datetime(2016, 2, 21, 0, 0), datetime.datetime(2016, 2, 22, 0, 0), datetime.datetime(2016, 2, 23, 0, 0), datetime.datetime(2016, 2, 24, 0, 0), datetime.datetime(2016, 2, 25, 0, 0), datetime.datetime(2016, 2, 26, 0, 0), datetime.datetime(2016, 2, 27, 0, 0), datetime.datetime(2016, 2, 28, 0, 0), datetime.datetime(2016, 2, 29, 0, 0), datetime.datetime(2016, 3, 1, 0, 0), datetime.datetime(2016, 3, 2, 0, 0), datetime.datetime(2016, 3, 3, 0, 0), datetime.datetime(2016, 3, 4, 0, 0), datetime.datetime(2016, 3, 5, 0, 0), datetime.datetime(2016, 3, 6, 0, 0), datetime.datetime(2016, 3, 7, 0, 0), datetime.datetime(2016, 3, 8, 0, 0), datetime.datetime(2016, 3, 9, 0, 0), datetime.datetime(2016, 3, 10, 0, 0), datetime.datetime(2016, 3, 11, 0, 0), datetime.datetime(2016, 3, 12, 0, 0), datetime.datetime(2016, 3, 13, 0, 0), datetime.datetime(2016, 3, 14, 0, 0), datetime.datetime(2016, 3, 15, 0, 0), datetime.datetime(2016, 3, 16, 0, 0), datetime.datetime(2016, 3, 17, 0, 0), datetime.datetime(2016, 3, 18, 0, 0), datetime.datetime(2016, 3, 19, 0, 0), datetime.datetime(2016, 3, 20, 0, 0), datetime.datetime(2016, 3, 21, 0, 0), datetime.datetime(2016, 3, 22, 0, 0), datetime.datetime(2016, 3, 23, 0, 0), datetime.datetime(2016, 3, 24, 0, 0), datetime.datetime(2016, 3, 25, 0, 0), datetime.datetime(2016, 3, 26, 0, 0), datetime.datetime(2016, 3, 27, 0, 0), datetime.datetime(2016, 3, 28, 0, 0), datetime.datetime(2016, 3, 29, 0, 0), datetime.datetime(2016, 3, 30, 0, 0), datetime.datetime(2016, 3, 31, 0, 0), datetime.datetime(2016, 4, 1, 0, 0), datetime.datetime(2016, 4, 2, 0, 0), datetime.datetime(2016, 4, 3, 0, 0), datetime.datetime(2016, 4, 4, 0, 0), datetime.datetime(2016, 4, 5, 0, 0), datetime.datetime(2016, 4, 6, 0, 0), datetime.datetime(2016, 4, 7, 0, 0), datetime.datetime(2016, 4, 8, 0, 0), datetime.datetime(2016, 4, 9, 0, 0), datetime.datetime(2016, 4, 10, 0, 0), datetime.datetime(2016, 4, 11, 0, 0), datetime.datetime(2016, 4, 12, 0, 0), datetime.datetime(2016, 4, 13, 0, 0), datetime.datetime(2016, 4, 14, 0, 0), datetime.datetime(2016, 4, 15, 0, 0), datetime.datetime(2016, 4, 16, 0, 0), datetime.datetime(2016, 4, 17, 0, 0), datetime.datetime(2016, 4, 18, 0, 0), datetime.datetime(2016, 4, 19, 0, 0), datetime.datetime(2016, 4, 20, 0, 0), datetime.datetime(2016, 4, 21, 0, 0), datetime.datetime(2016, 4, 22, 0, 0), datetime.datetime(2016, 4, 23, 0, 0), datetime.datetime(2016, 4, 24, 0, 0), datetime.datetime(2016, 4, 25, 0, 0), datetime.datetime(2016, 4, 26, 0, 0), datetime.datetime(2016, 4, 27, 0, 0), datetime.datetime(2016, 4, 28, 0, 0), datetime.datetime(2016, 4, 29, 0, 0), datetime.datetime(2016, 4, 30, 0, 0), datetime.datetime(2016, 5, 1, 0, 0), datetime.datetime(2016, 5, 2, 0, 0), datetime.datetime(2016, 5, 3, 0, 0), datetime.datetime(2016, 5, 4, 0, 0), datetime.datetime(2016, 5, 5, 0, 0), datetime.datetime(2016, 5, 6, 0, 0), datetime.datetime(2016, 5, 7, 0, 0), datetime.datetime(2016, 5, 8, 0, 0), datetime.datetime(2016, 5, 9, 0, 0), datetime.datetime(2016, 5, 10, 0, 0), datetime.datetime(2016, 5, 11, 0, 0), datetime.datetime(2016, 5, 12, 0, 0), datetime.datetime(2016, 5, 13, 0, 0), datetime.datetime(2016, 5, 14, 0, 0), datetime.datetime(2016, 5, 15, 0, 0), datetime.datetime(2016, 5, 16, 0, 0), datetime.datetime(2016, 5, 17, 0, 0), datetime.datetime(2016, 5, 18, 0, 0), datetime.datetime(2016, 5, 19, 0, 0), datetime.datetime(2016, 5, 20, 0, 0), datetime.datetime(2016, 5, 21, 0, 0), datetime.datetime(2016, 5, 22, 0, 0), datetime.datetime(2016, 5, 23, 0, 0), datetime.datetime(2016, 5, 24, 0, 0), datetime.datetime(2016, 5, 25, 0, 0), datetime.datetime(2016, 5, 26, 0, 0), datetime.datetime(2016, 5, 27, 0, 0), datetime.datetime(2016, 5, 28, 0, 0), datetime.datetime(2016, 5, 29, 0, 0), datetime.datetime(2016, 5, 30, 0, 0), datetime.datetime(2016, 5, 31, 0, 0), datetime.datetime(2016, 6, 1, 0, 0), datetime.datetime(2016, 6, 2, 0, 0), datetime.datetime(2016, 6, 3, 0, 0), datetime.datetime(2016, 6, 4, 0, 0), datetime.datetime(2016, 6, 5, 0, 0), datetime.datetime(2016, 6, 6, 0, 0), datetime.datetime(2016, 6, 7, 0, 0), datetime.datetime(2016, 6, 8, 0, 0), datetime.datetime(2016, 6, 9, 0, 0), datetime.datetime(2016, 6, 10, 0, 0), datetime.datetime(2016, 6, 11, 0, 0), datetime.datetime(2016, 6, 12, 0, 0), datetime.datetime(2016, 6, 13, 0, 0), datetime.datetime(2016, 6, 14, 0, 0), datetime.datetime(2016, 6, 15, 0, 0), datetime.datetime(2016, 6, 16, 0, 0), datetime.datetime(2016, 6, 17, 0, 0), datetime.datetime(2016, 6, 18, 0, 0), datetime.datetime(2016, 6, 19, 0, 0), datetime.datetime(2016, 6, 20, 0, 0), datetime.datetime(2016, 6, 21, 0, 0), datetime.datetime(2016, 6, 22, 0, 0), datetime.datetime(2016, 6, 23, 0, 0), datetime.datetime(2016, 6, 24, 0, 0), datetime.datetime(2016, 6, 25, 0, 0), datetime.datetime(2016, 6, 26, 0, 0), datetime.datetime(2016, 6, 27, 0, 0), datetime.datetime(2016, 6, 28, 0, 0), datetime.datetime(2016, 6, 29, 0, 0), datetime.datetime(2016, 6, 30, 0, 0)]

We read now the variables near the surface (t2, u10, v10, prate)


In [27]:
u10 = f.variables['u10']
v10 = f.variables['v10']
prate = f.variables['prate']
tmax = f.variables['tmax']
tmin = f.variables['tmin']

print u10.shape, prate.units


(182, 1, 73, 144) Kg/m^2/s

The values of the variables are then unpacked and the grid is fliped from 0 to 360 to -180 to 180


In [28]:
u10 = unpack(u10)
v10 = unpack(v10)
prate = unpack(prate)
tmax = unpack(tmax)
tmin = unpack(tmin)

u10, newlon = flip_grid(u10, longitude)
v10, newlon = flip_grid(v10, longitude)
prate, newlon = flip_grid(prate, longitude)
tmax, newlon = flip_grid(tmax, longitude)
tmin, newlon = flip_grid(tmin, longitude)

In [29]:
lons, lats = np.meshgrid(newlon, latitude)

We find then grid point of the latitude and longitude that we want


In [30]:
idx = find_nearest(-46.248889, -22.880833, lons, lats)

Import pandas to create a dataframe with dates as index and then save it to a csv file


In [31]:
import pandas as pd

tx2C = tmax - 273.15 # K to oC
tm2C = tmin - 273.15 # K to oC
acum = prate * 3600. * 24. # mm/s mm/day

data = {'u10': u10[:, 0, idx[0], idx[1]], 'v10': v10[:, 0, idx[0], idx[1]], 'prec': acum[:, idx[0], idx[1]], \
        'tmax': tx2C[:, 0, idx[0], idx[1]], 'tmin': tm2C[:, 0, idx[0], idx[1]]}

df = pd.DataFrame(data, index=dates)

print df


                 prec       tmax       tmin       u10       v10
2016-01-01  28.287354  27.130005  22.820007  1.259995 -0.330002
2016-01-02  19.500477  24.510010  22.410004 -2.180008  1.959991
2016-01-03  30.196798  23.690002  22.010010 -5.009995  3.709991
2016-01-04   2.488319  23.730011  20.540009 -6.699997  2.449997
2016-01-05   0.552944  24.720001  20.050018 -6.630005 -0.199997
2016-01-06   0.008630  26.760010  19.910004 -5.449997 -0.350006
2016-01-07   0.008630  28.420013  20.510010 -3.830002 -0.570007
2016-01-08   0.008630  28.949982  22.279999 -3.580002 -0.880005
2016-01-09  32.054382  27.160004  22.690002 -3.610001 -1.300003
2016-01-10   4.112631  25.760010  22.380005 -1.430008 -2.910004
2016-01-11   1.615662  26.800018  22.519989  0.800003 -4.610001
2016-01-12  35.959671  23.720001  22.279999 -0.479996 -4.009995
2016-01-13  59.382717  24.899994  22.490021 -1.210007 -0.509995
2016-01-14  38.119667  25.050018  22.110016 -2.419998  4.000000
2016-01-15  30.784304  24.670013  22.050018 -1.430008  1.369995
2016-01-16  19.846062  22.910004  20.680023 -2.330002  7.039993
2016-01-17   1.848954  22.010010  18.730011 -3.570007  5.819992
2016-01-18   9.287981  21.880005  19.889984 -5.130005  5.599991
2016-01-19  14.195518  21.639984  19.880005 -5.070007  4.529999
2016-01-20  20.675507  22.190002  19.630005 -4.770004  4.250000
2016-01-21   6.488639  22.490021  19.870026 -1.919998  6.110001
2016-01-22   0.224622  23.980011  18.970001 -2.300003  4.929993
2016-01-23   0.086381  24.459991  18.850006 -1.570007  3.119995
2016-01-24   0.008630  26.130005  18.290009 -2.470001 -0.380005
2016-01-25   7.819189  26.470001  21.070007 -1.630005 -1.520004
2016-01-26   7.421746  28.010010  21.760010  0.709991 -0.830002
2016-01-27  11.188795  28.399994  22.630005 -0.949997  1.300003
2016-01-28   3.784309  24.959991  22.139984 -3.180008  1.559998
2016-01-29  18.187187  24.800018  22.100006 -2.500000 -1.970001
2016-01-30   1.105908  27.709991  21.449982 -1.509995 -3.229996
...               ...        ...        ...       ...       ...
2016-06-01   0.034540  23.380005  19.250000 -0.770004 -0.490005
2016-06-02   1.054068  23.360016  20.290009  1.569992 -2.500000
2016-06-03   1.719343  22.490021  19.579987 -0.270004 -1.850006
2016-06-04   2.574719  21.699982  19.209991 -0.300003  0.100006
2016-06-05   4.553265  22.269989  20.329987 -0.779999 -1.279999
2016-06-06   2.626540  23.570007  20.300018  0.619995 -1.080002
2016-06-07   6.220787  22.300018  17.990021  0.979996  4.160004
2016-06-08   1.183679  19.010010  16.230011  0.160004  3.470001
2016-06-09   0.025910  18.570007  14.180023 -0.559998  3.300003
2016-06-10   0.008630  20.260010  13.389984 -0.910004  2.029999
2016-06-11   0.371513  18.209991  12.649994  2.619995  2.929993
2016-06-12   0.008630  16.880005   8.889984  1.800003  2.819992
2016-06-13   0.008630  16.320007   8.790009 -1.360001  3.130005
2016-06-14   0.008630  18.990021  12.740021 -4.250000 -0.279999
2016-06-15   0.008630  21.079987  14.940002 -2.330002 -1.320007
2016-06-16   0.008630  22.940002  13.670013 -1.380005 -1.649994
2016-06-17   0.008630  22.139984  14.190002 -0.770004 -3.649994
2016-06-18   0.008630  23.579987  13.829987 -0.559998 -1.110001
2016-06-19   0.008630  22.320007  10.699982 -1.570007 -0.190002
2016-06-20   0.017260  20.920013  10.040009 -1.449997  1.580002
2016-06-21   0.768956  19.459991  14.720001 -3.460007  2.669998
2016-06-22   0.069100  21.589996  15.320007 -2.190002 -2.589996
2016-06-23   0.008630  20.260010  15.850006  0.770004 -1.339996
2016-06-24   0.008630  19.940002  14.529999 -1.930008  2.750000
2016-06-25   0.069100  20.769989  15.240021 -4.869995  0.440002
2016-06-26   0.008630  20.760010  14.920013 -3.949997 -0.699997
2016-06-27   0.008630  21.000000  14.720001 -2.869995 -0.970001
2016-06-28   0.008630  21.410004  14.480011 -2.899994 -0.839996
2016-06-29   0.008630  22.430023  16.430023 -2.419998 -2.199997
2016-06-30   0.008630  22.490021  14.339996 -2.040009 -0.479996

[182 rows x 5 columns]

In [37]:
df.to_csv('./data/modeldata.csv', sep=';')