CoCalc Public FilesATMS-305-Content / Week-9 / Week-9-Handouts / March 15 in class exercise.ipynb
Author: Steve Nesbitt
Views : 53
Compute Environment: Ubuntu 18.04 (Deprecated)

# ATMS 305: March 15 in-class exercise

In this exercise, we're going to add to your xarray, matplotlib, and basemap skill set.

## Colormaps

Let's bring back our map of GFS 850 hPa temperatures at the first time in the model simulation.

In [2]:
%pylab inline
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import xarray as xr

#data prep
#grab coordinate values
time_1d = data.time # same as data['time']
lat_1d = data.lat # data['lat']
lon_1d = data.lon
#what are we doing here?
temp_850 = data['tmpprs'].sel(time=time_1d[0],lev=850.,method='nearest')
#make the lats and lons two dimensional like the data
lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)

Populating the interactive namespace from numpy and matplotlib

Now plot!

In [3]:
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)

map_fig.drawcoastlines()

map_fig.drawstates()

map_fig.drawcountries(linewidth=1.5)

plt.colorbar()

# # Make a title with the time value
plt.title('850 hPa temperature forecast ' + '(' + u'\u00b0' + 'C)' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()


Let's say that we didn't like the colormap in the above plot. matplotlib has many to choose from. Run the following code to display them all.

In [4]:
cmaps = [('Perceptually Uniform Sequential',
['viridis', 'inferno', 'plasma', 'magma']),
('Sequential',     ['Blues', 'BuGn', 'BuPu',
'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
'copper', 'gist_heat', 'gray', 'hot',
'pink', 'spring', 'summer', 'winter']),
('Diverging',      ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
'seismic']),
('Qualitative',    ['Accent', 'Dark2', 'Paired', 'Pastel1',
'Pastel2', 'Set1', 'Set2', 'Set3']),
('Miscellaneous',  ['gist_earth', 'terrain', 'ocean', 'gist_stern',
'brg', 'CMRmap', 'cubehelix',
'gnuplot', 'gnuplot2', 'gist_ncar',
'nipy_spectral', 'jet', 'rainbow',
'gist_rainbow', 'hsv', 'flag', 'prism'])]

nrows = max(len(cmap_list) for cmap_category, cmap_list in cmaps)
print(nrows)

max([len(cmap_list) for cmap_category, cmap_list in cmaps])
# [len(cmap_list) for cmap_category, cmap_list in cmaps]

18
18
In [ ]:
# Have colormaps separated into categories:
# http://matplotlib.org/examples/color/colormaps_reference.html
cmaps = [('Perceptually Uniform Sequential',
['viridis', 'inferno', 'plasma', 'magma']),
('Sequential',     ['Blues', 'BuGn', 'BuPu',
'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
'copper', 'gist_heat', 'gray', 'hot',
'pink', 'spring', 'summer', 'winter']),
('Diverging',      ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
'seismic']),
('Qualitative',    ['Accent', 'Dark2', 'Paired', 'Pastel1',
'Pastel2', 'Set1', 'Set2', 'Set3']),
('Miscellaneous',  ['gist_earth', 'terrain', 'ocean', 'gist_stern',
'brg', 'CMRmap', 'cubehelix',
'gnuplot', 'gnuplot2', 'gist_ncar',
'nipy_spectral', 'jet', 'rainbow',
'gist_rainbow', 'hsv', 'flag', 'prism'])]

nrows = max(len(cmap_list) for cmap_category, cmap_list in cmaps)

fig, axes = plt.subplots(nrows=nrows)
axes[0].set_title(cmap_category + ' colormaps', fontsize=14)

for ax, name in zip(axes, cmap_list):
pos = list(ax.get_position().bounds)
x_text = pos[0] - 0.01
y_text = pos[1] + pos[3]/2.
fig.text(x_text, y_text, name, va='center', ha='right', fontsize=10)

# Turn off *all* ticks & spines, not just the ones with colormaps.
for ax in axes:
ax.set_axis_off()

for cmap_category, cmap_list in cmaps:

plt.show()


If you want to select one from the available ones, you can choose one by entering the name of the colormap you want.

In [7]:
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='bwr')

map_fig.drawcoastlines()

map_fig.drawstates()

map_fig.drawcountries(linewidth=1.5)

plt.colorbar()

# Make a title with the time value
plt.title('850 hPa temperature forecast ' + '(' + u'\u00b0' + 'C)' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()


You can reverse the colormap simply by putting the suffix _r on the colortable.

In [ ]:
%pylab inline
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import xarray as xr
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('850 hPa temperature forecast ' + '(' + u'\u00b0' + 'C)' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()

Populating the interactive namespace from numpy and matplotlib

## Vectors, barbs, streamlines

In many fields, including atmospheric sciences, it is desirable to plot vector fields on a map or plot. Let's work on plotting winds at the 850 hPa level using vectors, barbs, and streamlines.

First, let's grab our vector components, which are u and v winds (east-west and north-south wind vector components, respectively).

In [ ]:
uwnd = data['ugrdprs'].sel(time=time_1d[0],lev=850.,method='nearest')
vwnd = data['vgrdprs'].sel(time=time_1d[0],lev=850.,method='nearest')

In [15]:
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

Q = map_fig.quiver(lon_2d, lat_2d, uwnd, vwnd,scale=500)

qk = plt.quiverkey(Q, 0.95, 0.98, 10, '10 m/s', labelpos='W')

map_fig.drawcoastlines()

map_fig.drawstates()

map_fig.drawcountries(linewidth=1.5)

# Make a title with the time value
plt.title('850 hPa wind forecast ' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()


Barbs - any meteorologist's dream! Do you know how to interpret them? What are the units (compare with above figure)?

In [16]:
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

Q = map_fig.barbs(lon_2d, lat_2d, uwnd, vwnd)

map_fig.drawcoastlines()

map_fig.drawstates()

map_fig.drawcountries(linewidth=1.5)

# Make a title with the time value
plt.title('850 hPa wind forecast ' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()


Finally, streamlines. For visualizing the flow. Might take a minute or three to compute.

In [17]:
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

Q = map_fig.streamplot(lon_2d, lat_2d, uwnd, vwnd, density=10.)

map_fig.drawcoastlines()

map_fig.drawstates()

map_fig.drawcountries(linewidth=1.5)

# Make a title with the time value
plt.title('850 hPa wind forecast ' + ' for '  + str(time_1d.values[0])[0:13] + 'UTC', fontsize=10)

plt.show()


## Subsetting & mathematical operations

Before we leave this topic, we should talk about some more advanced operations in xarray. The data model is very similar to pandas in that we can do mathematical operations on Datasets and variables easily. This can be useful for doing geospatial calculations.

Let's say we wanted to plot a time series of the mean and standard deviation of the 500 hPa temperature averaged over a box covering the eastern United States.

We can select a "box" of data using the slice command in the sel class of a Dataset.

In [18]:
box=data.sel(lon=slice(260,290),lat=slice(30,45),lev=500.)
box

<xarray.Dataset> Dimensions: (lat: 16, lon: 31, time: 33) Coordinates: * time (time) datetime64[ns] 2017-03-15 2017-03-15T12:00:00 ... lev float64 500.0 * lat (lat) float64 30.0 31.0 32.0 33.0 34.0 35.0 36.0 37.0 ... * lon (lon) float64 260.0 261.0 262.0 263.0 264.0 265.0 266.0 ... Data variables: absvprs (time, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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 00Z15mar2017, downloaded Mar 15 04:55 UTC Conventions: COARDS GrADS dataType: Grid history: Wed Mar 15 11:17:27 GMT 2017 : imported by GrADS Data Server 2.0

Let's do a spatial average over our box. We could average all of the fields in the same way, but that might take awhile:

box.mean(dim=['lat','lon'])


Instead, let's just do it over a specific variable we want.

In [19]:
tser_mean_850=box['tmpprs'].mean(dim=['lat','lon'])
tser_mean_850

<xarray.DataArray 'tmpprs' (time: 33)> array([ 250.61431829, 249.75867339, 251.20020562, 251.75121384, 253.16472266, 254.21068994, 254.66694007, 254.15302892, 254.89718259, 255.38972344, 256.68327147, 256.94355334, 257.04173897, 256.27843202, 255.28468221, 254.57036769, 255.67681946, 255.84798828, 255.98548846, 255.51472215, 255.95847271, 255.49314902, 255.60827055, 255.61351324, 256.00786772, 255.27339163, 255.82097306, 255.89173886, 256.95786833, 256.91089199, 257.98266678, 258.04214238, 258.59960175]) Coordinates: * time (time) datetime64[ns] 2017-03-15 2017-03-15T12:00:00 2017-03-16 ... lev float64 500.0
In [20]:
tser_std_850=box['tmpprs'].std(dim=['lat','lon'])
tser_std_850

<xarray.DataArray 'tmpprs' (time: 33)> array([ 7.86333065, 7.54476039, 6.39708162, 6.24138237, 6.03436783, 4.46145189, 3.30826857, 4.43383303, 5.73906782, 5.76284101, 4.60056897, 2.38810977, 2.39697762, 2.763669 , 4.43962394, 4.98373576, 4.44094062, 4.0185512 , 3.37050153, 2.78676554, 2.5750156 , 2.57954858, 2.62372798, 2.81670727, 3.03142215, 3.11426625, 3.45385701, 3.38491711, 2.92360837, 2.77942626, 2.83550817, 2.29379907, 2.71320019]) Coordinates: * time (time) datetime64[ns] 2017-03-15 2017-03-15T12:00:00 2017-03-16 ... lev float64 500.0

Now, let's plot: (we'll use the standard deviation as an error bar)

In [21]:
%pylab inline

plt.figure(figsize=(11,4))

plt.errorbar(tser_mean_850.time.values, tser_mean_850.values, yerr=tser_std_850.values)


Populating the interactive namespace from numpy and matplotlib
/projects/anaconda3/lib/python3.5/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['box'] %matplotlib prevents importing * from pylab and numpy "\n%matplotlib prevents importing * from pylab and numpy"
<Container object of 3 artists>

One more trick! Let's say we wanted to average over each day. No problem! Note that the freq keyword is the same as in pandas - so it should be familiar.

In [23]:
tser_dailymean_850 = tser_mean_850.resample(dim='time',how='mean', freq='D') # now monthly
tser_dailystd_850 = tser_mean_850.resample(dim='time',how='std', freq='D')


In [24]:
%pylab inline

plt.figure(figsize=(11,4))
#note change in time variable, since we resampled
plt.errorbar(tser_dailymean_850.time.values, tser_dailymean_850.values, yerr=tser_dailystd_850.values)