#Notes on BEST article on Watts data load("d:/climate/data/station/watts/station.tab") #good bad grouping station$rating=NA station$rating[station$RATING==1|station$RATING==2]="good" station$rating[station$RATING==4|station$RATING==5]="bad" station$rating=factor(station$rating) station$grid=factor(station$GRID) ################################### ##RECOLLATE TRENDS source("http://www.climateaudit.info/scripts/utilities.txt") trend=function(x) lm(x~I((1:length(x))/120))$coef[2] source("d:/climate/scripts/make.table.txt") download.file("http://www.climateaudit.info/data/station/ushcn/ushcn.collation.tab","temp",mode="wb") load("temp") X=ushcn.collation$raw; id=dimnames(X)[[2]] #calculate anomalies 1961-1990 X=X[,paste(station$id)] X=window(X,start=1961) for(i in 1:ncol(X)) X[,i]=anom(X[,i]) N=ncol(X) #calculate trends for post-1950 x=apply( window(X,start=1950),2,trend) station$traw=station$trend=x x=apply( window(X,start=1979),2,trend) station$t79=station$trend=x ###END RECOLLATE TRENDS ############################## ########################### ## COUNTS AND INFORMATION # By 5 classes count=tapply(!is.na(station$RATING),station$RATING,sum) count["sum"]=sum(count) Qual=data.frame(count) Qual$rohde=c(15,73,216,627,78,1221-1009,1221) Qual$roh_trend=c(.391,.534,.243,.373,.51,NA,.32) Qual$trend=round(c(tapply(station$trend,station$RATING,mean,na.rm=T),mean(station$trend,na.rm=T)),2) Qual$t79=round(c(tapply(station$t79,station$RATING,mean,na.rm=T),mean(station$t79,na.rm=T)),2) row.names(Qual)=paste("#", row.names(Qual) ) Qual # count rohde roh_trend trend t79 # 1 52 15 0.391 0.15 0.26 # 2 133 73 0.534 0.15 0.20 # 3 339 216 0.243 0.15 0.26 # 4 375 627 0.373 0.15 0.28 # 5 131 78 0.510 0.18 0.24 # x 191 212 NA 0.10 0.22 # sum 1221 1221 0.320 0.14 0.25 make.table(Qual) # by good/bad tapply(!is.na(station$watts),station$watts,sum); #CRS HYGR MMTS OTHER ##409 49 751 12 tapply(!is.na(station$TYPE),station$TYPE,sum); # R S U #820 288 113 count=tapply(!is.na(station$rating),station$rating,sum) count["sum"]=sum(count) Qb=data.frame(count) Qb$rohde =c(705,88,705+88) Qb$roh_trend=c(.388,.509,NA) Qb$trend=c(round(tapply(station$trend,station$rating,mean,na.rm=T),2),NA) Qb$t79=round(c(tapply(station$t79,station$rating,mean,na.rm=T),mean(station$t79,na.rm=T)),2) row.names(Qb)=paste("#", row.names(Qb) ) temp=!is.na(station$rating)&station$TYPE=="R" Qb$rural=c(tapply(!is.na(station$TYPE[temp]),station$rating[temp],sum),NA) temp=!is.na(station$rating)&station$TYPE=="U" Qb$urban=c(tapply(!is.na(station$TYPE[temp]),station$rating[temp],sum),NA) temp=!is.na(station$rating)&station$watts=="CRS" Qb$CRS=c(tapply(!is.na(station$TYPE[temp]),station$rating[temp],sum),NA) temp=!is.na(station$rating)&station$watts=="MMTS" Qb$MMTS=c(tapply(!is.na(station$TYPE[temp]),station$rating[temp],sum),NA) Qb[3,5:8]=apply(Qb[1:2,5:8],2,sum) Qb # count rohde roh_trend trend rural urban CRS MMTS t79 # bad 506 705 0.388 0.16 337 38 146 344 0.27 # good 185 88 0.509 0.15 119 28 78 79 0.22 # sum 691 793 NA NA 456 66 224 423 0.25 ###################### ##stratify trends on 1030 classified stations library(nlme) #R vs Urban station$trend=station$t79 Work=groupedData(trend~id|TYPE/rating/grid,station[!is.na(station$trend)&!station$RATING=="x",c("trend","id","grid","rating","TYPE","lat","long")]) fm=lmList(trend~1|TYPE/rating/grid,Work[!is.na(Work$rating),]) #691 Y=ranef(fm) #164 combos fixef(fm) # 0.1690481 Y=data.frame(work=row.names(Y),trend=Y[,1]) x=Y$work; x=strsplit(gsub("/","%",x),"%") Y$rural=sapply(x, function(A) A[[1]]) Y$good=sapply(x, function(A) A[[2]]) Y$grid=sapply(x, function(A) A[[3]]) Y=Y[,c(3,4,5,2)] #Y=Y[!Y$rural=="S",] Y$count=factor(Y$grid) levels(Y$count)=tapply(!is.na(Y$grid),Y$grid,sum) Y$combo=paste(Y$rural,Y$good,sep="_") tapply(!is.na(Y$combo),Y$combo,sum) #R_bad R_good U_bad U_good # 41 35 17 19 round(fixef(fm)+tapply(Y$trend,Y$combo,mean),2) # R_bad R_good S_bad S_good U_bad U_good #0.25 0.17 0.31 0.24 0.30 0.31 #1950 on # R_bad R_good S_bad S_good U_bad U_good # 0.14 0.12 0.20 0.20 0.16 0.24 ##INSTR:MMTS aND crs temp= !station$watts=="HYGR" &!station$watts=="OTHER" Work=groupedData(trend~id|watts/rating/grid,station[!is.na(station$trend)&!station$RATING=="x"&temp,c("trend","id","grid","rating","TYPE","watts")]) Work$combo=paste(Work$watts,Work$TYPE,sep="_") fm=lmList(trend~1|watts/rating/grid,Work[!is.na(Work$rating),]) fixef(fm) # 0.1666171 Y=ranef(fm) Y=data.frame(work=row.names(Y),trend=Y[,1]) x=Y$work; x=strsplit(gsub("/","%",x),"%") Y$watts=sapply(x, function(A) A[[1]]) Y$good=sapply(x, function(A) A[[2]]) Y$grid=sapply(x, function(A) A[[3]]) Y=Y[,c(3,4,5,2)] Y$combo=paste(Y$watts,Y$good,sep="_") tapply(!is.na(Y$combo),Y$combo,sum) # CRS_bad CRS_good MMTS_bad MMTS_good # 33 34 40 29 round(fixef(fm)+ tapply(Y$trend,Y$combo,mean),2) #1950 on # CRS_bad CRS_good MMTS_bad MMTS_good # 0.18 0.21 0.17 0.09 #1979 #CRS_bad CRS_good MMTS_bad MMTS_good # 0.30 0.25 0.29 0.14 #BOTH fm=lmList(trend~1|combo/rating/grid,Work[!is.na(Work$rating),]) fixef(fm) Y=ranef(fm) Y=data.frame(work=row.names(Y),trend=Y[,1]) x=Y$work; x=strsplit(gsub("/","%",x),"%") Y$watts=sapply(x, function(A) A[[1]]) Y$good=sapply(x, function(A) A[[2]]) Y$grid=sapply(x, function(A) A[[3]]) Y=Y[,c(3,4,5,2)] Y$combo=paste(Y$watts,Y$good,sep="_") tapply(!is.na(Y$combo),Y$combo,sum) # CRS_R_bad CRS_R_good CRS_S_bad CRS_S_good CRS_U_bad CRS_U_good MMTS_R_bad MMTS_R_good MMTS_S_bad # 31 26 16 13 8 10 38 28 28 #MMTS_S_good MMTS_U_bad MMTS_U_good # 9 12 1 round(fixef(fm)+tapply(Y$trend,Y$combo,mean),2) #1950 # CRS_R_bad CRS_R_good CRS_S_bad CRS_S_good CRS_U_bad CRS_U_good MMTS_R_bad MMTS_R_good MMTS_S_bad # 0.16 0.19 0.21 0.18 0.34 0.25 0.14 0.07 0.20 #MMTS_S_good MMTS_U_bad MMTS_U_good # 0.11 0.14 0.35 #1979 # CRS_R_bad CRS_R_good CRS_S_bad CRS_S_good CRS_U_bad CRS_U_good MMTS_R_bad MMTS_R_good MMTS_S_bad # 0.26 0.25 0.35 0.17 0.42 0.27 0.25 0.11 0.33 #MMTS_S_good MMTS_U_bad MMTS_U_good #Qb 0.23 0.28 0.58