Exploring Climate Data: Past and Future

Roland Viger, Rich Signell, USGS

First presented at the 2012 Unidata Workshop: Navigating Earth System Science Data, 9-13 July.

What if you were watching Ken Burns's Dustbowl, saw the striking image below, and wondered: "How much precipitation there really was back in the dustbowl years?" How easy is it to access and manipulate climate data in a scientific analysis? Here we'll show some powerful tools that make it easy.

In [1]:
from IPython.core.display import Image
Image('http://www-tc.pbs.org/kenburns/dustbowl/media/photos/s2571-lg.jpg')
Out[1]:

Above:Dust storm hits Hooker, OK, June 4, 1937.

To find out how much rainfall was there during the dust bowl years, we can use the USGS/CIDA GeoDataPortal (GDP) which can compute statistics of a gridded field within specified shapes, such as county outlines. Hooker is in Texas County, Oklahoma, so here we use the GDP to compute a historical time series of mean precipitation in Texas County using the PRISM dataset. We then compare to climate forecast projections to see if similar droughts are predicted to occur in the future, and what the impact of different climate scenarios might be.

In [2]:
import numpy as np
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import urllib
import os
from IPython.core.display import HTML
In [3]:
cd /usgs/data1/notebook
/usgs/data1/notebook
In [4]:
import pyGDP
import numpy as np
import matplotlib.dates as mdates
In [5]:
# Now, if we used GDP on the browser, we would navigate through something liks this:
HTML('<iframe src=http://screencast.com/t/K7KTcaFrSUc width=800 height=600></iframe>')
Out[5]:
In [6]:
# Create a pyGDP object
myGDP = pyGDP.pyGDPwebProcessing()
In [7]:
# Let's see what shapefiles are already available on the GDP server
# this changes with time, since uploaded shapefiles are kept for a few days
shapefiles = myGDP.getShapefiles()
print 'Available Shapefiles:'
for s in shapefiles:
    print s
Available Shapefiles:
draw:mexico
upload:lochvale_basin
draw:Chico
sample:simplified_HUC8s
sample:Alaska
sample:CSC_Boundaries
upload:OKCNTYD
upload:rgn3
upload:sleepers_basin
sample:CONUS_States
draw:test_daymet
upload:troutlake_basin
sample:FWS_LCC
upload:panola_basin
upload:NACP
In [8]:
# Is our shapefile there already?
# If not, upload it. 
OKshapeFile = 'upload:OKCNTYD'
if not OKshapeFile in shapefiles:
    shpfile = myGDP.uploadShapeFile('OKCNTYD.zip')
In [9]:
# Let's check the attributes of the shapefile
attributes = myGDP.getAttributes(OKshapeFile)
print "Shapefile attributes:"
for a in attributes:
    print a
Shapefile attributes:
OBJECTID_1
OBJECTID
RECTYPE
VERSION
REVISION
MODDATE
POLYID
FEATUREID
CNTRLONG
CNTRLAT
DESCRIP
STFIPS
Shape_area
Shape_len
In [10]:
# In this particular example, we are interested in attribute = 'DESCRIP', 
# which provides the County names for Oklahoma
user_attribute = 'DESCRIP'
values = myGDP.getValues(OKshapeFile, user_attribute)
print "Shapefile attribute values:"
for v in values:
    print v
Shapefile attribute values:
Adair
Alfalfa
Atoka
Beaver
Beckham
Blaine
Bryan
Caddo
Canadian
Carter
Cherokee
Choctaw
Cimarron
Cleveland
Coal
Comanche
Cotton
Craig
Creek
Custer
Delaware
Dewey
Ellis
Garfield
Garvin
Grady
Grant
Greer
Harmon
Harper
Haskell
Hughes
Jackson
Jefferson
Johnston
Kay
Kingfisher
Kiowa
Latimer
Le Flore
Lincoln
Logan
Love
Major
Marshall
Mayes
McClain
McCurtain
McIntosh
Murray
Muskogee
Noble
Nowata
Okfuskee
Oklahoma
Okmulgee
Osage
Ottawa
Pawnee
Payne
Pittsburg
Pontotoc
Pottawatomie
Pushmataha
Roger Mills
Rogers
Seminole
Sequoyah
Stephens
Texas
Tillman
Tulsa
Wagoner
Washington
Washita
Woods
Woodward
In [11]:
# we want Texas County, Oklahoma, which is where Hooker is located
user_value = 'Texas'

#Once we have our shapefile, attribute, and value, we pick a dataset we are interested in.f
# Let's see what datasets are available. 
dataSets = myGDP.getDataSetURI()
print "Available gridded datasets"
for d in dataSets:
    print d
The dataSetURI outputs a select few URIs and may not work with the specific shapefile you are providing.
To ensure compatibility, we recommend selecting a dataSetURI that is specific to the shapefile.
Or you may utilize the web gdp @ http://cida.usgs.gov/gdp/ to get a dataSet matching your specified shapefile.

Available gridded datasets
http://regclim.coas.oregonstate.edu:8080/thredds/dodsC/regcmdata/NCEP/merged/monthly/RegCM3_A2_monthly_merged_NCEP.ncml
dods://igsarm-cida-thredds1.er.usgs.gov:8080/thredds/dodsC/dcp/conus_grid.w_meta.ncml
http://cida.usgs.gov/qa/thredds/dodsC/prism
dods://igsarm-cida-thredds1.er.usgs.gov:8080/thredds/dodsC/maurer/maurer_brekke_w_meta.ncml
dods://igsarm-cida-thredds1.er.usgs.gov:8080/thredds/dodsC/dcp/alaska_grid.w_meta.ncml
dods://igsarm-cida-thredds1.er.usgs.gov:8080/thredds/dodsC/gmo/GMO_w_meta.ncml
In [12]:
# If you choose a DAP URL, use the "dods:" prefix, even
# if the list above has a "http:" prefix.
# For example:  dods://cida.usgs.gov/qa/thredds/dodsC/prism
# Let's see what data variables are in our dataset
dataSetURI = 'dods://cida.usgs.gov/qa/thredds/dodsC/prism'
dataTypes = myGDP.getDataType(dataSetURI)
print "Available variables:"
for d in dataTypes:
    print d
Available variables:
ppt
tmx
tmn
In [13]:
# Let's see what the available time range is for our data variable
user_dataType = 'ppt'  # precip
timeRange = myGDP.getTimeRange(dataSetURI, user_dataType)
for t in timeRange:
    print t
1895-01-01T00:00:00Z
2012-08-01T00:00:00Z
In [14]:
timeBegin = '1900-01-01T00:00:00Z'
timeEnd   = '2012-08-01T00:00:00Z'
In [15]:
# Once we have our shapefile, attribute, value, dataset, datatype, and timerange as inputs, we can go ahead
# and submit our request.
name1='gdp_texas_county_prism.csv'
if not os.path.exists(name1):
    url_csv = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, dataSetURI, user_dataType, 
          timeBegin, timeEnd, user_attribute, user_value, delim='COMMA', stat='MEAN' )
    f = urllib.urlretrieve(url_csv,name1)
In [16]:
# load historical PRISM precip
jd,precip=np.loadtxt(name1,unpack=True,skiprows=3,delimiter=',', 
                     converters={0: mdates.strpdate2num('%Y-%m-%dT%H:%M:%SZ')})
In [17]:
def boxfilt(data,boxwidth):
    from scipy import signal
    import numpy as np
    weights=signal.get_window('boxcar',boxwidth)
    dataf=np.convolve(data,weights/boxwidth,mode='same')
    dataf=np.ma.array(dataf)
    dataf[:boxwidth/2]=np.nan
    dataf[-boxwidth/2:]=np.nan
    dataf=np.ma.masked_where(dataf==np.nan,dataf)
    return dataf
In [18]:
# PRISM data is monthly:  filter over 36 months
plp=boxfilt(precip,36)

fig=plt.figure(figsize=(12,2), dpi=80) 
ax1 = fig.add_subplot(111)
g1=ax1.plot_date(jd,plp,fmt='b-')
g2=ax1.plot_date(jd,0*jd+np.mean(precip),fmt='k-')
fig.autofmt_xdate()
plt.title('Average Precip for Texas County, Oklahoma, calculated via GDP using PRISM data ')
plt.grid()
In [19]:
HTML('<iframe src=http://www.ipcc.ch/publications_and_data/ar4/wg1/en/spmsspm-projections-of.html width=900 height=350></iframe>')
Out[19]:
In [20]:
hayhoe_URI ='dods://cida-eros-thredds1.er.usgs.gov:8082/thredds/dodsC/dcp/conus_grid.w_meta.ncml'
dataType = 'ccsm3_a1fi_pr'
timeRange = myGDP.getTimeRange(hayhoe_URI, dataType)
In [21]:
timeRange
Out[21]:
['1960-01-01T00:00:00Z', '2099-12-31T00:00:00Z']
In [22]:
# retrieve the CCSM3 model A1FI "Business-as-Usual" scenario:
name2='gdp_texas_county_ccsm_a1f1.csv'
if not os.path.exists(name1):
    result2 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, hayhoe_URI, dataType,
            timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )
    f = urllib.urlretrieve(result2,name2)                                       
In [23]:
# now retrieve the CCSM3 model B1 "Eco-Friendly" scenario:
name3='gdp_texas_county_ccsm_b1.csv'
if not os.path.exists(name1):
    dataType = 'ccsm3_b1_pr'
    result3 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, hayhoe_URI, dataType,
            timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )
    f = urllib.urlretrieve(result3,name3)
In [24]:
# Load the GDP result for: CCSM A1FI "Business-as-Usual" scenario:
jd_a1f1,precip_a1f1 = np.loadtxt(name2,unpack=True,skiprows=3,
    delimiter=',',converters={0: mdates.strpdate2num('%Y-%m-%dT%H:%M:%SZ')}) 

# Load the GDP result for:  CCSM B1 "Eco-Friendly" scenario:
jd_b1,precip_b1     = np.loadtxt(name3,unpack=True,skiprows=3,
    delimiter=',',converters={0: mdates.strpdate2num('%Y-%m-%dT%H:%M:%SZ')}) 
In [25]:
# Hayhoe climate downscaling is hourly: filter over 1080 days (36 months)
plp_a1f1=boxfilt(precip_a1f1,1080)
plp_b1=boxfilt(precip_b1,1080)
#plp_a1b_c=boxfilt(precip_a1b_c,36)
In [26]:
fig=plt.figure(figsize=(15,3), dpi=80) 
ax1 = fig.add_subplot(111)
fac=30. # convert from mm/day to mm/month (approx)
# plot A1F1 scenario
g1=ax1.plot_date(jd_a1f1,plp_a1f1*fac,fmt='b-')
# plot B1 scenario 
g2=ax1.plot_date(jd_b1,plp_b1*fac,fmt='g-')
# plot PRISM data
g3=ax1.plot_date(jd,plp,fmt='r-')  # for some reason when I add this the labels get borked
ax1.xaxis.set_major_locator(mdates.YearLocator(10,month=1,day=1))
ylabel('mm/month')
plt.title('Average Precip for Texas County, Oklahoma, calculated via GDP using Hayhoe Downscaled GCM ')
grid()
legend(('A1FI','B1','PRISM Data'),loc='upper left')
Out[26]:
<matplotlib.legend.Legend at 0x457ab10>

As we can see from the above plot, the CCSM model is not doing very well simulating the precipitation in Texas County, OK during the period when the simulation and data overlap (1960-present). This makes us less confident about the future precipitation simulations, and suggests we might need to try some different climate models and learn a bit more about climate simulations. When we do learn more, we find out that models have known biases in certain regions.

Now just to show that we can access more than climate model time series, let's extract precipitation data from a dry winter (1936-1937) and a normal winter (2009-2010) for Texas County and look at the spatial patterns.

We'll use the netCDF4-Python library, which allows us to open OPeNDAP datasets just as if they were local NetCDF files.

In [27]:
import netCDF4
url='http://cida.usgs.gov/qa/thredds/dodsC/prism'
box = [-102,36.5,-101,37]  # Bounding box for Texas County, Oklahoma
In [28]:
# define a mean precipitation function, here hard-wired for the PRISM data
def mean_precip(nc,bbox=None,start=None,stop=None):
    lon=nc.variables['lon'][:]
    lat=nc.variables['lat'][:]
    tindex0=netCDF4.date2index(start,nc.variables['time'],select='nearest')
    tindex1=netCDF4.date2index(stop,nc.variables['time'],select='nearest')
    bi=(lon>=box[0])&(lon<=box[2])
    bj=(lat>=box[1])&(lat<=box[3])
    p=nc.variables['ppt'][tindex0:tindex1,bj,bi]
    latmin=np.min(lat[bj])
    p=np.mean(p,axis=0)
    lon=lon[bi]
    lat=lat[bj]
    return p,lon,lat
In [29]:
nc = netCDF4.Dataset(url)
p,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(1936,11,1,0,0),
                        stop=datetime.datetime(1937,4,1,0,0))
p2,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(2009,11,1,0,0),
                       stop=datetime.datetime(2010,4,1,0,0))
latmin = np.min(lat)
In [30]:
fig = plt.figure(figsize=(12,5), dpi=80) 
ax = fig.add_axes([0.1, 0.15, 0.3, 0.8])
pc = ax.pcolormesh(lon, lat, p, cmap=plt.cm.jet_r)
ax.set_aspect(1.0/np.cos(latmin * np.pi / 180.0))
plt.title('Precip in Texas County, Oklahoma: Winter 1936-1937')

cbax = fig.add_axes([0.45, 0.3, 0.03, 0.4])
cb = plt.colorbar(pc, cax=cbax,  orientation='vertical')
cb.set_label('Precip (mm/month)')

ax2 = fig.add_axes([0.6, 0.15, 0.3, 0.8])
pc2 = ax2.pcolormesh(lon, lat, p2, cmap=plt.cm.jet_r)
ax2.set_aspect(1.0/np.cos(latmin * np.pi / 180.0))
plt.title('Precip in Texas County, Oklahoma: Winter 2009-2010')

cbax2 = fig.add_axes([0.95, 0.3, 0.03, 0.4])
cb2 = plt.colorbar(pc2, cax=cbax2,  orientation='vertical')
cb2.set_label('Precip (mm/month)')


plt.show()

From the above patterns, we can see that it's significantly drier in the northwestern part of the county in both years. We can also see that the maximum precip in 1936-1937 is less than the minimum precipitation in 2009-2010. We can see just how much each part of the county was drier by doing the different plot below.

In [31]:
fig=plt.figure(figsize=(12,5), dpi=80) 
ax3 = fig.add_axes([0.1, 0.15, 0.3, 0.8])
pc3 = ax3.pcolormesh(lon, lat, p2-p, cmap=plt.cm.jet_r)
ax3.set_aspect(1.0/np.cos(latmin * np.pi / 180.0))
plt.title('Precip in Texas County, Oklahoma: Difference 2010-1937')

cbax3 = fig.add_axes([0.45, 0.3, 0.03, 0.4])
cb3 = plt.colorbar(pc3, cax=cbax3,  orientation='vertical')
cb3.set_label('Precip (mm/month)')

The above plot shows that relative to 2010, the drought during 1937 had the biggest different in precip in the northeastern part of the county.

Hopefully this demo inspires other investigation of historical and projected climate data using the GDP and Python.