###EMANUEL 2005 FIGURE 1 #LANDSEA 2005 PDI DEFINITION (Caption to Figure 1) #PDI takes into account frequency, duration #and intensity of tropical cyclones by cubing the winds during the lifetime of the systems while they are of #at least tropical-storm force (18 m s-1) and summing them up for the year. Values shown are multiplied by #10-6 in units of m3 s-3. ##LOAD HURRICANE DATA cwd.data<-"http://www.climateaudit.org/data/hurricane/unisys" Track<-read.table(file.path(cwd.data,"Track.ATL.txt"),sep="\t",header=TRUE) dim(Track) # 38651 9 hurricane<-read.table(file.path(cwd.data,"hurricane.ATL.txt"),sep="\t",header=TRUE) dim(hurricane) # 1362 11 year0<-range(hurricane$year);year0<-year0[1]:year0[2];K<-length(year0) #METRIC CONVERSION #conversion from knots/hour to m sec-1 #1 knot = 1.15 mph #1 mi = 1609 m #1 knot = 1.15* 1609/3600 = 0.5139861 m sec-1 Track$wind.ms<- 0.5139861 * Track$wind # m sec-1 #BIN-PIN FUNCTION f<-function(x) { f<-filter(x,c(.25,.5,.25));f[is.na(f)]<-x[is.na(f)];f} #Emanuel function g<-function(x) { g<-filter(x,c(.25,.5,.25));g} #usual (1,2,1) smooth ##PDI USING ALL DATA IN METRIC AND QUARTER-DAYS hurricane$pdi<-tapply(Track$wind.ms^3,Track$id,sum,na.rm=T)/4 # divide by 4 since quarter-days chron.pdi<- c(tapply(hurricane$pdi,hurricane$year,sum))# range(chron.pdi*10^-6) # 0.04492386 8.05746141 range(4* chron.pdi*10^-6);# 0.1796954 32.2298456 #close to Landsea #therefore do without quarter-days hurricane$pdi<-tapply(Track$wind.ms^3,Track$id,sum,na.rm=T) chron.pdi<- c(tapply(hurricane$pdi,hurricane$year,sum))# #PDI according to Landsea description temp<-(Track$wind.ms>18)&!is.na(Track$wind.ms);sum(temp) # 27851 pdi.test<-tapply(Track$wind.ms[temp]^3,Track$id[temp],sum) hurricane$pdi.landsea<-0 hurricane$pdi.landsea[as.numeric(names(pdi.test))]<-pdi.test chron.pdi.landsea<-tapply(hurricane$pdi.landsea,hurricane$year,sum,na.rm=TRUE); range(chron.pdi.landsea*10^-6) # 0.1279780 31.8757674 #this more or less matches the scale #adjustments not considered so far #the wind-minimum limitation makes little difference to calculation #COMPARE TWO INDICES cor(chron.pdi,chron.pdi.landsea) #[1] 0.9995951 ##LANDSEA FIGURE 1 layout(1) par(mar=c(3,3,1,1)) plot(1949:2004,chron.pdi.landsea[(1949:2004)-1850]*10^-6,col="red",type="h",lwd=6,xlab="",ylab="",axes=FALSE,ylim=c(0,32),xlim=c(1947,2007)) axis(side=1);axis(side=2,las=1);box() lines(1949:2004,g(g(chron.pdi.landsea[(1949:2004)-1850]*10^-6)) ,lwd=2) lines(2005:2006,chron.pdi.landsea[(2005:2006)-1850]*10^-6,col="cyan",type="h",lwd=6) ##EMANUEL FIGURE 1 SCALE #supposedly multiplied by 2.1*10^-12 (!?! range(chron.pdi.landsea*2.1*10^-12)# 2.687539e-07 6.693911e-05 #this doesn't match the scale at all ##PLOT INCREMENTS AS ON BLOG y<-chron.pdi.landsea*10^-6 nf<-layout(array(1:3,dim=c(1,3)),widths=c(1.2,1.1,1.1)) par(mar=c(3,3,1,1)) plot(1949:2004,y[(1949:2004)-1850],type="h",col="grey80",lwd=5,xlim=c(1945,2006),ylab="",las=1,ylim=c(0,32)) lines(1949:2004,f(f(y[1:(K-2)]))[(1949:2004)-1850],type="l",lwd=2,lty=2) lines(1949:2004,g(g(y[1:(K-2)]))[(1949:2004)-1850],lwd=2,lty=4,col="red") year0<-1949:2004;fm<-lm(y[(1949:2004)-1850]~year0);lines(1949:2004,fm$fitted.values,lwd=2) par(mar=c(3,1,1,1)) #extend to 2005 plot(1949:2005,y[(1949:2005)-1850],type="h",col="grey80",lwd=5,xlim=c(1945,2006),ylab="",las=1,axes=FALSE,ylim=c(0,32)) axis(side=2,labels=FALSE);axis(side=1);box() lines(1949:2005,f(f(y[1:(K-1)]))[(1949:2005)-1850],type="l",lwd=2,lty=2) lines(1949:2005,g(g(y[1:(K-1)]))[(1949:2005)-1850],lwd=2,lty=3,col="red") year0<-1949:2005;fm<-lm(y[(1949:2005)-1850]~year0);lines(1949:2005,fm$fitted.values,lwd=2) par(mar=c(3,1,1,1)) #extend to 2006 plot(1949:2006,y[(1949:2006)-1850],type="h",col="grey80",lwd=5,xlim=c(1945,2006),ylab="",las=1,axes=FALSE,ylim=c(0,32)) axis(side=2,labels=FALSE);axis(side=1);box() lines(2006:2006,y[2006-1850],col="cyan",type="h",lwd=5) lines(1949:2006,g(g(y[1:K]))[(1949:2006)-1850],lwd=2,lty=3,col="red") lines(1949:2006,f(f(y[1:K]))[(1949:2006)-1850],lwd=2,lty=2) year0<-1949:2006; fm<-lm(y[(1949:2006)-1850]~year0);lines(1949:2006,fm$fitted.values,lwd=2)