##COMPARE NOAA TO GRIDDED ################### ##README ############### #ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/grid/README_GRID_TEMP.txt #GHCN homogeneity adjusted data was the primary source for developing the gridded #fields. In grid boxes without homogeneity adjusted data, GHCN raw data was used to provide #additional coverage when possible. Each month of data consists of 2592 gridded data points #produced on a 5 X 5 degree basis for the entire globe (72 longitude X 36 latitude grid boxes). ##NEWS #http://www.noaanews.noaa.gov/stories2007/s2819.htm #http://www.ncdc.noaa.gov/oa/climate/research/2007/feb/feb07.html ######################## ###MAKE GRIDDED FILE FROM MANUALLY UNZIPPED FILE ############################## ##download and unzip ## ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/anom/anom-grid2-1880-current.dat.gz ## to ## d:/climate/data/noaa/anom-grid2-1880-current.dat url<-file.path("d:/climate/data/noaa","anom-grid2-1880-current.dat") fred<-readLines(url) #startsin 1880 Jan N<-length(fred);N temp<-!is.na(match(1:N,seq(1,N-216,217))) fred<-fred[!temp]# 329616 noaa<-cbind(as.numeric(substr(fred,1,6)), as.numeric(substr(fred,7,12)),as.numeric(substr(fred,13,18)),as.numeric(substr(fred,19,24)), as.numeric(substr(fred,25,30)), as.numeric(substr(fred,31,36)),as.numeric(substr(fred,37,42)),as.numeric(substr(fred,43,48)), as.numeric(substr(fred,49,54)), as.numeric(substr(fred,55,60)),as.numeric(substr(fred,61,66)),as.numeric(substr(fred,67,72)) ) rm(fred) temp<-(noaa== -32768);noaa[temp]<-NA dim(noaa) #329616 12 # 329616*12/2592 =1526 ##old<-noaa noaa<- c(t(noaa)) noaa<-array(noaa,dim=c(2592,length(noaa)/2592 ) );dim(noaa)# 2592 1526 noaa<-t(noaa) #1526 2592 ##save(noaa,file="d:/climate/data/noaa/noaa.tab") ############################# ###MAKE STATION FILE FROM MANUALLY UNZPPED FILE (GHCN V2 ADJ) ################################################ ##download and unzip from ## http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean_adj.Z ## to ## d:/climate/data/jones/ghcn/v2.mean_adj ##unadjusted from http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean.Z ## there is also a "raw" version v2a<-readLines("d:/climate/data/jones/ghcn/v2.mean_adj");length(v2a) # 420447 versus 590152 id<-as.numeric(substr(v2a,1,12)) wmo<-as.numeric(substr(v2a,4,8)) year<-as.numeric(substr(v2a,13,16)) v2d<-cbind(as.numeric(substr(v2a,17,21)),as.numeric(substr(v2a,22,26)),as.numeric(substr(v2a,27,31)) , as.numeric(substr(v2a,32,36)),as.numeric(substr(v2a,37,41)),as.numeric(substr(v2a,42,46)) , as.numeric(substr(v2a,47,51)),as.numeric(substr(v2a,52,56)),as.numeric(substr(v2a,57,61)) , as.numeric(substr(v2a,62,66)),as.numeric(substr(v2a,67,71)),as.numeric(substr(v2a,72,76)) ) v2d<-cbind(year,v2d,id,wmo) rm(id);rm(wmo);rm(year);rm(v2a) temp<-(v2d== -9999);sum(temp) v2d[temp]<-NA dim(v2d)# 420447 15 month0<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec") dimnames(v2d)[[2]][2:13]<-month0 #save(v2d,file="d:/climate/data/jones/ghcn/v2adj.tab") ################## ##TEST WITH STATION 29612 ########################### #load("d:/climate/data/noaa/noaa.tab") source("http://www.climateaudit.org/scripts/gridcell/load.stations.txt") id0<-29612 #Barabinsk, Russia info<-stations[stations$wmo==id0,];info # id site wmo version long pop lat #1771 222296120000 BARABINSK 29612 0 78.37 37 55.33 jones(info$lat,info$long) #484 stations$jones<-jones(stations$lat,stations$long) test<-tapply(!is.na(stations$id),factor(stations$jones,levels=1:2592),sum) test[is.na(test)]=0 hist(test,breaks=c(0:10,seq(15,105,5))) #collect stations in gridbox #this averages all measurements and doesn't (at present) try to first average a wmo station and then a second station #rm(v2d);method<-"raw";load("d:/climate/data/jones/ghcn/v2d.tab") #this does the 2nd version method<-"adj"; load("d:/climate/data/jones/ghcn/v2adj.tab") #table v2d id1<-stations$wmo[stations$jones== jones(info$lat,info$long)] temp<-!is.na(match(v2d[,"wmo"],id1));sum(temp) K<-length(unique(v2d[temp,"id"])); if (K>0) { tsite<-NULL for (j in 1:K) { temp1<- (v2d[temp,"id"] == unique(v2d[temp,"id"])[j]) tt<-fill.array(v2d[temp,1:13][temp1,]) tsite<-ts.union(tsite,ts(as.matrix(c(t(tt))),start=c(min(as.numeric(row.names(tt))),1),freq=12) ) } tsite<-round(cbind(tsite,apply(tsite,1,mean,na.rm=TRUE))/10,2) dimnames(tsite)[[2]]<-c(unique(v2d[temp,"id"]),"avg") time0<-unlist(time(tsite)) temp<- (time0>=1961)&(time0<1991);sum(temp) N<-nrow(tsite)/12 m0<- apply( t( array(tsite[temp,"avg"],dim=c(12,N))) ,2,mean,na.rm=TRUE) months<-rep(1:12,N) tsite.anom<-round( (tsite-m0[months]),2) } gridcell3<-ts(noaa[,jones(info$lat,info$long)],start=c(1880,1),freq=12)/100 if(method=="raw"){tsite.anom.raw<-tsite.anom;tsite.raw<-tsite} if(method=="adj"){tsite.anom.adj<-tsite.anom;tsite.adj<-tsite} combine<-ts.union(gridcell3,tsite.anom.raw[,"avg"],tsite.anom.adj[,"avg"]); a<-tsp(combine) dimnames(combine)[[2]]<-c("grid","ghcn.raw","ghcn.adj") ### ## grid ghcn.raw ghcn.adj ## grid 1.0000000 0.9992223 0.9922494 ## ghcn.raw 0.9992223 1.0000000 0.9944624 ## ghcn.adj 0.9922494 0.9944624 1.0000000 #FIGURE TO COMPARE GRID AND RAW GHCN png("d:/climate/images/gridcell/barabinsk1.png",width=640,height=480) nf<-layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) par(mar=c(0,4,1,1)) plot(seq(a[1],a[2],1/a[3]),combine[,1],type="l",xlab="",ylab="Grid",axes=FALSE) axis(side=1,labels=FALSE);axis(side=2,las=1);box() par(mar=c(3,4,0,1)) plot(seq(a[1],a[2],1/a[3]),combine[,2],type="l",xlab="",ylab="GHCN",axes=FALSE) axis(side=1);axis(side=2,las=1);box() dev.off() #GHCN minus ADJUSTED #png("d:/climate/images/gridcell/barabinsk.noaa.1.png",width=600,height=480) plot(seq(a[1],a[2],1/a[3]),combine[,3],type="l",col="grey70") lines(seq(a[1],a[2],1/a[3]),filter(combine[,3],truncated.gauss.weights(25)),col="red") text(1940,.65,"NOAA Gridcell minus GHCN v2", font=2,pos=4) #dev.off() #GHCN minus RAW #png("d:/climate/images/gridcell/barabinsk.noaa.2.png",width=600,height=480) #plot(seq(a[1],a[2],1/a[3]),combine[,3],type="l",col="grey70") #lines(seq(a[1],a[2],1/a[3]),filter(combine[,3],truncated.gauss.weights(25)),col="red") #text(1940,.65,"NOAA Gridcell minus GHCN v2 ('Raw')", font=2,pos=4) #dev.off() #FIGURE TO COMPARE RAW AND ADJ GHCNGHCN png("d:/climate/images/gridcell/barabinsk2.png",width=640,height=300) nf<-layout(1);# par(mar=c(3,4,1,1)) plot(seq(a[1],a[2],1/a[3]),combine[,2]-combine[,3],type="l",xlab="",ylab="",axes=FALSE,col="grey70") lines(seq(a[1],a[2],1/a[3]),filter(combine[,2]-combine[,3],truncated.gauss.weights(25)),col="red") axis(side=1);axis(side=2,las=1);box() text(1940,.5,"GHCN Raw minus Adjusted", font=2,pos=4) dev.off()