# 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)