Open in CoCalc with one click!

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=xr.open_dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_1p00/gfs20170315/gfs_1p00_00z") #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) gradient = np.linspace(0, 1, 256) gradient = np.vstack((gradient, gradient)) def plot_color_gradients(cmap_category, cmap_list, nrows): fig, axes = plt.subplots(nrows=nrows) fig.subplots_adjust(top=0.8, bottom=0.01, left=0.2, right=1.5) axes[0].set_title(cmap_category + ' colormaps', fontsize=14) for ax, name in zip(axes, cmap_list): ax.imshow(gradient, aspect='auto', cmap=plt.get_cmap(name)) 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: plot_color_gradients(cmap_category, cmap_list, nrows) 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) #add a length key 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) #Add titles!
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) #Add titles!
Populating the interactive namespace from numpy and matplotlib
<Container object of 3 artists>
In [ ]: