First presented at the 2012 Unidata Workshop: Navigating Earth System Science Data, 9-13 July.
What if you were watching Ken Burns's "The Dust Bowl", 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.
from IPython.core.display import Image
Image('http://www-tc.pbs.org/kenburns/dustbowl/media/photos/s2571-lg.jpg')
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.
import numpy as np
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import urllib
import os
from IPython.core.display import HTML
import time
import pandas as pd
cd /usgs/data2/notebook
/usgs/data2/notebook
import pyGDP
import numpy as np
import matplotlib.dates as mdates
One way to interface with the GDP is with the interactive web interface, shown below. In this interface, you can upload a shapefile or draw on the screen to define a polygon region, then you specify the statistics and datasets you want to use via dropdown menus.
HTML('<iframe src=http://screencast.com/t/K7KTcaFrSUc width=800 height=600></iframe>')
Here we use the python interface to the GDP, called PyGDP, which allows for scripting. You can get the code and documentation at https://github.com/USGS-CIDA/pyGDP.
# Create a pyGDP object
myGDP = pyGDP.pyGDPwebProcessing()
# 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: derivative:CONUS_States draw:LA_Extended sample:simplified_HUC8s derivative:Level_III_Ecoregions draw:test draw:Siso derivative:FWS_LCC draw:KS sample:Alaska sample:CSC_Boundaries upload:OKCNTYD upload:Wllstn_grid derivative:US_Counties sample:CONUS_states sample:nps_boundary_2013 sample:FWS_LCC derivative:NCA_Regions derivative:wbdhu8_alb_simp upload:WllstnGrid draw:tests3
# Is our shapefile there already?
# If not, upload it.
OKshapeFile = 'upload:OKCNTYD'
if not OKshapeFile in shapefiles:
shpfile = myGDP.uploadShapeFile('OKCNTYD.zip')
# 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 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
# we want Texas County, Oklahoma, which is where Hooker is located
user_value = 'Texas'
# Let's see what gridded datasets are available for the GDP to operate on
dataSets = myGDP.getDataSetURI()
print "Available gridded datasets:"
for d in dataSets:
print d[0]
Available gridded datasets: title Conterminous U.S. actual evapotranspiration data Great Lakes Coastal Forecasting System Nowcast/Lake Michigan Nowcast History 2D Gridded Observed Meteorological Data, 1950-1999 Urban Growth Projection for Southeast Regional Assessment Project Great Lakes Coastal Forecasting System Nowcast/Lake Ontario Nowcast History 2D Great Lakes Coastal Forecasting System Nowcast/Lake Huron Nowcast History 2D Global Land Data Assimilation System ICLUS v1.3 Housing Density Classification Projections for the Conterminous USA Sea Level Rise Projections for DSL-SAMBI ICLUS v1.3 Raw Housing Density and Impervious Surface Projections for the Conterminous USA Thermosteric Sea Level Rise Projections with Parametric Uncertainty Great Lakes Coastal Forecasting System Nowcast/Lake Superior Nowcast History 2D Urban Growth Projection for DSL-SAMBI Great Lakes Coastal Forecasting System Nowcast/Lake Erie Nowcast History 2D CIG Northern US Rockies and Pacific Northwest Statistical Downscaling Daily Statistically Downscaled Climate Projections for the US and southern Canada east of the Rocky Mountains. Gridded Observed Meteorological Data: 1949-2010 Half degree-Alaska Daily Downscaled Climate Projections by Katharine Hayhoe Eighth degree-CONUS Daily Downscaled Climate Projections by Katharine Hayhoe Mean Monthly Evaporation Atlas for the Contiguous 48 United States (1956-1970) Parameter-elevation Regressions on Independent Slopes Model Monthly Climate Data for the Continental United States. Real Time River Forecasting Center Quantitative Precipitation Estimate River Forecasting Center Quantitative Precipitation Estimate Archive University of Idaho Daily Meteorological data for continental US Bias Corrected Constructed Analogs V2 Daily Climate Projections Bias Corrected Spatially Downscaled Monthly Climate Predictions USGS Dynamical Downscaled Regional Climate 800m Downscaled NEX CMIP5 Climate Projections for the Continental US Normal Incident Solar Radiation Atlas North American Land Data Assimilation System Phase 2 Daymet Daily surface weather on a 1km grid for North America, 1980-2012 Columbia River Basin Daily MACA-VIC Results Western US Hydroclimate Scenarios Project Observations and Statistically Downscaled Data TopoWx: Topoclimatic Daily Air Temperature Dataset for the Conterminous United States Western US Hydroclimate Scenarios Project Dynamically Downscaled Data Bias Corrected Spatially Downscaled Monthly CMIP5 Climate Projections California Basin Characterization Model Downscaled Climate and Hydrology United States Quantitative Precipitation Archive
dataSets[0][0]
'title'
df = pd.DataFrame(dataSets[1:],columns=['title','abstract','urls'])
df.head()
title | abstract | urls | |
---|---|---|---|
0 | Conterminous U.S. actual evapotranspiration data | Actual ET (ETa) is produced using the operatio... | [dods://cida.usgs.gov/thredds/dodsC/ssebopeta/... |
1 | Great Lakes Coastal Forecasting System Nowcast... | Great Lakes Coastal Forecasting System Nowcast... | [dods://michigan.glin.net:8080/thredds/dodsC/g... |
2 | Gridded Observed Meteorological Data, 1950-1999 | These daily gridded observations at 1/8 degre... | [dods://cida.usgs.gov/thredds/dodsC/gmo/GMO_w_... |
3 | Urban Growth Projection for Southeast Regional... | This dataset represents the extent of urbaniza... | [dods://cida.usgs.gov/thredds/dodsC/serap/sera... |
4 | Great Lakes Coastal Forecasting System Nowcast... | Great Lakes Coastal Forecasting System Nowcast... | [dods://michigan.glin.net:8080/thredds/dodsC/g... |
print df['title']
0 Conterminous U.S. actual evapotranspiration data 1 Great Lakes Coastal Forecasting System Nowcast... 2 Gridded Observed Meteorological Data, 1950-1999 3 Urban Growth Projection for Southeast Regional... 4 Great Lakes Coastal Forecasting System Nowcast... 5 Great Lakes Coastal Forecasting System Nowcast... 6 Global Land Data Assimilation System 7 ICLUS v1.3 Housing Density Classification Proj... 8 Sea Level Rise Projections for DSL-SAMBI 9 ICLUS v1.3 Raw Housing Density and Impervious ... 10 Thermosteric Sea Level Rise Projections with P... 11 Great Lakes Coastal Forecasting System Nowcast... 12 Urban Growth Projection for DSL-SAMBI 13 Great Lakes Coastal Forecasting System Nowcast... 14 CIG Northern US Rockies and Pacific Northwest ... 15 Daily Statistically Downscaled Climate Project... 16 Gridded Observed Meteorological Data: 1949-2010 17 Half degree-Alaska Daily Downscaled Climate Pr... 18 Eighth degree-CONUS Daily Downscaled Climate P... 19 Mean Monthly Evaporation Atlas for the Contigu... 20 Parameter-elevation Regressions on Independent... 21 Real Time River Forecasting Center Quantitativ... 22 River Forecasting Center Quantitative Precipit... 23 University of Idaho Daily Meteorological data ... 24 Bias Corrected Constructed Analogs V2 Daily Cl... 25 Bias Corrected Spatially Downscaled Monthly Cl... 26 USGS Dynamical Downscaled Regional Climate 27 800m Downscaled NEX CMIP5 Climate Projections ... 28 Normal Incident Solar Radiation Atlas 29 North American Land Data Assimilation System P... 30 Daymet Daily surface weather on a 1km grid for... 31 Columbia River Basin Daily MACA-VIC Results 32 Western US Hydroclimate Scenarios Project Obse... 33 TopoWx: Topoclimatic Daily Air Temperature Dat... 34 Western US Hydroclimate Scenarios Project Dyna... 35 Bias Corrected Spatially Downscaled Monthly CM... 36 California Basin Characterization Model Downsc... 37 United States Quantitative Precipitation Archive Name: title, dtype: object
df.ix[20].urls
['dods://cida.usgs.gov/thredds/dodsC/prism']
# 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/thredds/dodsC/prism'
dataTypes = myGDP.getDataType(dataSetURI)
print "Available variables:"
for d in dataTypes:
print d
Available variables: ppt tmx tmn
# 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 2013-02-01T00:00:00Z
timeBegin = '1900-01-01T00:00:00Z'
timeEnd = '2012-08-01T00:00:00Z'
# 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)
# load historical PRISM precip
jd,precip=np.loadtxt(name1,unpack=True,skiprows=3,delimiter=',',
converters={0: mdates.strpdate2num('%Y-%m-%dT%H:%M:%SZ')})
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
# 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()
HTML('<iframe src=http://www.ipcc.ch/publications_and_data/ar4/wg1/en/spmsspm-projections-of.html width=900 height=350></iframe>')
#hayhoe_URI ='dods://cida-eros-thredds1.er.usgs.gov:8082/thredds/dodsC/dcp/conus_grid.w_meta.ncml'
hayhoe_URI ='dods://cida.usgs.gov/thredds/dodsC/dcp/conus'
timeRange = myGDP.getTimeRange(hayhoe_URI, dataType)
timeRange
[]
# retrieve the CCSM3 model A1FI "Business-as-Usual" scenario:
name2='gdp_texas_county_ccsm_a1fi.csv'
if not os.path.exists(name2):
dataType = 'ccsm-a1fi-pr-NAm-grid'
result2 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, hayhoe_URI, dataType,
timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )
f = urllib.urlretrieve(result2,name2)
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-163-ff87e1a31549> in <module>() 4 dataType = 'ccsm-a1fi-pr-NAm-grid' 5 result2 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, hayhoe_URI, dataType, ----> 6 timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' ) 7 f = urllib.urlretrieve(result2,name2) IndexError: list index out of range
# now retrieve the CCSM3 model B1 "Eco-Friendly" scenario:
time0=time.time();
name3='gdp_texas_county_ccsm_b1.csv'
if not os.path.exists(name3):
dataType = 'ccsm-b1-pr-NAm-grid'
result3 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, hayhoe_URI, dataType,
timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )
f = urllib.urlretrieve(result3,name3)
print('elapsed time=%d s' % (time.time()-time0))
# 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')})
# 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)
fig=plt.figure(figsize=(15,3), dpi=80)
ax1 = fig.add_subplot(111)
fac=30. # convert from mm/day to mm/month (approx)
# plot A1FI 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')
<matplotlib.legend.Legend at 0x7f808d1dcfd0>
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.
import netCDF4
url='http://cida.usgs.gov/thredds/dodsC/prism'
box = [-102,36.5,-100.95,37] # Bounding box for Texas County, Oklahoma
#box = [-104,36.,-100,39.0] # Bounding box for larger dust bowl region
# 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
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)
# look at March 1935, just before black sunday on April 14, 1935
nc = netCDF4.Dataset(url)
p,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(1935,3,1,0,0),
stop=datetime.datetime(1935,4,1,0,0))
p2,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(2009,3,1,0,0),
stop=datetime.datetime(2009,4,1,0,0))
latmin = np.min(lat)
import cartopy.crs as ccrs
import cartopy.feature as cfeature
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
ax = plt.axes(projection=ccrs.PlateCarree())
pc = ax.pcolormesh(lon, lat, p, cmap=plt.cm.jet_r)
ax.add_feature(states_provinces,edgecolor='gray')
ax.text(-101,36.86,'Hooker',align='left')
ax.plot(-101,36.86,'o')
cb = plt.colorbar(pc, cax=cbax, orientation='vertical')
cb.set_label('Precip (mm/month)')
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.
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.