############################################################################### # trb 0.07 : 2010 07 18 : Performance Refactor (22 seconds with data read) # trb 0.06 : 2010 07 17 : Performance Refactor (3x faster, still 8hrs for all) # trb 0.05 : 2010 07 14 : Gridding by Latitude and Longitude # trb 0.04 : 2010 07 13 : Proportional Latitudinal Bands # trb 0.03 : 2010 07 12 : Weighted Latitudinal Bands # trb 0.02 : 2010 07 11 : Simple Latitudinal Bands # trb 0.01 : 2010 07 10 : GHCN Simple Mean Average # http://rhinohide.wordpress.com ############################################################################### ver <- "0.07" read <- FALSE # set to FALSE if this is your first time through read <- TRUE # set to TRUE if you have previous data to read ############################################################################### # Read Data Files ############################################################################### # ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean # ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/v2.temperature.inv ############################################################################### # this flag is to indicate if we must read for the first time # or if we can read the data from previous runs if (!read) { #------------------------------------------------------------------------------ # Temperature File (GHCN formatted monthly means) #------------------------------------------------------------------------------ # mm for 'monthly means' mmwids=c(11,1,4,5,5,5,5,5,5,5,5,5,5,5,5) mmcols=c("id","dupe","year","jan","feb","mar","apr","may","jun", "jul","aug","sep","oct","nov","dec") # read the file mm <- read.fwf("../data/v2.mean", widths=mmwids, col.names=mmcols, header=FALSE, na.strings="-9999") #------------------------------------------------------------------------------ # Inventory File #------------------------------------------------------------------------------ # there are five occurrances of the pound hash in v2.temperature.inv # R scripts chokes on reading those, remove them first invcols <- c("id","name","lat","lon","alt","elev","ru","pop","terr","arid", "coast","cdist","airport","adist","vegetation","abc","bright") invwids <- c(11,32,6,8,5,5,1,5,2,2,2,2,1,2,16,1,5) inv <- read.fwf("../data/v2.temperature.inv",widths=invwids,col.names=invcols) ############################################################################### # Set up arrays for latitudinal banding ############################################################################### hindex=new.env(hash=TRUE) hlat=new.env(hash=TRUE) hlon=new.env(hash=TRUE) for (i in 1:nrow(inv)) { hindex[[as.character(inv[i,1])]] <- i hlat[[as.character(inv[i,1])]] <- inv[i,3] hlon[[as.character(inv[i,1])]] <- inv[i,4] } mmlat <- rep(NA,nrow(mm)) mmlon <- rep(NA,nrow(mm)) for (i in 1:nrow(mm)) { mmlat[i] <- hlat[[as.character(mm[i,1])]] mmlon[i] <- hlon[[as.character(mm[i,1])]] } # add inv lat lon mm1 <- cbind(mm,mmlat,mmlon) # general formula for a set number of latitude bands of equal area nbands <- 36 bands <- rep(NA,nbands) b = 90 # start with santa bands[1]=b bands[nbands+1]=-90 # rounding errors at nbands ~ 32 for (i in 1:(nbands-1)) { a = asin(sin(b*pi/180)-2/nbands) b = a * 180 / pi bands[i+1] = b } bands # 'bars' are the longitudinal edges of the grid # calculate the bars long now so we can just lookup later nbars <- nbands * 2 # this could be set to any number bars <- matrix(NA,nbands+1,nbars+1) bars[,1] = -180 # start in Howland Island and move east ;-) for (i in 1:(nbands+1)) { for (j in 1:nbars) { bars[i,j+1] = 360*(j/nbars)-180 } } bars[1,] bars[nbands+1,] #------------------------------------------------------------------------------ # assign each record to a grid point #------------------------------------------------------------------------------ mmlat <- rep(NA,nrow(mm1)) mmlon <- rep(NA,nrow(mm1)) for (i in 1:nrow(mm1)) { mmlat[i] = sum(bands>mm1[i,16]) mmlon[i] = sum(bars[mmlat[i],]= year0 & mm2[,3]<=year1 mm2 <- mm2[idx,] mm2[1,] #------------------------------------------------------------------------------ # clean up #------------------------------------------------------------------------------ mm <- NULL mm1 <- NULL mmlat <- NULL mmlon <- NULL inv <- NULL hlat <- NULL hlon <- NULL #------------------------------------------------------------------------------ # write the data #------------------------------------------------------------------------------ write.table(mm2,".Rmm2",row.names=F,col.names=T) write.table(bands,".Rbands",row.names=F,col.names=F) write.table(bars,".Rbars",row.names=F,col.names=F) write.table(c(year0,year1),".Ryears",row.names=F,col.names=F) } else { #------------------------------------------------------------------------------ # read the data #------------------------------------------------------------------------------ mm2 <- read.table(".Rmm2",header=T) bands <- read.table(".Rbands")[,1] nbands <- length(bands) -1 bars <- read.table(".Rbars") nbars <- ncol(bars) -1 yrs <- read.table(".Ryears")[,1] year0 <- yrs[1] year1 <- yrs[2] range <- year1 - year0 + 1 } # end if for data ############################################################################### # Calculate Annual Mean Average ############################################################################### #------------------------------------------------------------------------------ # initialize some vectors to hold annual results #------------------------------------------------------------------------------ lyear <- rep(NA, range) lmean <- matrix(NA,nbands,range) lcount <- matrix(NA,nbands,range) gmean <- array(NA, dim=c(nbands,nbars,range)) gcount <- array(NA, dim=c(nbands,nbars,range)) gmean <- tapply(mm2[,4], list(mm2[,1],mm2[,2],mm2[,3]), mean, na.rm=T) #------------------------------------------------------------------------------ # loop through the range of years #------------------------------------------------------------------------------ for (y in 1:range) { # loop through each year yr <- year0 + y - 1 for (i in 1:nbands) { # loop through each lat band lmean[i,y] = mean(gmean[i,,y],na.rm=T) #lcount[i,y] = 0 # TO-DO: station count per grid } lyear[y]=yr # keep track of the ranges of years } #------------------------------------------------------------------------------ # calculate the global land temperature #------------------------------------------------------------------------------ # since each band is equal area, take a simple mean ameans = colMeans(lmean,na.rm=T) mean(ameans, na.rm=T) # 15.40352 ############################################################################### # Display ############################################################################### # plot some stuff # do this nband times for (i in 1:nbands) { # had some issues when trying to plot vectors with no data if ( length(lmean[i,!is.na(lmean[i,])]) > 0 ) { img=paste("trb-",ver,"-zone-",i,".png", sep="") tit=paste("Latitudes",bands[i],"to",bands[i+1], sep=" ") png(img, width=600, height=150) par(mar=c(2, 2, 2, 1)) plot(lmean[i,] ~ lyear, ylab="Temp (C)", xlab="", type="l", main=tit, cex.main=0.8) } } # to avoid the 0 for 2010, drop last value img=paste("trb-",ver,"-zone-all.png", sep="") subt=paste("http://rhinohide.wordpress.com (ver ",ver," )",sep="") png(img, width=600, height=150) par(mar=c(2, 2, 2, 1)) plot (ameans[1:(length(ameans)-1)] ~ lyear[1:(length(ameans)-1)], type="l", xlab="", ylab="Temp (C)", main="GHCN Gridded Averaged Temps", sub=subt)