######################### #UTILITY FUNCTIONS TO DEAL WITH TIME SERIES ISSUES CHARACTERISTIC OF TEMPERATURE ########################## #################### ###REQUIRES #################### #package ncdf #"d:/climate/data/gridcell/hadcrut3/HadCRUT3.nc" #5x5 monthly gridded from CRU #downloaded and unzipped from http://hadobs.metoffice.com/hadcrut3/data/HadCRUT3.nc #"d:/climate/data/gridcell/giss/calannual-1880-2005.nc" #2x2 degree annual gridded from GISS # ftp://data.giss.nasa.gov/pub/gistemp/netcdf/calannual-1880-2005.nc #"d:/climate/data/gridcell/ghcn/v2d.tab" #v2.mean unsipped from GHCN # http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean.Z ## array arranged by year (row) and month (cols) with possibly missing years ##in put array is 13 cols with first column year month0<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec") #JONES LOOKUP jones<-function(lat,long) { i<-18-floor(lat/5); j<- 36+ceiling(long/5); jones<-72*(i-1) + j jones } fill.array<-function(X,method="column") { if (method=="row.name") {xyear<-as.numeric(row.names(X));X<-X[,c(13,1:12)] } if(class(X)=="matrix") {xyear<-X[,1];X<-X[,2:13]} if (class(X)=="data.frame") {xyear<-X$year;X<-as.matrix(X[,2:13])} years<-min(xyear,na.rm=TRUE):max(xyear,na.rm=TRUE) fill.array<-array(NA,dim=c(length(years),12)) temp<-!is.na( match(years,xyear) ) fill.array[temp,]<-X row.names(fill.array)<-years dimnames(fill.array)[[2]]<-month0 fill.array } ts.array<-function(x) { year1<- floor(tsp(x)[1]):floor(tsp(x)[2]) month1<- round(12*(tsp(x)[1:2]%%1),0) x<-c(rep(NA,month1[1]),x ,rep(NA, 11-month1[2])) ts.array<-t(array(x,dim=c(12,length(x)/12)) ) row.names(ts.array)<-year1 ts.array } #X is an array or time series; M is minimum to return an average #allows for time series; array with 13 cols, year in col 1; array with 12 cols ts.annavg<-function(X,M=6) { if(class(X)=="ts") {start0<-floor(tsp(X)[1]);X<-ts.array(X)} else if( ((class(X)=="matrix")|(class(X)=="data.frame") ) &(ncol(X)==12)) start0<-as.numeric(row.names(X)[1]) else {X start0<-X[1,1]; X<-X[,2:13]} count<-apply(!is.na(X),1,sum);temp<-(count1) { tsite<-NULL for (j in 1:M) { temp1<-(id[temp]==levels(factor(id[temp]))[j]) year1<- c(min(v2d[temp,"year"][temp1]),max(v2d[temp,"year"][temp1])) tt<-array(NA,dim=c(year1[2]-year1[1]+1,12)) temp2<-!is.na(match (year1[1]:year1[2],v2d[temp,"year"][temp1])) tt[temp2,]<-v2d[temp,2:13][temp1,] tt<-ts( c(t(tt)),start=c(year1[1],1),freq=12) tsite<-ts.union(tsite,tt)} tt<-round(ts( apply(tsite,1,mean,na.rm=TRUE), start=c(tsp(tsite)[1],1),freq=12),0) } if(M==1) { year1<- c(min(v2d[temp,"year"]),max(v2d[temp,"year"])) tt<-array(NA,dim=c(year1[2]-year1[1]+1,12)) temp2<-!is.na(match (year1[1]:year1[2],v2d[temp,"year"])) tt[temp2,]<-v2d[temp,2:13] tt<-ts( tt,start=c(year1[1],1),freq=12) } if(M==0) tt<-NA tt } ###################### ####GRIDDED VERSIONS #1. HadCRUT3 #2. GISS 2007 #3. HadCRU2 #4. NOAA. ##### #TO DO #4. CRUTEM3 #5. Crutem2 #6. HadCRu1 ####################### ################## ##LOAD HADCRUT3 ###### #function to return annualized time series from HadCRU3 for lat long #requires HAdCru3 be in place - ncdf object instr #returns monthly series squared off to year with NAs ##library(ncdf) ##url<-"d:/climate/data/gridcell/hadcrut3/HadCRUT3.nc" #previously downloaded and unzipped from http://hadobs.metoffice.com/hadcrut3/data/HadCRUT3.nc #test <- v$var[[1]] #test$size # 72 36 1 1883 #test$dim[[1]]$vals #[1] -177.5 -172.5 -167.5 -162.5 -157.5 -152.5 -147 #test$dim[[2]]$vals # [1] -87.5 -82.5 -77.5 -72.5 -67.5 -62.5 -57.5 read.hadcru3<-function (lat,long,url="d:/climate/data/gridcell/hadcrut3/HadCRUT3.nc"){ v<-open.ncdf(url) instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #dim(instr)# [1] 72 36 1883 #this is organized in 72 longitudes from -177.5 to 177.5 and 36 latitudes from -87.5 to 87.5 j<- 19+ floor(lat/5 +.01); i<- 37+ floor(long/5 +.01); gridcell3<-instr[i,j,] index<-length(gridcell3)%%12 gridcell3<-c(gridcell3,rep(NA,12-index) ) close.ncdf(v) read.hadcru3<-ts(gridcell3,start=c(1850,1),freq=12) read.hadcru3 } ##station.hadcru3<-read.hadcru3(lat=52.5,long=2.5) #returns monthly data ###################### # LOAD GISS GRIDDED ################# ##this is ANNUAL ONLY #this already downloaded manually #requires library(ncdf) #already loaded #loc<-file.path("d:/climate/data/gridcell/giss","calannual-1880-2005.nc") #v<-open.ncdf(loc) #length 10 #[1] "id" "ndims" "natts" "unlimdimid" # [5] "filename" "varid2Rindex" "writable" "dim" # [9] "nvars" "var" #test <- v$var[[1]] #[1] "id" "name" "ndims" "natts" "size" "prec" "dimids" #[8] "units" "longname" "dims" "dim" "varsize" "unlim" "missval" #[15] "hasAddOffset" "hasScaleFact" #test$size # 180 90 126 #test$dim[[1]]$vals #[1] -179 -177 -175 .... #test$dim[[2]]$vals # [1] -89 -87 -85 ... #test$dim[[3]]$vals #1880 2005 read.giss.grid<-function(lat,long,url="d:/climate/data/gridcell/giss/calannual-1880-2005.nc" ) { i<- floor(long/2 +.0001) +91; j<- floor(lat/2 +0.001) +46 v<-open.ncdf(url) #length 10 giss <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #dim(giss)# [1] 180 90 126 #this is organized in 180 longitudes from -179 to 179 and 90 latitudes from -89 to 89 read.giss.grid<-ts( giss[i,j,],start=1880) close.ncdf(v) read.giss.grid } ##station.gissgrid<-read.giss.grid(lat=57.5,long=77.5) ###HadCRU2 read.hadcru2<-function (lat,long,url="d:climate/data/gridcell/hadcrut2/hadcruv.tab"){ load(url) k<-jones(lat,long) gridcell3<-v[,k] index<-length(gridcell3)%%12 gridcell3<-ts (c(gridcell3,rep(NA,12-index) ),start=c(1856,1),freq=12) rm(v) read.hadcru3<-gridcell3/100 read.hadcru3 } ##station.hadcru2<-read.hadcru2(57.5,77.5) ########################## ####GLOBAL STATIONS #1. GHCN v2 2007 #monthly #2. GISS #monthly #3. GSN #daily ##TO DO #4. Jones 1991 #5. Jones 1994 #6. Jones 1996 #7. GHCN v1 ######################## #in each case, a list is returned #1 - monthly "raw" time series; #2 - monthly anomaly time series using calculated 1961-1990 normal if available, otherwise GHCN normal #3 - normal or "GHCN" if this is used #for daily data, M1= 14 readings are required to take monthly average (adjustable parameter) ########################## ########### #1. GHCN V2 VERSION AS ANOMALY #download zipped file from # http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean.Z #script to collate this data into R-table to save time since used often #v2<-readLines("d:/climate/data/gridcell/ghcn/v2.mean");length(v2) # 590152 #id<-as.numeric(substr(v2,1,12)) #wmo<-as.numeric(substr(v2,4,8)) #year<-as.numeric(substr(v2,13,16)) #JFM<-cbind(as.numeric(substr(v2,17,21)),as.numeric(substr(v2,22,26)),as.numeric(substr(v2,27,31)) ) #AMJ<-cbind(as.numeric(substr(v2,32,36)),as.numeric(substr(v2,37,41)),as.numeric(substr(v2,42,46)) ) #JAS<-cbind(as.numeric(substr(v2,47,51)),as.numeric(substr(v2,52,56)),as.numeric(substr(v2,57,61)) ) #OND<- cbind(as.numeric(substr(v2,62,66)),as.numeric(substr(v2,67,71)),as.numeric(substr(v2,72,76)) ) #v2d<-cbind(year,JFM,AMJ,JAS,OND,id,wmo) #rm(JFM);rm(AMJ);rm(JAS);rm(OND);rm(id);rm(wmo);rm(year) #temp<-(v2d== -9999);sum(temp) #v2d[temp]<-NA #dim(v2d)# 590152 15 #save(v2d,file="d:/climate/data/gridcell/ghcn/v2d.tab") ###adjusted version http://www1.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean_adj.Z # v2a<-readLines("d:/climate/data/gridcell/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)) #v2adj<-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)) ) #v2adj<-cbind(year,v2adj,id,wmo) #rm(id);rm(wmo);rm(year);rm(v2a) #temp<-(v2adj== -9999);sum(temp) #v2adj[temp]<-NA #dim(v2adj)# 420447 15 #month0<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec") #dimnames(v2adj)[[2]][2:13]<-month0 #save(v2adj,file="d:/climate/data/jones/ghcn/v2adj.tab") #from calculated anomaly object #load("d:/climate/data/jones/ghcn/v2d.tab") # #590015 rows - large data set with calculated anomalies #function to return time series of wm0 versions; requires v3 (anomaly version)( #this could be back-integrated to work with non-normalized data #data is monthly #this returns M time series (monthly or annual) #requires v3 to be present #allows for ragged representation #starting point is matrix row- years; cols -month; possibly missing years #the data comes as year-month arrays with possibly missing years #monthly or annual #returns combined or separate #required read.ghcn<-function(id0,v2=v2d) { #test<-try(dim(v2d));if(class(test)=="try-error") load (file.path( url,"v2d.tab") ) #checks to see if data set is already loaded temp<-(v2[,"wmo"]==id0); read.ghcnv2<-NA K<-length(unique(v2[temp,"id"])); if (K>0) { tsite<-NULL for (j in 1:K) { temp1<- (v2[temp,"id"] == unique(v2[temp,"id"])[j]) tt<-fill.array(v2[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(v2[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) read.ghcn<-list(tsite,tsite.anom,m0 )# raw monthly collation and total; names(read.ghcn)<-c("raw","anom","normal") } read.ghcn } # station.ghcn<-read.ghcnv2(id0) #station.normal<-station.ghcn[[3]] ############################ ############################## ###GISS VERSION - CURRENT # also see scripts at http://x256.org/~hb/giss.r2 # data sets 0 ;0- raw uncombined; 1- combined; 2 -adjusteed download_html <- function(url) { download.file(url, "temp.html"); html_handle <- file("temp.html", "rt"); html_data <- readLines(html_handle); close(html_handle); unlink("temp.html"); return(html_data); } url_escape <- function(string) { string <- gsub(" ", "%20", string) string <- gsub("#", "%23", string) string <- gsub("&", "%26", string) string <- gsub("=", "%3D", string) string <- gsub("\\?", "%3F", string) return(string); } # station_name can either be a name, or a 5-digit WMO number. # If there is no match, NA is returned. # This should work as long as NASA don't change the URLs they use for providing data. download_giss_data <- function(station_name, match_num=1,dset=0) { giss_url <- paste("http://data.giss.nasa.gov/cgi-bin/gistemp/findstation.py?datatype=gistemp&data_set=",dset,"&name=", url_escape(station_name), "&submit=1", sep = ""); my_data <- download_html(giss_url); urls <- grep("gistemp_station.py", my_data, value=1); # if the station_name is actually a WMO number, throw away any URLs that match the number in the wrong place of the extended number scheme. if( length(grep("^[0-9]{5}$", station_name)) != 0 ) { urls <- grep(paste("\\?id=[0-9]{3}", station_name, sep=""), urls, value=1); } if( length(urls) == 0 || match_num > length(urls) ) return(NA); K<-length(urls) url <- paste("http://data.giss.nasa.gov", gsub("^.*\"(/cgi-bin/gistemp/gistemp_station.py[^\"]*)\".*$", "\\1", urls), sep=""); #varied from Nicholas to leave as vector id<-substr(url,65,76) urls.new<-rep(NA,K);url.new<-urls.new; download_data<-rep(list(NA),K) for (j in 1:K){ my_data <- download_html(url[j]); urls.new[j] <- grep("monthly data as text", my_data, value=1); url.new[j] <- gsub("^.*\"(http://data.giss.nasa.gov/work/gistemp/STATIONS/[^\"]*)\".*$", "\\1", urls.new[j]); my_handle <- url(url.new[j], "rt"); download_data[[j]] <- read.table(my_handle, skip=1);#instead of NIcholas method names(download_data[[j]])[1]<-"year" close(my_handle); } names(download_data)<-id download_giss_data<-rep(list(NA),3);names(download_giss_data)<-c("raw","anom","normal") tsite<-NULL for (j in 1:K) { temp<-(download_data[[j]] == 999.9); download_data[[j]][temp]<-NA tt<-fill.array(download_data[[j]][,1:13]) #if(method=="annual") tsite<-ts.union(tsite,ts.annavg(tt,M=6)) tsite<-ts.union(tsite,ts(as.matrix(c(t(tt))),start=c(min(as.numeric(row.names(tt))),1),freq=12) )#if(method=="monthly") } tsite<-cbind(tsite,apply(tsite,1,mean,na.rm=TRUE)) dimnames(tsite)[[2]]<-c(id,"avg") download_giss_data[[1]]<-tsite years<-floor(tsp(tsite)[1]):floor(tsp(tsite)[2]) N<-nrow(tsite)/12 chron<- t( array(tsite[,"avg"],dim=c(12,N))) temp<-!is.na(match(years,1961:1990)) m0<- apply(chron[temp,],2,mean,na.rm=TRUE) months<-rep(1:12,N) names(m0)<-month0 download_giss_data[[3]]<-m0 tsite.anom<-round(tsite-m0[months],1) download_giss_data[[2]]<-tsite.anom; download_giss_data } # station.giss<-download_giss_data(id0) ########################## ################################# #### GSN ##function to return annualized anomaly series from GSN ##daily historical information ##individualized files which may be recent subst only, in which case GHCN normal used #e.g Barabinsk 3332 67 ##http://www1.ncdc.noaa.gov/pub/data/ghcn/daily/gsn/22200029612.dly read.gsn<-function(id0,M1=14) { read.gsn<-rep(list(NA),3) names(read.gsn)<-c("raw","anom","m0") url<-"http://www1.ncdc.noaa.gov/pub/data/ghcn/daily/gsn" loc<-file.path(url,paste("222000",id0,".dly",sep="")) width0<-c(6,5,4,2,4,rep(c(5,3),31)) test<-try( read.fwf(loc,widths=width0,nrow=-1,fill=TRUE)) if (class(test)=="try-error") read.gsn<-NA else { test<-test[,c(1:5,seq(6,66,2))] temp3<- (test== -9999);test[temp3]<-NA levels(test[,5]) # [1] "PRCP" "TMAX" "TMIN" #dim(test)# 1052 36 names(test)[1:5]<-c("code","wmo","year","month","type") test<-test[!(test$type=="PRCP"),] # 1995 2007 test<-test[test$year<2007,] #not enough 2007 data yet test$count<-apply(!is.na(test[,6:36]),1,sum) test$avg<-NA; test$avg[test$count>M1]<-apply(test[test$count>M1,6:36],1,mean,na.rm=TRUE) test$time<-100*test$year+test$month g<-function(x,y){g<-100*x+y;g} years<- min(test$year):max(test$year) gsn<-sort(outer(years,1:12,g) ); gsn<-data.frame(gsn);names(gsn)<-"time" temp1<-(test$type=="TMIN") temp<-!is.na(match(gsn$time,test$time[temp1])) gsn$tmin<-NA; gsn$tmin[temp]<-test$avg[temp1]/10 temp<-!is.na(match(gsn$time,test$time[!temp1])) gsn$tmax<-NA; gsn$tmax[temp]<-test$avg[!temp1]/10 gsn$avg<-(gsn$tmax+gsn$tmin)/2 if (sum(!is.na(match(1961:1990,years))) >20 ) { temp<- !is.na(match(floor(gsn$time/100), 1961:1990) ) m0<- apply( t( array(gsn$avg[temp],dim=c(12,30))) ,2,mean,na.rm=TRUE) #times 100 } else {m0<-station.normal} gsn$anom<- NA; if (sum (!is.na(m0))==12) gsn$anom<- round( (gsn$avg -m0[gsn$time%%100]),2) read.gsn[[1]]<-ts(round(gsn$avg,1),start=c(years[1],1),freq=12) read.gsn[[2]]<- ts(gsn$anom,start=c(years[1],1),freq=12) if (sum(!is.na(match(1961:1990,years))) >20 ) read.gsn[[3]]<-m0 else read.gsn[[3]]<-"GHCN" } read.gsn } # station.gsn<-read.gsn(id0) ################################# ############ ##### ##SPECIALIZED ARCHIVES ##### ##RUSSIA #1. meteo.ru #2. unaami #3. G02141 - Radionov #TO DO #4. NDP048 #5. NDP040 #6. CHINA NDP039 ############# ###1. METEO.RU #35663 1945 9 17 0 50 0 115 0 184 0 0 2 0 #35663 1945 9 18 0 70 0 135 0 213 0 0 2 0 #35663 1945 9 19 0 87 0 148 0 216 0 0 2 0 #35663 1945 9 20 0 94 0 137 0 184 0 0 2 0 #NOTE: possible inconsistency in TMIN, TMAX as many values have TMIN and TMID, but no TMAX #avg calculated as average of TMAX and TMIN where available #reference: http://meteo.ru/english/data_temperat_precipitation/structure.php #daily data in long table #the years are sometimes missing and this needs to be allowed for #on MArch 14, 2007, changed from V4 to V6; also a column added #V6 version - last column alphanumeruc #25563 1898 1 1 9 999.9 9 999.9 9 999.9 9 999.9 9 9 9999 #25563 1898 1 2 9 999.9 9 999.9 9 999.9 9 999.9 9 9 9999 #25563 1898 1 3 9 999.9 9 999.9 9 999.9 9 999.9 9 9 9999 read.meteo<-function(id0) { ##id0<-35663 name0<-paste("F",id0,"v6.zip",sep="") #formerly V4.ZIP url<-"ftp://ftp.meteo.ru/okldata" test<-try(download.file( file.path(url,name0),name0, "auto", quiet = FALSE, cacheOK = TRUE, mode = "wb")); if(class(test)=="try-error") read.meteo<-test else { read.meteo<-rep(list(NA),3) names(read.meteo)<-c("raw","anom","normal") data_handle <- unz(name0,paste("F",id0,"v6",sep="") ,"r") #formerly V4.M02 Data <- scan(data_handle, what = list(c(rep(0,14),"")), comment.char = "\x1a"); #Data <- read.table(data_handle, colClasses=c(rep("numeric",14),"character")) #Data <- read.fwf(data_handle, widths=c(5,5,3,3,2,6,2,6,2,6,2,6,2,2,5)) close(data_handle) Data<-unlist(Data); N<-length(Data);temp<-!is.na(match(1:N, seq(15,N,15) )) Data<-Data[!temp] Data<-as.numeric(Data) Data<-t(array(Data,dim=c(14,length(Data)/14))) Data<-Data[,c(1:4,6,8,10)] Data<-data.frame(Data) names(Data)<-c("id","year","month","day","tmin","tmid","tmax") temp<-(Data== 999.9);Data[temp]<-NA Data$time<-100*Data$year+Data$month Data$tmean<-apply(Data[,c("tmin","tmax")],1,mean) #requires both to be present year1<-min(Data$year):max(Data$year) g<-function(x,y){g<-100*x+y;g} time0<-sort(c(outer(year1,1:12,g))) chron<-tapply(Data$tmean,factor(Data$time,levels=time0),mean,na.rm=TRUE) chron<-t( array(chron,dim=c(12,length(chron)/12) ) ) row.names(chron)<-year1;dimnames(chron)[[2]]<-month0 read.meteo[[1]]<- ts(c(t(chron)),start=c(year1[1],1),freq=12) temp<- !is.na(match(as.numeric(row.names(chron)),1961:1990) ) if(sum(temp)>20) {m0<-apply(chron[temp,],2,mean,na.rm=TRUE);read.meteo[[3]]<-m0} else m0<-station.normal test<-scale(chron,center=m0,scale=FALSE) read.meteo[[2]]<-ts( c(t(test)),start=c(min(year1),1),freq=12 ) } read.meteo } ##station.meteo<-read.meteo(id0) #ftp://ftp.meteo.ru/okldata/F25563v6.ZIP ##################### #2. UNAAMI ### #selection of data described at url # [1] "//20744,52.70,72.37,15,222,MALYE KARMAKULY,TUNDRA" # [2] "-10.3,-21.3,-11.6, -8.8, 1.7, 2.6, 6.3, 5.5, 4.0, -3.2,-10.4, -9.4,//1897" # [3] "-18.1,-17.2,-15.5,-10.3, -5.3, 1.0, 8.7, 99.9, 99.9, -4.1, -9.3,-18.9,//1898" #monthly raw read.unaami<-function(id0){ url<-"http://www.unaami.noaa.gov/analyses/sat" loc<-file.path(url, paste("sat",id0,".d",sep="")) test<-try(read.table(loc,skip=1,sep=",")) if(class(test)=="try-error") read.unaami<-NA else { read.unaami<-rep(list(NA),3) temp3<-(test== 99.9) test[temp3]<-NA x<-as.character(test[,13]) test[,13]<-as.numeric(substr(x,3,nchar(x)) ) read.unaami[[1]]<-ts(c(t(test[,1:12])),start=c(test[1,13],1),freq=12) temp<-!is.na(match(test[,13],1961:1990)) m0<-apply(test[temp,1:12],2,mean,na.rm=TRUE);names(m0)<-month0 read.unaami[[3]]<-m0 test[,1:12]<-scale(test[,1:12],center=m0,scale=FALSE) read.unaami[[2]]<-ts(c(t(test[,1:12])),start=c(test[1,13],1),freq=12) } read.unaami } #station.unaami<-read.unaami(20744) #none for id0 ########################### ##3. G02141 by RADIUNOV ## #this has a missing call structure with mix of name short forms and wmo numbers; requires a manually collated list that I've archived at CA #monthly data in long form #temperature in 9th column # [1] "25744 1961 1 -1 -1 1 62.48 166.22 -18.8 1011.7 999.9 999.9 7.2 4.1 81 999.99 999.99 9999.9 36.0 999.99 999.99 kamenskoe" # [2] "25744 1961 2 -1 -1 1 62.48 166.22 -20.0 1016.9 999.9 999.9 5.7 4.3 77 999.99 999.99 9999.9 25.0 999.99 999.99 kamenskoe" # [3] "25744 1961 3 -1 -1 1 62.48 166.22 -27.8 1022.8 999.9 999.9 5.2 2.3 76 999.99 999.99 9999.9 17.0 999.99 999.99 kamenskoe" read.radionov<-function(id0) { url<-"http://www.climateaudit.org/data/station" url<-"d:/climate/data/station/russia" radionov<-scan(file.path(url,"radionov.list.txt"),skip=1,what="") h<-unlist(strsplit(radionov, "\\.")) h<-array(h,dim=c(2,length(h)/2)); h<-t(h); h<-h[,2] temp1<- (as.numeric(h)==id0) if(sum(temp1)==0) read.radionov<-NA else { read.radionov<-rep(list(NA),3) names(read.radionov)<-c("raw","anom","normal") loc<-file.path("ftp://sidads.colorado.edu/pub/DATASETS/NOAA/G02141",paste("uni",radionov[temp1],"dat",sep=".")) test<-read.fwf(loc,widths=c(5,5,3,26,6)) test<-test[,c(1:3,5)] names(test)<-c("wmo","year","month","mean") tt<-t(array(test$mean,dim=c(12,nrow(test)/12) )) years<-min(test$year):max(test$year);row.names(tt)<-years read.radionov[[1]]<-ts(c(t(tt)),start=c(years[1],1),freq=12) temp<-!is.na(match(years,1961:1990)) m0<-apply(tt[temp,],2,mean,na.rm=TRUE) read.radionov$normal<-m0 tt<-scale(tt,center=m0,scale=FALSE) read.radionov[[2]]<-ts(c(t(tt)),start=c(years[1],1),freq=12) } read.radionov } #station.radionov<-read.radionov(25744)