# 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
draw:Chico
sample:simplified_HUC8s
sample:CSC_Boundaries
sample:CONUS_States
draw:test_daymet
sample:FWS_LCC

In [8]:
# Is our shapefile there already?
if not OKshapeFile in shapefiles:

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:
Alfalfa
Atoka
Beaver
Beckham
Blaine
Bryan
Carter
Cherokee
Choctaw
Cimarron
Cleveland
Coal
Comanche
Cotton
Craig
Creek
Custer
Delaware
Dewey
Ellis
Garfield
Garvin
Grant
Greer
Harmon
Harper
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/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
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
return dataf

In [18]:
# PRISM data is monthly:  filter over 36 months
plp=boxfilt(precip,36)

fig=plt.figure(figsize=(12,2), dpi=80)
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:
delimiter=',',converters={0: mdates.strpdate2num('%Y-%m-%dT%H:%M:%SZ')})

# Load the GDP result for:  CCSM B1 "Eco-Friendly" scenario:
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)
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.