In [1]:
import pygrib
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
In [2]:
grbs = pygrib.open('data/gfs.t00z.pgrb2.0p25.f006')
The file can be iterated to show all of the grib messages. In each message there are informations like name, units, type of level, level, valid date and forecast date
In [3]:
for grb in grbs:
print grb
Each of those messages have a set of attributes, or keys, which can be used for searching
In [4]:
for key in grb.keys():
print key
Those keys can be used to filter variables. If I want the 10 metre above the ground wind, I can use the name key, the typeOfLevel and level. This will return a list of grib messages that have been found. Using the [0] selects the first element in that list
In [5]:
u10_grb = grbs.select(name='10 metre U wind component', typeOfLevel='heightAboveGround', level=10)[0]
v10_grb = grbs.select(name='10 metre V wind component', typeOfLevel='heightAboveGround', level=10)[0]
t2_grb = grbs.select(name='2 metre temperature', typeOfLevel='heightAboveGround', level=2)[0]
The key values can be used in my object to return a numpy array of data
In [6]:
u10 = u10_grb.values
v10 = v10_grb.values
t2 = t2_grb.values
print t2.shape, u10.shape, v10.shape
Using grib_message.latslons will return a 2D numpy array of latitudes and longitudes
In [7]:
lats, lons = u10_grb.latlons()
print lats.shape, lons.shape, lats.min(), lats.max(), lons.min(), lons.max()
In [8]:
print lats
print lons
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 [9]:
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 functions to make the temperature and wind plot.
In [10]:
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 = 20
cs = m.contourf(lons, lats, t2)
qv = m.quiver(lons[::skip, ::skip], lats[::skip, ::skip], u10[::skip, ::skip], v10[::skip, ::skip])
There are many different projections to choose from in Basemap:
In [11]:
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 = 20
x, y = m(lons, lats)
cs = m.contourf(x, y, t2)
qv = m.quiver(x[::skip, ::skip], y[::skip, ::skip], u10[::skip, ::skip], v10[::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 = 20
x, y = m(lons, lats)
cs = m.contourf(x, y, t2)
qv = m.quiver(x[::skip, ::skip], y[::skip, ::skip], u10[::skip, ::skip], v10[::skip, ::skip])
As my data is sorted from 0 to 360 in longitude, you need to reorder them from -180 to 180. To do so we can create a functions called flip_grid that takes the variable that we want to reorder and my longitude 2D array. An if statment can be used to check if the variable is 2D or 3D.
In [12]:
def flip_grid(var, lons):
fltr = lons[0] >= 180
newlons = np.concatenate(((lons - 360)[:, fltr], lons[:, ~fltr]), axis=-1)
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)
return newvar, newlons
The flip_grid function is used to reorder u10, v10 and t2 variables, and the data is ploted
In [13]:
u10, newlons = flip_grid(u10, lons)
v10, newlons = flip_grid(v10, lons)
t2, newlons = flip_grid(t2, lons)
In [14]:
fig = plt.figure(figsize=(16,35))
m = Basemap(llcrnrlon=newlons.min(), llcrnrlat=lats.min(), urcrnrlon=newlons.max(), urcrnrlat=lats.max())
m.drawcoastlines()
m.drawcountries()
skip = 20
cs = m.contourf(newlons, lats, t2)
qv = m.quiver(newlons[::skip, ::skip], lats[::skip, ::skip], u10[::skip, ::skip], v10[::skip, ::skip])
To make a more regional plot, we can set different boundary lattitudes and longitudes. Also, we are able to use most of the matplotlib tools to make a more sofisticated plot
In [15]:
fig = plt.figure(figsize=(12,25))
m = Basemap(llcrnrlon=-100, llcrnrlat=-60, urcrnrlon=-20, urcrnrlat=10, resolution='i')
parallels = np.arange(-90, 90, 10)
meridians = np.arange(0, 360, 10)
m.drawcoastlines()
m.drawcountries()
m.drawstates() # bad idea!!! >.<
m.drawmeridians(meridians, labels=[0, 0, 0, 1])
m.drawparallels(parallels, labels=[1, 0, 0, 0])
skip = 10
plt.title('Temp. a 2 metros e Vento a 10 metros: %s' % u10_grb.validDate)
cs = m.contourf(newlons, lats, t2 - 273.15, cmap='RdYlBu_r')
cbar = m.colorbar(cs)
cbar.set_label(r'Temperatura ($^{o}$C)')
qv = m.quiver(newlons[::skip, ::skip], lats[::skip, ::skip], u10[::skip, ::skip], v10[::skip, ::skip], scale=300)
qk = plt.quiverkey(qv, 0.9, 1.015, 20, '20 m/s', labelpos='W')
We can use the functions readshapefile to load a shapefile with a updated version of the political boundaries of the brazilian states. (For more examples on how to use it, check the following link: http://basemaptutorial.readthedocs.io/en/latest/shapefile.html )
In [16]:
fig = plt.figure(figsize=(12,25))
m = Basemap(llcrnrlon=-100, llcrnrlat=-60, urcrnrlon=-20, urcrnrlat=10, resolution='i')
parallels = np.arange(-90, 90, 10)
meridians = np.arange(0, 360, 10)
m.drawcoastlines()
m.drawcountries()
m.readshapefile('/home/rafaelca/estados_2010_shapefile/estados_2010', 'estado_2010') # much better now :D
m.drawmeridians(meridians, labels=[0, 0, 0, 1])
m.drawparallels(parallels, labels=[1, 0, 0, 0])
skip = 10
plt.title('Temp. a 2 metros e Vento a 10 metros: %s' % u10_grb.validDate)
cs = m.contourf(newlons, lats, t2 - 273.15, cmap='RdYlBu_r')
cbar = m.colorbar(cs)
cbar.set_label(r'Temperatura ($^{o}$C)')
qv = m.quiver(newlons[::skip, ::skip], lats[::skip, ::skip], u10[::skip, ::skip], v10[::skip, ::skip], scale=300)
qk = plt.quiverkey(qv, 0.9, 1.015, 20, '20 m/s', labelpos='W')
To make a cross section of the zonal wind, for example, we need a list of the variable u for each isobaric level
In [17]:
list_of_messages = grbs.select(shortName='u', typeOfLevel='isobaricInhPa')
for message in list_of_messages:
print message
A variable up is then initialized as an empty list so it can be filled with the 2D array of u for each of the isobaric levels. A second variable lev is also created to be filled with the isobaric levels
In [18]:
up = []
lev = []
for grb in list_of_messages:
up.append(grb.values)
lev.append(grb.level)
up = np.array(up)
lev = np.array(lev)
print lev, lev.shape, up.shape
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())
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, newlons, lats)
print newlons[idx[0], idx[1]], lats[idx[0], idx[1]]
Selecting the index idx[1], which means the longitude, we can make a plot for up[:, :, idx[1]] which is the zonal wind for all levels, all latitudes in longitude nearest to -46.616667
In [21]:
fig = plt.figure(figsize=(20,8))
ax = fig.add_subplot(111)
newlats, levs = np.meshgrid(lats[:, idx[1]], lev)
cs = ax.contourf(newlats, levs, up[:, :, idx[1]], cmap='PiYG')
cbar = plt.colorbar(cs)
ax.invert_yaxis()
To read multiple grib files we need to make a for loop to save the variables of interest in a numpy array. The following steps are followed below:
In [23]:
u10 = []
v10 = []
dates = []
import os
list_of_files = os.listdir('data')
for gribfile in sorted(list_of_files):
if gribfile.startswith('gfs.t00z.pgrb2.0p25'):
infile = os.path.join('data', gribfile)
print 'Open file: %s' % infile
grbs = pygrib.open(infile)
u10_grb = grbs.select(name='10 metre U wind component', typeOfLevel='heightAboveGround', level=10)[0]
v10_grb = grbs.select(name='10 metre V wind component', typeOfLevel='heightAboveGround', level=10)[0]
u10.append(u10_grb.values)
v10.append(v10_grb.values)
dates.append(u10_grb.validDate)
grbs.close()
lats, lons = u10_grb.latlons()
u10 = np.array(u10)
v10 = np.array(v10)
Intead of using:
list_of_files = os.listdir('.')
for infile in sorted(list_of_files):
if infile.startswith('gfs.t00z.pgrb2.0p25'):We can use the glob function for glob library:
list_of_files = glob('gfs.t00z.pgrb2.0p25*')
for infile in sorted(list_of_files):
In [24]:
u10 = []
v10 = []
dates = []
from glob import glob
list_of_files = glob('data/gfs.t00z.pgrb2.0p25*')
for infile in sorted(list_of_files):
print 'Open file: %s' % infile
grbs = pygrib.open(infile)
u10_grb = grbs.select(name='10 metre U wind component', typeOfLevel='heightAboveGround', level=10)[0]
v10_grb = grbs.select(name='10 metre V wind component', typeOfLevel='heightAboveGround', level=10)[0]
u10.append(u10_grb.values)
v10.append(v10_grb.values)
dates.append(u10_grb.validDate)
grbs.close()
lats, lons = u10_grb.latlons()
u10 = np.array(u10)
v10 = np.array(v10)
Now the same steps that were performed earlier are repeated now:
In [25]:
u10, newlons = flip_grid(u10, lons)
v10, newlons = flip_grid(v10, lons)
In [26]:
idx = find_nearest(-46.616667, -23.650000, newlons, lats)
Import pandas to create a dataframe with dates as index and then save it to a csv file
In [27]:
import pandas as pd
data = {'u10': u10[:, idx[0], idx[1]], 'v10': v10[:, idx[0], idx[1]]}
df = pd.DataFrame(data, index=dates)
print df
In [28]:
df.to_csv('data.csv', sep=';')