#http://cdiac.ornl.gov/ftp/ndp048/f38974.dat #http://cdiac.ornl.gov/ftp/ndp048/f38987.dat ndp40_colnames=c("WMO","Year","Month","Day","TFLAG","TMIN","QTMIN","TMID","QTMID","TMAX","QTMAX","R","CR","QR","DATA_FLAGS") ndp40_colwids=c(6,6,4,4,2,6,2,6,2,6,2,6,2,2,5) ndp40_38974 <- read.fwf("f38974.dat", header=F, na.strings=" 999.9", widths=ndp40_colwids, col.names=ndp40_colnames) ndp40_38987 <- read.fwf("f38987.dat", header=F, na.strings=" 999.9", widths=ndp40_colwids, col.names=ndp40_colnames) #http://cdiac.ornl.gov/ftp/ndp048/38974.dat.Z #http://cdiac.ornl.gov/ftp/ndp048/38987.dat.Z ndp48_colnames=c("WMO", "Year", "Month", "Day", "Hour", "Stuff1", "Temp", "TQC", "Stuff2") ndp48_colwids=c(5,4,2,2,2,74,4,1,34) ndp48_38974 <- read.fwf("38974.dat", header=F, na.strings="9999", widths=ndp48_colwids, col.names=ndp48_colnames) ndp48_38987 <- read.fwf("38987.dat", header=F, na.strings="9999", widths=ndp48_colwids, col.names=ndp48_colnames) for (i in 1:4) { if (i == 1) { # GHCN 229389740011 synop <- ndp40_38974 dupe <- 1 file <- "ndp40_38974.txt" } else if (i == 2) { # GHCN 229389740021 synop <- ndp40_38987 dupe <- 1 file <- "ndp40_38987.txt" } else if (i == 3) { # GHCN 229389740011 synop <- ndp48_38974[,c(1,2,3,4,5,7)] dupe <- 1 file <- "ndp48_38974.txt" } else if ( i == 4 ) { # GHCN 229389740021 synop <- ndp48_38987[,c(1,2,3,4,5,7)] dupe <- 2 file <- "gsod_38987.txt" } if (i == 3 || i == 4 ) { # ndp048 id <- synop[1,1] t1 <- synop # TOB - shift the temps down one 'slot' t2 <- t1[,6] t1[2:nrow(t1),6] <- t2[1:(length(t2)-1)] # trim data to 1936-1989 t1 <- t1[2:nrow(t1),] mask <- t1[,2] >= 1936 & t1[,2] <= 1989 t1 <- t1[mask,] range <- max(t1$Year) - min(t1$Year) + 1 d_mean <- array(NA,dim=c(range,12,31)) m_mean <- matrix(NA,range,12) d_mean <- tapply(t1[,6], list(t1[,2],t1[,3],t1[,4]), mean, na.rm=T) } else if (i == 1 || i == 2 ) { #ndp040 id <- synop[1,1] t1 <- synop # trim data to 1936-1989 t1 <- t1[2:nrow(t1),] mask <- t1[,2] >= 1936 & t1[,2] <= 1989 t1 <- t1[mask,] range <- max(t1$Year) - min(t1$Year) + 1 d_mean <- array(NA,dim=c(range,12,31)) m_mean <- matrix(NA,range,12) d_mean <- tapply(t1[,8], list(t1[,2],t1[,3],t1[,4]), mean, na.rm=T) d_mean <- d_mean * 10 } for (y in 1:range) { for (m in 1:12) { t = mean(d_mean[y,m,], na.rm=T) m_mean[y,m] <- as.integer( t ) if ( (! is.na(t)) ) { if (t < 0) { m_mean[y,m] <- as.integer( t - 0.5) } else { m_mean[y,m] <- as.integer( t + 0.5) } } } } cc <- 229 stid <- paste(cc,id,dupe,sep="") stnid <- as.numeric(rep(stid,range)) yrs <- seq(1936:1989) table <- cbind(stnid,yrs+1935,m_mean) write.table(table, file, row.names=F, col.names=F) } ndp40_38974a <- read.table("ndp40_38974.txt",header=F) # 229389740011 ndp40_38987a <- read.table("ndp40_38987.txt",header=F) # 229389740021 ndp48_38974a <- read.table("ndp48_38974.txt",header=F) # 229389740011 ndp48_38987a <- read.table("ndp48_38987.txt",header=F) # 229389740021 ghcn_389740011 <- read.table("389740011-ghcn.txt",header=F, na.strings="-9999") ghcn_389740020 <- read.table("389740020-ghcn.txt",header=F, na.strings="-9999") ghcn_389740021 <- read.table("389740021-ghcn.txt",header=F, na.strings="-9999") sum(abs(ghcn_389740021[,3:14] - ndp40_38987a[,3:14]),na.rm=T) # 0 sum(abs(ghcn_389740021[,3:14] - ndp48_38987a[,3:14]),na.rm=T) # 92 sum(abs(ghcn_389740011[,3:14] - ndp40_38974a[,3:14]),na.rm=T) # 3 sum(abs(ghcn_389740011[,3:14] - ndp48_38974a[,3:14]),na.rm=T) # 113 #---------- library(reshape) # reshape ghcn colnames(ghcn_389740011) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") colnames(ghcn_389740021) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") # #g <- ghcn_389740011[,2:14] g <- ghcn_389740021[,2:14] # gm <- melt(g, id.var="Year", variable_name="Month") frac <- (unclass(gm$Month)-0.5)/12 g2 <- cbind(frac+gm[,1],gm[,3]) g3 <- g2[order(g2[,1]), ] # reshape ndp colnames(ndp40_38974a) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") colnames(ndp48_38974a) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") colnames(ndp40_38987a) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") colnames(ndp48_38987a) <- c("ID","Year", "Jan","Feb","Mar", "Apr", "May", "Jun", "Jul","Aug","Sep","Oct","Nov","Dec") # n <- ndp40_38974a[,2:14] #n <- ndp48_38974a[,2:14] #n <- ndp40_38987a[,2:14] #n <- ndp48_38987a[,2:14] # nm <- melt(n, id.var="Year", variable_name="Month") frac <- (unclass(nm$Month)-0.5)/12 n2 <- cbind(frac+nm[,1],nm[,3]) n3 <- n2[order(n2[,1]), ] #plot(g3[,2]/10 ~ g3[,1],type="l", main="GHCN 38970021", xlab="Years", ylab="Deg C") plot((g3[,2]-n3[,2]) ~ n2[,1],type="l", main="GHCN 38970021 - NDP040 38974", xlab="Years", ylab="1/10th Deg C")