ATMS 305: More working with matplotlib and basemap¶
Since basemap is part of matplotlib, we can use the mapping capabilities of basemap along with the plotting capabilities of matplotlib for ultimate in geo-geekiness. This capability is very useful for visualizing geoscience data.
In [2]:
%matplotlib inline
import xarray as xr
Now, let's open a file with xarray, and get the variables available. You can change the time and the run as you wish in the URL. These files are hosted at NOAA/NWS/NCEP, and updated in real time.
In [4]:
data=xr.open_dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_1p00/gfs20170329/gfs_1p00_00z")
data['capesfc']
Out[4]:
Conveniently store the coordinate variables:
In [7]:
time_1d = data.time
lat_1d = data.lat
lon_1d = data.lon
pres_1d = data.lev
time_1d = data['time']
In [8]:
time_1d
Out[8]:
Let's grab temperature on isobaric surfaces (temperature at a constant pressure value in the atmosphere). What are the dimensions?
In [9]:
temperature=data['tmpprs']
temperature = data.tmpprs
temperature
Out[9]:
Let's grab the map at the initial time, and at 850 hPa.
In [10]:
temp_850=temperature.sel(time=time_1d[0],lev=850.,method='nearest')
temp_850
Out[10]:
In [11]:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.pcolormesh(lon_1d,lat_1d,temp_850)
plt.xlabel('Longitude (deg)')
plt.ylabel('Latitude (deg)')
title = 'GFS Temperature - 850 hPa at {0}'.format(str(time_1d.values[0])[:-13]) # remove the useless 0s
plt.title(title)
cb = plt.colorbar()
cb.set_label('Temperature (K)')
plt.tight_layout()
Now let's use Basemap to generate a contour map over the US!
In [23]:
# Step 1: import Basemap
from mpl_toolkits.basemap import Basemap
# Step 2: define map
map_fig = Basemap(llcrnrlat=22.,llcrnrlon=360-125.,urcrnrlat=50.,urcrnrlon=360.-65.)
# Step 2.5: data prep
lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
# Step 3: plot the map
map_fig.pcolormesh(lon_1d,lat_1d,temp_850,cmap='jet')
#Step 4: make it pretty
map_fig.drawcoastlines()
map_fig.drawstates()
map_fig.drawcountries()
plt.colorbar()
plt.show()
In [24]:
from mpl_toolkits.basemap import Basemap
fig=plt.figure(figsize=(11,8.5)) #width, height in inches
map_fig = Basemap(llcrnrlon=260.,llcrnrlat=30.,urcrnrlon=280.,urcrnrlat=50.,
resolution='l', projection='cyl')
# #need to mesh the latitudes and longitudes
lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
#store the plotting object into CS, and then add contour labels
#also convert to degrees C
CS=map_fig.contour(lon_2d, lat_2d, temp_850-273.15)
plt.clabel(CS,inline=1, fontsize=14, fmt='%1.0f')
map_fig.drawcoastlines()
map_fig.drawstates()
map_fig.drawcountries(linewidth=1.5)
# Make a title with the time value
# plt.title('Temperature forecast ' + '(' + u'\u00b0' + 'C)' + ' for ' + str(time_1d.values[0])[0:13] + 'UTC', fontsize=20)
plt.title('Temperature Forecast ({0} C) for {1}'.format(u'\u00b0', str(time_1d.values[0])[0:13]))
plt.show()
Let's do a filled contour instead!
In [30]:
%matplotlib inline
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(11,8.5)) #width, height in inches
map_fig = Basemap(llcrnrlon=260.,llcrnrlat=30.,urcrnrlon=280.,urcrnrlat=50.,
resolution='l', projection='cyl')
#use contourf to fill the contours
CS=map_fig.contourf(lon_2d, lat_2d, temp_850-273.15,cmap='RdBu_r')
map_fig.drawcoastlines()
map_fig.drawstates()
map_fig.drawcountries(linewidth=1.5)
plt.colorbar()
# Make a title with the time value
plt.title('Temperature forecast ' + '(' + u'\u00b0' + 'F)' + ' for ' + str(time_1d.values[0])[0:13] + 'UTC', fontsize=20)
plt.show()
Let's change the contour interval and the colormap!
In [26]:
%matplotlib inline
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(11,8.5)) #width, height in inches
map = Basemap(llcrnrlon=260.,llcrnrlat=30.,urcrnrlon=280.,urcrnrlat=50.,
resolution='l', projection='cyl')
#specify levels, and change the colortable
#see here for more color tables http://matplotlib.org/examples/color/colormaps_reference.html
#add _r to reverse color table
levels=np.linspace(-32,32,100)
map.contourf(lon_2d, lat_2d, temp_850-273.15,levels,cmap='seismic')
map.drawcoastlines()
map.drawstates()
map.drawcountries(linewidth=1.5)
plt.colorbar()
# Make a title with the time value
plt.title('Temperature forecast ' + '(' + u'\u00b0' + 'F)' + ' for ' + str(time_1d.values[0])[0:13] + 'UTC', fontsize=20)
plt.show()
Lines and contours. Why not?
In [27]:
%matplotlib inline
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(11,8.5)) #width, height in inches
map = Basemap(llcrnrlon=260.,llcrnrlat=30.,urcrnrlon=280.,urcrnrlat=50.,
resolution='l', projection='cyl')
#use pcolormesh to plot a blocky representation where data is defined,
#note the change in how to specify values
PM=map.pcolormesh(lon_2d, lat_2d, temp_850-273.15,cmap='seismic',vmin=-32,vmax=32)
CS=map.contour(lon_2d, lat_2d, temp_850-273.15,colors='k')
plt.clabel(CS,inline=0, fontsize=14, fmt='%1.0f')
map.drawcoastlines()
map.drawstates()
map.drawcountries(linewidth=1.5)
plt.colorbar(PM)
# Make a title with the time value
plt.title('Temperature forecast ' + '(' + u'\u00b0' + 'F)' + ' for ' + str(time_1d.values[0])[0:13] + 'UTC', fontsize=20)
plt.show()
Cross-sections¶
We can also draw vertical cross sections using the vertical coordinate as the y-axis and latitude and longitude as the x-axis.
In [29]:
time_1d = data.time
lat_1d = data.lat
lon_1d = data.lon
lev_1d = data.lev
#Get the temperature in degrees Kelvin
temp_xsect=data['tmpprs'].sel(time=time_1d[0],lon=-88+360,method='nearest')
#Get the east-west wind in m/s
uwnd_xsect=data['ugrdprs'].sel(time=time_1d[0],lon=-88+360,method='nearest')
temp_xsect
Out[29]:
Here, we don't need a map. We can just use matplotlib. Lets contourf the temperature and contour the u-wind. Meteorologists - see the thermal wind in action!
In [32]:
fig,ax = plt.subplots(1,1)
# color map
levels=np.arange(-100,40,5)
CS=ax.contourf(lat_1d,lev_1d,temp_xsect-273.15,levels)
ax.set_ylim([1000,100])
plt.colorbar(CS)
# lines on the map
levels2=np.arange(-200,200,20)
CS2=ax.contour(lat_1d,lev_1d,uwnd_xsect,levels2,colors='k')
ax.clabel(CS2,inline=1, fontsize=12, fmt='%1.0f')
plt.xlabel('Latitude')
plt.ylabel('Pressue (hPa)')
plt.show()
Time-height cross sections¶
Here, we want values at a specific point over time.
In [33]:
time_1d = data.time.values
lat_1d = data.lat
lon_1d = data.lon
lev_1d = data.lev
#Get the relative humidity in %
rh_xsect=data['rhprs'].sel(lat=40.,lon=-88+360,method='nearest')
#Get the absolute vorticity in s**-1
absvprs_xsect=data['absvprs'].sel(lat=40.,lon=-88+360,method='nearest')
In [34]:
fig,ax = plt.subplots(1,1, figsize=(10, 10))
levels=np.arange(0,35e-5,5e-5)
CS=ax.contourf(lev_1d,time_1d,absvprs_xsect,levels)
plt.colorbar(CS)
levels2=np.arange(0,110,20)
CS2=ax.contour(lev_1d,time_1d,rh_xsect,levels2,colors='k')
ax.clabel(CS2,inline=1, fontsize=12, fmt='%1.0f')
plt.show()
In [35]:
fig,ax = plt.subplots(1,1)
fig.set_size_inches([15,5])
levels=np.arange(0,35e-5,5e-5)
CS=ax.contourf(time_1d,lev_1d,np.transpose(absvprs_xsect.values),levels)
ax.set_ylim([1000,100])
plt.colorbar(CS)
levels2=np.arange(0,110,20)
CS2=ax.contour(time_1d,lev_1d,np.transpose(rh_xsect.values),levels2,colors='k')
ax.clabel(CS2,inline=1, fontsize=12, fmt='%1.0f')
plt.show()
In [ ]:
plt.contour() #shift+tab gives you information on the command