CoCalc Public FilesATMS-305-Content / Week-9 / Week 9 Exercise: March 13.ipynbOpen with one click!
Authors: Anna Nesbitt, Steve Nesbitt
Views : 23
Compute Environment: Ubuntu 18.04 (Deprecated)

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 [1]:
%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 [2]:
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:

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

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

In [5]:
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]
In [40]:
%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()

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

In [59]:
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 [62]:
%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()

Let's change the contour interval and the colormap!

In [65]:
%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 [72]:
%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 [73]:
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!

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

Time-height cross sections

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

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