Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 100
Kernel: Python 3 (Anaconda)

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.

%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.

data=xr.open_dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_1p00/gfs20170312/gfs_1p00_12z")
data
<xarray.Dataset> Dimensions: (lat: 181, lev: 31, lon: 360, time: 33) Coordinates: * time (time) datetime64[ns] 2017-03-12T12:00:00 2017-03-13 ... * lev (lev) float64 1e+03 975.0 950.0 925.0 900.0 850.0 800.0 ... * lat (lat) float64 -90.0 -89.0 -88.0 -87.0 -86.0 -85.0 -84.0 ... * lon (lon) float64 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 ... Data variables: absvprs (time, lev, lat, lon) float64 ... no4lftxsfc (time, lat, lon) float64 ... no5wavh500mb (time, lat, lon) float64 ... acpcpsfc (time, lat, lon) float64 ... albdosfc (time, lat, lon) float64 ... apcpsfc (time, lat, lon) float64 ... capesfc (time, lat, lon) float64 ... cape180_0mb (time, lat, lon) float64 ... cape255_0mb (time, lat, lon) float64 ... cfrzrsfc (time, lat, lon) float64 ... cicepsfc (time, lat, lon) float64 ... cinsfc (time, lat, lon) float64 ... cin180_0mb (time, lat, lon) float64 ... cin255_0mb (time, lat, lon) float64 ... clwmrprs (time, lev, lat, lon) float64 ... cpofpsfc (time, lat, lon) float64 ... cpratsfc (time, lat, lon) float64 ... crainsfc (time, lat, lon) float64 ... csnowsfc (time, lat, lon) float64 ... cwatclm (time, lat, lon) float64 ... cworkclm (time, lat, lon) float64 ... dlwrfsfc (time, lat, lon) float64 ... dpt2m (time, lat, lon) float64 ... dswrfsfc (time, lat, lon) float64 ... fldcpsfc (time, lat, lon) float64 ... gfluxsfc (time, lat, lon) float64 ... gustsfc (time, lat, lon) float64 ... hgtsfc (time, lat, lon) float64 ... hgtprs (time, lev, lat, lon) float64 ... hgt2pv (time, lat, lon) float64 ... hgtneg2pv (time, lat, lon) float64 ... hgttop0c (time, lat, lon) float64 ... hgt0c (time, lat, lon) float64 ... hgtmwl (time, lat, lon) float64 ... hgttrop (time, lat, lon) float64 ... hindexsfc (time, lat, lon) float64 ... hlcy3000_0m (time, lat, lon) float64 ... hpblsfc (time, lat, lon) float64 ... icahtmwl (time, lat, lon) float64 ... icahttrop (time, lat, lon) float64 ... icecsfc (time, lat, lon) float64 ... icsevprs (time, lev, lat, lon) float64 ... landsfc (time, lat, lon) float64 ... lftxsfc (time, lat, lon) float64 ... lhtflsfc (time, lat, lon) float64 ... msletmsl (time, lat, lon) float64 ... o3mrprs (time, lev, lat, lon) float64 ... pevprsfc (time, lat, lon) float64 ... plpl255_0mb (time, lat, lon) float64 ... potsig995 (time, lat, lon) float64 ... pratesfc (time, lat, lon) float64 ... preslclb (time, lat, lon) float64 ... preslclt (time, lat, lon) float64 ... presmclb (time, lat, lon) float64 ... presmclt (time, lat, lon) float64 ... preshclb (time, lat, lon) float64 ... preshclt (time, lat, lon) float64 ... pressfc (time, lat, lon) float64 ... pres80m (time, lat, lon) float64 ... pres2pv (time, lat, lon) float64 ... presneg2pv (time, lat, lon) float64 ... prescclb (time, lat, lon) float64 ... prescclt (time, lat, lon) float64 ... presmwl (time, lat, lon) float64 ... prestrop (time, lat, lon) float64 ... prmslmsl (time, lat, lon) float64 ... pwatclm (time, lat, lon) float64 ... rhprs (time, lev, lat, lon) float64 ... rh2m (time, lat, lon) float64 ... rhsg330_1000 (time, lat, lon) float64 ... rhsg440_1000 (time, lat, lon) float64 ... rhsg720_940 (time, lat, lon) float64 ... rhsg440_720 (time, lat, lon) float64 ... rhsig995 (time, lat, lon) float64 ... rh30_0mb (time, lat, lon) float64 ... rhclm (time, lat, lon) float64 ... rhtop0c (time, lat, lon) float64 ... rh0c (time, lat, lon) float64 ... shtflsfc (time, lat, lon) float64 ... snodsfc (time, lat, lon) float64 ... soilw0_10cm (time, lat, lon) float64 ... soilw10_40cm (time, lat, lon) float64 ... soilw40_100cm (time, lat, lon) float64 ... soilw100_200cm (time, lat, lon) float64 ... spfh2m (time, lat, lon) float64 ... spfh80m (time, lat, lon) float64 ... spfh30_0mb (time, lat, lon) float64 ... sunsdsfc (time, lat, lon) float64 ... tcdcclm (time, lat, lon) float64 ... tcdcblcll (time, lat, lon) float64 ... tcdclcll (time, lat, lon) float64 ... tcdcmcll (time, lat, lon) float64 ... tcdchcll (time, lat, lon) float64 ... tcdcccll (time, lat, lon) float64 ... tmax2m (time, lat, lon) float64 ... tmin2m (time, lat, lon) float64 ... tmplclt (time, lat, lon) float64 ... tmpmclt (time, lat, lon) float64 ... tmphclt (time, lat, lon) float64 ... tmpsfc (time, lat, lon) float64 ... tmpprs (time, lev, lat, lon) float64 ... tmp_1829m (time, lat, lon) float64 ... tmp_2743m (time, lat, lon) float64 ... tmp_3658m (time, lat, lon) float64 ... tmp2m (time, lat, lon) float64 ... tmp80m (time, lat, lon) float64 ... tmp100m (time, lat, lon) float64 ... tmpsig995 (time, lat, lon) float64 ... tmp30_0mb (time, lat, lon) float64 ... tmp2pv (time, lat, lon) float64 ... tmpneg2pv (time, lat, lon) float64 ... tmpmwl (time, lat, lon) float64 ... tmptrop (time, lat, lon) float64 ... tozneclm (time, lat, lon) float64 ... tsoil0_10cm (time, lat, lon) float64 ... tsoil10_40cm (time, lat, lon) float64 ... tsoil40_100cm (time, lat, lon) float64 ... tsoil100_200cm (time, lat, lon) float64 ... ugwdsfc (time, lat, lon) float64 ... uflxsfc (time, lat, lon) float64 ... ugrdprs (time, lev, lat, lon) float64 ... ugrd_1829m (time, lat, lon) float64 ... ugrd_2743m (time, lat, lon) float64 ... ugrd_3658m (time, lat, lon) float64 ... ugrd10m (time, lat, lon) float64 ... ugrd80m (time, lat, lon) float64 ... ugrd100m (time, lat, lon) float64 ... ugrdsig995 (time, lat, lon) float64 ... ugrd30_0mb (time, lat, lon) float64 ... ugrd2pv (time, lat, lon) float64 ... ugrdneg2pv (time, lat, lon) float64 ... ugrdpbl (time, lat, lon) float64 ... ugrdmwl (time, lat, lon) float64 ... ugrdtrop (time, lat, lon) float64 ... ulwrfsfc (time, lat, lon) float64 ... ulwrftoa (time, lat, lon) float64 ... ustm6000_0m (time, lat, lon) float64 ... uswrfsfc (time, lat, lon) float64 ... uswrftoa (time, lat, lon) float64 ... vgwdsfc (time, lat, lon) float64 ... vflxsfc (time, lat, lon) float64 ... vgrdprs (time, lev, lat, lon) float64 ... vgrd_1829m (time, lat, lon) float64 ... vgrd_2743m (time, lat, lon) float64 ... vgrd_3658m (time, lat, lon) float64 ... vgrd10m (time, lat, lon) float64 ... vgrd80m (time, lat, lon) float64 ... vgrd100m (time, lat, lon) float64 ... vgrdsig995 (time, lat, lon) float64 ... vgrd30_0mb (time, lat, lon) float64 ... vgrd2pv (time, lat, lon) float64 ... vgrdneg2pv (time, lat, lon) float64 ... vgrdpbl (time, lat, lon) float64 ... vgrdmwl (time, lat, lon) float64 ... vgrdtrop (time, lat, lon) float64 ... vratepbl (time, lat, lon) float64 ... vstm6000_0m (time, lat, lon) float64 ... vvelprs (time, lev, lat, lon) float64 ... vvelsig995 (time, lat, lon) float64 ... vwsh2pv (time, lat, lon) float64 ... vwshneg2pv (time, lat, lon) float64 ... vwshtrop (time, lat, lon) float64 ... watrsfc (time, lat, lon) float64 ... weasdsfc (time, lat, lon) float64 ... wiltsfc (time, lat, lon) float64 ... var00212m (time, lat, lon) float64 ... Attributes: title: GFS 1.0 deg starting from 12Z12mar2017, downloaded Mar 12 16:54 UTC Conventions: COARDS GrADS dataType: Grid history: Sun Mar 12 18:17:19 UTC 2017 : imported by GrADS Data Server 2.0

Conveniently store the coordinate variables:

time_1d = data.time
lat_1d = data.lat
lon_1d = data.lon

time_1d = data['time']
<xarray.DataArray 'time' (time: 33)> array(['2017-03-12T12:00:00.000000000', '2017-03-13T00:00:00.000000000', '2017-03-13T12:00:00.000000000', '2017-03-14T00:00:00.000000000', '2017-03-14T12:00:00.000000000', '2017-03-15T00:00:00.000000000', '2017-03-15T12:00:00.000000000', '2017-03-16T00:00:00.000000000', '2017-03-16T12:00:00.000000000', '2017-03-17T00:00:00.000000000', '2017-03-17T12:00:00.000000000', '2017-03-18T00:00:00.000000000', '2017-03-18T12:00:00.000000000', '2017-03-19T00:00:00.000000000', '2017-03-19T12:00:00.000000000', '2017-03-20T00:00:00.000000000', '2017-03-20T12:00:00.000000000', '2017-03-21T00:00:00.000000000', '2017-03-21T12:00:00.000000000', '2017-03-22T00:00:00.000000000', '2017-03-22T12:00:00.000000000', '2017-03-23T00:00:00.000000000', '2017-03-23T12:00:00.000000000', '2017-03-24T00:00:00.000000000', '2017-03-24T12:00:00.000000000', '2017-03-25T00:00:00.000000000', '2017-03-25T12:00:00.000000000', '2017-03-26T00:00:00.000000000', '2017-03-26T12:00:00.000000000', '2017-03-27T00:00:00.000000000', '2017-03-27T12:00:00.000000000', '2017-03-28T00:00:00.000000000', '2017-03-28T12:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 2017-03-12T12:00:00 2017-03-13 ... Attributes: grads_dim: t grads_mapping: linear grads_size: 33 grads_min: 12z12mar2017 grads_step: 12hr long_name: time minimum: 12z12mar2017 maximum: 12z28mar2017 resolution: 0.5

Let's grab temperature on isobaric surfaces (temperature at a constant pressure value in the atmosphere). What are the dimensions?

temperature=data['tmpprs']
temperature = data.tmpprs
temperature
<xarray.DataArray 'tmpprs' (time: 33, lev: 31, lat: 181, lon: 360)> [66658680 values with dtype=float64] Coordinates: * time (time) datetime64[ns] 2017-03-12T12:00:00 2017-03-13 ... * lev (lev) float64 1e+03 975.0 950.0 925.0 900.0 850.0 800.0 750.0 ... * lat (lat) float64 -90.0 -89.0 -88.0 -87.0 -86.0 -85.0 -84.0 -83.0 ... * lon (lon) float64 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 ... Attributes: long_name: ** (1000 975 950 925 900.. 7 5 3 2 1) temperature [k]

Let's grab the map at the initial time, and at 850 hPa.

temp_850=temperature.sel(time=time_1d[0],lev=850.,method='nearest')
temp_850
<xarray.DataArray 'tmpprs' (lat: 181, lon: 360)> array([[ 235.69999695, 235.69999695, 235.69999695, ..., 235.69999695, 235.69999695, 235.69999695], [ 235.30000305, 235.30000305, 235.30000305, ..., 235.19999695, 235.19999695, 235.19999695], [ 236.30000305, 236.30000305, 236.30000305, ..., 236.1000061 , 236.19999695, 236.19999695], ..., [ 253.40000916, 253.40000916, 253.40000916, ..., 253.40000916, 253.40000916, 253.40000916], [ 253.80000305, 253.80000305, 253.80000305, ..., 253.80000305, 253.80000305, 253.80000305], [ 253.1000061 , 253.1000061 , 253.1000061 , ..., 253.1000061 , 253.1000061 , 253.1000061 ]]) Coordinates: time datetime64[ns] 2017-03-12T12:00:00 lev float64 850.0 * lat (lat) float64 -90.0 -89.0 -88.0 -87.0 -86.0 -85.0 -84.0 -83.0 ... * lon (lon) float64 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 ... Attributes: long_name: ** (1000 975 950 925 900.. 7 5 3 2 1) temperature [k] 
%matplotlib inline
import numpy as np
import pandas as pd 

import matplotlib.pyplot as plt

plt.pcolor(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()
Image in a Jupyter notebook

Now let's use Basemap to generate a contour map over the US!

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()
Image in a Jupyter notebook

Let's do a filled contour instead!

%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='h', 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()
Image in a Jupyter notebook

Let's change the contour interval and the colormap!

%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()
Image in a Jupyter notebook

Lines and contours. Why not?

%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()
Image in a Jupyter notebook

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.

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')

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!

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.show()
Image in a Jupyter notebook

Time-height cross sections

Here, we want values at a specific point over time.

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')

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()
Image in a Jupyter notebook
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()
/projects/anaconda3/lib/python3.5/site-packages/matplotlib/contour.py:370: RuntimeWarning: invalid value encountered in true_divide dist = np.add.reduce(([(abs(s)[i] / L[i]) for i in range(xsize)]), -1)
Image in a Jupyter notebook