# 2/7/2009 Nicolas Nierenberg
# Analysis code for ECMWF humidity data
# As much as I would like to have made this "turnkey" the data has to be requested
# from ECMWF using a form. The file "ECMWF/output.nc" is a result of this request for the ERA-40
# data for the years 1958-2001 at all pressure levels. This is the full set of complete
# years. Also for all grid cells.
library(ncdf)
# Following are some general routines that I wrote. Probably should put them in a separate
# file but I haven't sorted how to do that yet.
# This is some weird code to get a p-value for a specific variable in a regression.
# Best way I could figure out how to do it. I didn't use this below as I decided to
# just present the confidence intervals.
get.pvalue = function(model,variable) {
sum.model=summary(model)
cf=coef(sum.model)
cf[variable,"Pr(>|t|)"]
}
# This function probably already exists in R but I couldn't find it.
latToRads = function(lat) (abs(lat)) / (180/pi)
# This computes a weighted average of a vector of values given a vector of latitudes.
# The idea is that the area of a grid cell varies as the cosine of the latitude.
wAvg = function(x,lat){
sum(x*cos(latToRads(lat)))/sum(cos(latToRads(lat)))
}
# Computes a vector of weights, but I didn't use it here.
weights= function(lat){
cos(latToRads(lat))/sum(cos(latToRads(lat)))
}
# Following is the code that is actually focused on the analysis of the humidity data.
# open the ncdf file from ECMWF and read the specific humidity and relative humidity
handle=open.ncdf("ECMWF/output.nc")
ECMWF.q=get.var.ncdf(handle,varid="q")
ECMWF.r=get.var.ncdf(handle,varid="r")
latitude=get.var.ncdf(handle,varid="latitude")
q.avgs = function(year.range,lat,press){
q=numeric()
for(year in year.range) {
startmonth=((year-1)*12)+1
q[year]=sum(ECMWF.q[,lat,press,startmonth:startmonth+11])/12
}
q
}
# Compute annual weighted average values for q and r for a specific pressure, set of years, and set of lattitudes. Returns a data frame with years, average q, and average r.
qr.trend = function(press,year.range,lat.range){
lats = get.var.ncdf(handle,"latitude")
d = data.frame(year=year.range,avg.q=numeric(length(year.range)))
q = numeric(length(lats))
if(!is.na(ECMWF.r))
r = numeric(length(lats))
for(year in 1:length(year.range)){
startmonth=((year-1)*12)+1
for(lat in lat.range){
q[lat] = sum(ECMWF.q[,lat,press,startmonth:startmonth+11])/12
if(!is.na(ECMWF.r))
r[lat] = sum(ECMWF.r[,lat,press,startmonth:startmonth+11])/12
}
d$avg.q[year]=wAvg(q[lat.range],lats[lat.range])
if(!is.na(ECMWF.r))
d$avg.r[year]=wAvg(r[lat.range],lats[lat.range])
}
d
}
# Compute the trends and 95% confidence intervals for q and r for the years and latitudes specified. Returned
# as a data table. Uses qr.trend to compute the yearly averages. It uses the function lm to regress
# the values against the year and uses the year coefficient as the trend. The function confint uses
# the results of the regression to compute the 95% confidence interval.
qr.trends = function(year.range,lat.range) {
press=get.var.ncdf(handle,varid="levelist")
len=length(press)
t=data.frame(press=press)
for(p in 1:length(press)){
trend=qr.trend(p,year.range,lat.range)
res = lm(avg.q~year,data=trend)
t$trend.q[p]=res$coeff[2]
interval=confint(res,"year")
t$lo.q[p]=interval[1,1]
t$hi.q[p]=interval[1,2]
if(!is.na(ECMWF.r)){
res = lm(avg.r~year,data=trend)
t$trend.r[p]=res$coeff[2]
interval=confint(res,"year")
t$lo.r[p]=interval[1,1]
t$hi.r[p]=interval[1,2]
}
}
t
}
# global trends for the entire period.
options(digits=2,scipen=1)
qr.trends(1:44,1:73)
# Compute the NH mid lats 20N to 50N for the entire period
qr.trends(1:44,17:29)
# Compute the NF mid lats using 30N to 55N (Just a different definition of mid lats)
qr.trends(1:44,15:25)
# Compute the SH mid lats 20S to 50S for the entire period
qr.trends(1:44,45:57)
# Compute mid lats and tropics for the entire period
qr.trends(1:44,17:57)
# The following are for the satellite era of 1979 to 2001
# Compute the NH mid lats 20N to 50N for the entire period
qr.trends(22:44,17:29)
# Compute the NF mid lats using 30N to 55N (Just a different definition of mid lats)
qr.trends(22:44,15:25)
# Analysis of ERA Interim/By reading it into an object of the same name I can
# use the functions defined above. I just put the file in a different folder when
# I downloaded it from the ECMWF server. To save space I just analyzed q in this
# dataset.
handle=open.ncdf("ECMWF-ERA_Interim/output.nc")
ECMWF.q=get.var.ncdf(handle,varid="q")
ECMWF.r=NA
# Compute the NH Trend for the 19 year period.
qr.trends(1:19,28:48)
# Compute all mid lat trends for the 19 year period.
qr.trends(1:19,28:95)