######################### #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=12) { 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<-(count0) x<-c(x, rep(NA,12-m)) years<-length(x)/12 x<-array(x,dim=c(12,years)) annavg<-apply(x,2,mean,na.rm=TRUE) return(annavg) } # Create a directory and subdirectories if they do not exist. createsubdirs<-function(dirpath1) { if( ! file.exists(dirpath1) ) { if( ! file.exists(dirname(dirpath1)) ) { # If parent directory is missing, request its creation. createsubdirs( dirname(dirpath1) ) } # the parent directory exists but not this dir dir.create( dirpath1 ) } } station.normal<-function(id0) { temp<-(wmo==id0); #v2d[temp,] M<-length(levels(factor(id[temp])));#M #4 if (M>1) { 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. CRUTEM3 #3. HadSST #4. HadCRU2 #5. GISS #6. NOAA. #7. OLD HADCRU (MBH98) #8. NDP020 #9. Kaplan SST ##### #TO DO #5. Crutem2 #6. HadCRu1 ##1. HADCRU3 ######### #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 #basis 1961-1990 ##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 #http://hadobs.metoffice.com/hadcrut3/data/HadCRUT3.nc download.hadcrut3=function() { download.file("http://hadobs.metoffice.com/hadcrut3/data/HadCRUT3.nc","temp.nc",mode="wb") handle<-open.ncdf("temp.nc") download.hadcrut3 <- get.var.ncdf( handle, handle$var[[1]]) # 1850 to Oct 2008 #1906 close(handle) download.hadcrut3 } 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 #2. CRUTEM 3 # v$var$temp$dim[[1]] longiture -177.5 to 177.5 # v$var$temp$dim[[2]] latitude -87.5 to 87.5 # v$var$temp$dim[[4]] months since 1850 in Julian days #http://hadobs.metoffice.com/crutem3/data/CRUTEM3.nc download.crutem3=function() { download.file("http://hadobs.metoffice.com/crutem3/data/CRUTEM3.nc","temp.nc",mode="wb") handle<-open.ncdf("temp.nc") download.crutem3 <- get.var.ncdf( handle, handle$var[[1]]) # 1850 to Oct 2008 #1906 close(handle) download.crutem3 } read.crutem3<-function (lat,long,url="d:/climate/data/gridcell/crutem3/crutem3.nc"){ v<-open.ncdf(url) instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #dim(instr)# [1] 72 36 1885 #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,] temp=(gridcell3< -99999);gridcell3[temp]=NA index<-length(gridcell3)%%12 gridcell3<-c(gridcell3,rep(NA,12-index) ) close.ncdf(v) read.crutem3<-ts(gridcell3,start=c(1850,1),freq=12) read.crutem3 } ##3. HADSST2 #1850 on #Gridded global SST anomalies from 1856-present derived from UK Met Office SST data #which has had "sophisticated" statistical techniques applied to it to fill in gaps #http://web1.cdc.noaa.gov/cdc/data.noaa.oisst.v2.html #ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc ##ftp://ftp.cdc.noaa.gov/Datasets/kaplan_sst/sst.mon.anom.nc #downloaded from page http://hadobs.metoffice.com/hadsst2/data/download.html #format http://hadobs.metoffice.com/hadsst2/data/Read_instructions_sst.txt #dateline to datelone; S to N #http://hadobs.metoffice.com/hadsst2/data/HadSST2_1850on.nc.gz" #requires library(ncdf) #url="d:/climate/data/gridcell/hadsst2/HadSST2_1850on.nc.gz" #url="http://hadobs.metoffice.com/hadsst2/data/HadSST2_1850on.nc.gz" #example 723A (18 03.079N,57 36.561 E; water depth 807.8 m) #lat= 18+03.079/60;long= 57+36.561/60; read.hadsst2<-function (lat,long,url="d:/climate/data/gridcell/hadsst2/HadSST2_1850on.nc"){ v<-open.ncdf(url) instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 #dim(instr)# [1] 72 36 1885 #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,]; temp=(gridcell3< -99999);gridcell3[temp]=NA index<-length(gridcell3)%%12 gridcell3<-c(gridcell3,rep(NA,12-index) ) close.ncdf(v) read.hadsst2<-ts(gridcell3,start=c(1850,1),freq=12) read.hadsst2 } #cal5annual-1880-2005.nc #4 #HadCRU2 ########## read.hadcru2<-function (lat,long,url="d:climate/data/gridcell/hadcrut2/hadcrut2.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) ## 5. 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) ###6. NOAA ########## #ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/anom/anom-grid2-1880-current.dat.gz ##ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/grid/grid_mean_temp_1880_current.dat.gz #readme: ftp://ftp.ncdc.noaa.gov/pub/data/er-ghcn-sst/er-ghcn-sst.readme.txt anomaly reference not specified #readme: http://www.ncdc.noaa.gov/oa/climate/research/sst/ersstv3.php 1971-2000 reference download.noaa=function(dset="land"){ if(dset=="land") { ### STALE: download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/anom/anom-grid2-1880-current.dat.gz", "anom-grid2-1880-current.dat.gz",mode="wb"); download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/grid/grid_1880_2008.dat.gz", "anom.gz",mode="wb"); # download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/grid/grid_1880_2008.dat.gz", "anom.gz",mode="wb"); # download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/anom/anom-grid2-1880-current.dat.gz", "anom.gz",mode="wb"); # download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/grid/grid_mean_temp_1880_current.dat.gz", "anom.gz",mode="wb"); handle <- gzfile("anom.gz"); Data <- readLines(handle); close(handle); length(Data) # 335916 #12 per line #months on lines 1, 218 #217 lines per month/year # #216*12=2592 gridcells per month # 331142 lines in file N<-length(Data);N; temp<-!is.na(match(1:N,seq(1,N-216,217))) Data<-Data[!temp];N=length(Data);N #334268 noaa=array(NA,dim=c(N,12)) widths0=cbind(seq(1,67,6),seq(6,72,6)) f=function(i) as.numeric(substr(Data,widths0[i,1],widths0[i,2])) for (i in 1:12) noaa[,i]=f(i) temp =(noaa== -32768);noaa[temp]<-NA dim(noaa) #329616 12 # 329616*12/2592 =1526 #sum(!is.na(noaa)) # 776560 noaa= c(t(noaa)) noaa= t( array(noaa,dim=c(2592,length(noaa)/2592 )) );dim(noaa)# 1527 2592 noaa=ts(noaa,start=c(1880,1),freq=12)/100;tsp(noaa) # 1880.000 2008.917 12.000 download.noaa=noaa } if(dset=="ersst") { download.file("ftp://ftp.ncdc.noaa.gov/pub/data/er-ghcn-sst/temp5d6c.v2.adj.1880.last.dat.gz", "d:/climate/data/gridcell/noaa/anom.gz",mode="wb"); handle <- gzfile("d:/climate/data/gridcell/noaa/anom.gz"); Data <- scan(handle);N=length(Data) Data=t( array(Data,dim=c(2594,N/2594)) ) Data=ts(Data[,3:2594],start=c(1880,1),freq=12) temp=(Data== -999.9);Data[temp]=NA download.noaa=Data } download.noaa } read.noaa<-function (lat,long,noaa=noaa,url="d:climate/data/gridcell/noaa/noaa.tab"){ #load(url) k<-jones(lat,long) gridcell3<-noaa[,k] index<-length(gridcell3)%%12 gridcell3<-ts (c(gridcell3,rep(NA,12-index) ),start=c(1880,1),freq=12) rm(noaa) read.noaa<-gridcell3/100 read.noaa } ##station.noaa<-read.noaa(57.5,77.5) ##7. OLD HADCRU (MBH98) ############## ##version used in MBH98 d:climate/data/MBH98.2004/INSTRUMENTAL/anomalies-new.dat read.hadcru.old<-function (lat,long,url="d:/climate/data/MBH98.2004/v.combined.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.harcru.old<-gridcell3/100 read.hadcru.old } #8. NDP020 ######### read.ndp020<-function (lat,long,url="d:climate/data/gridcell/ndp020/v.ndp020.grid.tab"){ load(url) k<-jones.ndp020(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.ndp020<-gridcell3/100 read.ndp020 } ##station.noaa<-read.noaa(57.5,77.5) #9. #KAPLAN SST V2REFEREENCES #DOWNLOAD and save in directory kaplan.data ##library(ncdf) #Spatial coverage: # 87.5S-87.5N (by 5.0), 2.5E-357.5E (by 5.0) #Temporal coverage: # Monthly anomalies from 01/1856 to 08/2005 (for now) kaplan<-function(lat,long){ if(long<0) long<-360+long; kaplan<-cbind( 1+floor(long/5),18-floor(-lat/5)) kaplan} ##kaplan(lat,long) #simple function to give index from lat-long read.kaplansst=function(lat,long) { kaplan.data<-"d:/climate/data/gridcell/kaplan.sst.v2" Kaplan<-open.ncdf( file.path(kaplan.data,"sst.mon.anom.nc")) #ex.nc #kaplan.v2.cdf v1 <- Kaplan$var[[1]] sst <- get.var.ncdf( Kaplan, v1 ) #dim(sst)# 72 36 1808 index=kaplan(lat,long) read.kaplansst= ts( sst[index[1],index[2],],start=c(1856,1),freq=12) read.kaplansst } #http://web1.cdc.noaa.gov/cdc/data.noaa.oisst.v2.html #ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc ##ftp://ftp.cdc.noaa.gov/Datasets/kaplan_sst/sst.mon.anom.nc #10. MSU Gridded download.msu=function(year,dset="t2lt") { url<-file.path("http://www.atmos.uah.edu/data/msu",dset) ##16 values per line 648*16= 10348 values = 4 * 2592 long=seq(-178.75,178.75,2.5);lat=seq(-88.75,88.75,2.5) loc<-file.path(url,paste("tltmonamg.",year,"_5.2",sep="")) fred=readLines(loc) n=nchar(fred) info=fred[n==48];fred=fred[n==80] writeLines(fred,"temp.dat") test=scan("temp.dat",n=-1) # 103680 N=length(test) download.msu= array(test,dim=c(144,72,N/10368))/100 ;##dim(144x72x12); #Eto W; S to N download.msu } #144 72 10 ######################### ######################### ########################## ####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 #v2d - raw; v2adj - adjusted read.ghcn<-function(id0,name0="v2d") { load (file.path( "d:/climate/data/station/ghcn",paste(name0,"tab",sep=".")) ) #checks to see if data set is already loaded temp<-(v2d[,"wmo"]==id0); read.ghcn<-NA 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) read.ghcn<-list(tsite,tsite.anom,m0 )# raw monthly collation and total; names(read.ghcn)<-c("raw","anom","normal") } read.ghcn } #station identification is 5-digit WMO plus 4-digit id read.ghcn.station<-function(id0,dset="raw") { if(dset=="raw") X=v2d; if(dset=="adj") X=v2adj if(class(id0)=="character") id0=as.numeric(id0) read.ghcn.station<-NA if(dset=="raw") temp<-(floor(X[,"id"]/10)==id0); if(dset=="adj") { #if(nchar(id0)==9) temp<- (floor(as.numeric(substr(X[,"id"],4,12)) /10)==id0); #if(nchar(id0)==11) temp<- (floor(as.numeric(substr(X[,"id"],1,12)) /10)==id0); temp=grep(id0,X$id) } K<-length(unique(X[temp,"id"])); if (K>0) { tsite<-NULL for (j in 1:K) { temp1<- (X[temp,"id"] == unique(X[temp,"id"])[j]) tt<-fill.array(X[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<-ts( round(apply(tsite,1,mean,na.rm=TRUE)/10,2),start=c(tsp(tsite)[1],1),freq=12) time0<-unlist(time(tsite)) temp<- (time0>=1961)&(time0<1991);sum(temp) N<-length(tsite)/12 m0<- apply( t( array(tsite,dim=c(12,N))) ,2,mean,na.rm=TRUE) months<-rep(1:12,N) tsite.anom<-round( (tsite-m0[months]),2) read.ghcn.station<-list(tsite,tsite.anom,m0 )# raw monthly collation and total; names(read.ghcn.station)<-c("raw","anom","normal") } read.ghcn.station } # 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.old <- 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 } download_giss_data <- function(station_name, match_num=1,dset=1) { 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<-NA 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") } dimnames(tsite)[[2]]<-id download_giss_data<-tsite download_giss_data } read.giss <- 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]])=c("year","jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec","DJF","MAM","JJA","SON","ANN") temp=(download_data[[j]] == 999.9);download_data[[j]][temp]=NA close(my_handle); } names(download_data)<-id read.giss=download_data read.giss } # station.giss<-download_giss_data(id0) read.giss.station <- function(id0, 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(id0), "&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(urls) == 0 || match_num > length(urls) ) return(NA); temp=( substr(urls,56,63)==id0); url <- paste("http://data.giss.nasa.gov", gsub("^.*\"(/cgi-bin/gistemp/gistemp_station.py[^\"]*)\".*$", "\\1", urls[temp]), sep=""); #varied from Nicholas to leave as vector my_data <- download_html(url); urls.new <- grep("monthly data as text", my_data, value=1); url.new <- gsub("^.*\"(http://data.giss.nasa.gov/work/gistemp/STATIONS/[^\"]*)\".*$", "\\1", urls.new); my_handle <- url(url.new, "rt"); download_data <- read.table(my_handle, skip=1);#instead of NIcholas method names(download_data)[1]<-"year" close(my_handle); temp<-(download_data == 999.9); download_data[temp]<-NA tt<-fill.array(download_data[,1:13]) read.giss.station<-ts(as.matrix(c(t(tt))),start=c(min(as.numeric(row.names(tt))),1),freq=12) read.giss.station } ################################# #### 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,ccode=425,iii="000",url="http://www1.ncdc.noaa.gov/pub/data/ghcn/daily/gsn",option="monthly") { read.gsn<-rep(list(NA),3) names(read.gsn)<-c("raw","anom","m0") loc<-file.path(url,paste(ccode,iii,id0,".dly",sep="")) width0<-c(6,5,4,2,4,rep(c(5,3),31)) test<-try(download.file(loc,"temp.dat") ) if (class(test)=="try-error") read.gsn<-NA else { test=read.fwf("temp.dat",widths=width0,nrow=-1,fill=TRUE) unlink("temp.dat") 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$time<-100*test$year+test$month if (option=="monthly"){ 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) 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 gsn$mincount<-NA; gsn$mincount[temp]<-test$count[temp1] 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 gsn$maxcount<-NA; gsn$maxcount[temp]<-test$count[!temp1] 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" } if (option=="daily"){ 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") x=c(t(outer(test$time[temp1], 1:31,g))) y=c(t( test[temp1,6:36])) tmin=cbind(x,y) tmin=tmin[!is.na(y),] temp2<-(test$type=="TMAX") x=c(t(outer(test$time[temp2], 1:31,g))) y=c(t( test[temp2,6:36])) tmax=cbind(x,y) tmax=tmax[!is.na(y),] time0=sort(unique(c(tmin[,"x"],tmax[,"x"]))) read.gsn=data.frame(time0) read.gsn$tmin=NA;read.gsn$tmax=NA test=match(read.gsn$time0,tmin[,1]) read.gsn$tmin[!is.na(test)]=tmin[test[!is.na(test)],2 ] test=match(read.gsn$time0,tmax[,1]) read.gsn$tmax[!is.na(test)]=tmax[test[!is.na(test)],2 ] } } #try-error read.gsn } ##station.gsn<-read.gsn(id0) ################################# ##JONES 1994 #dowloaded from Phil Jones #10010 710 84 9JAN MAYEN NORWAY 19211990 101921 #10050 780 -142 9ISFJORD RADIO NORWAY 19121979 101912 # [1] " 10010 710 84 9JAN MAYEN NORWAY 19211990 101921" # [2] "1921 -44 -71 -68 -43 -8 22 47 58 27 -20 -21 -40 -13" # [3] "1922 -10 -18 -62 -38 -16 28 47 62 27 -1 -38 -27 -4" # 10010 710 84 9JAN MAYEN NORWA 19911996 101991 # 10250 695 -190 11TROMO/SKATTO NORWA 19911996 101991 # 10280 745 -190 14BJORNOYA NORWA 19911996 101991 read.jones94<-function(id0) { station.id<-read.fwf("d:/climate/data/station/jones94/inv94pdj.txt",widths= c(6) ) station.id[,1]<-as.numeric(levels(station.id[,1]))[station.id[,1]] station.id<-floor(station.id[,1]/10);station.id<-station.id[!is.na(station.id)] if( sum(station.id==id0)==0) read.jones94<-NA else { fred<-readLines("d:/climate/data/station/jones94/cruwlda2.dat") temp<-(substr(fred,1,5)==id0); x<-nchar(fred); #unique(x) #[1] 69 70 45 49 51 #values in 70; stations in other temp<-(x==70);sum(!temp)#288 N<-length(fred);index<-(1:N)[!temp] station.info<-fred[index] M<-length(index); k<-(1:M) [as.numeric(substr(station.info,1,5))==id0] fred<- fred[(index[k]+1):(index[k+1]-1)] x<-cbind(as.numeric(substr(fred,1,4)),as.numeric(substr(fred,5,9)),as.numeric(substr(fred,10,14)), as.numeric(substr(fred,15,19)),as.numeric(substr(fred,20,24)),as.numeric(substr(fred,25,29)), as.numeric(substr(fred,30,34)),as.numeric(substr(fred,35,39)),as.numeric(substr(fred,40,44)), as.numeric(substr(fred,45,49)),as.numeric(substr(fred,50,54)),as.numeric(substr(fred,55,59)), as.numeric(substr(fred,60,64)), as.numeric(substr(fred,65,70))) temp<-(x== -999);x[temp]<-NA } } ############ ##### ##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(as.matrix(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) ###USHCN #ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/hcn_doe_mean_data.Z #filnet="flag" - GISS version removing NAs download.ushcn=function() { library(uncompress) download.file("ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/hcn_doe_mean_data.Z", "temp.Z",mode="wb") handle <- file("temp.Z", "rb") Data <- readBin(handle, "raw", 9999999) close(handle) uncomp_data <- uncompress(Data) download.ushcn<- strsplit(rawToChar(uncomp_data), "\n") download.ushcn=unlist(download.ushcn) download.ushcn } read.ushcn=function(usid,version="new",units="C",filnet_flag="default"){ class(try(ushcn)) hcnid=as.numeric(substr(ushcn,1,6)) temp=(hcnid== usid);print(sum(temp))#207 #554 if(version=="new") widths0=c(6,6,1,1,rep(c(6,4),13) ) if(version=="shap") widths0=c(6,6,1,1,rep(c(6,4),13) ) if(version=="mmts") widths0=c(6,6,1,1,rep(c(6,4),13) ) if(version=="old") widths0=c(6,6,1,1,rep(c(5,4),13) ) if(version=="urban") widths0=c(6,5,1,6,rep(7,12) ) M=length(widths0) what0=rep(FALSE,M);what0[c(1,2,seq(5,29,2))]=TRUE if(version=="urban") what0=c(TRUE,TRUE,FALSE,rep(TRUE,13)) widths0=cumsum(widths0);M=length(widths0) widths0=cbind( 1+ c(0,widths0[1:(M-1)]),widths0) X=array(NA,dim=c(sum(temp),M) ); X=data.frame(X) ff=function(i){ ff=substr(ushcn[temp],widths0[i,1],widths0[i,2]);ff} for(i in 1:M){X[,i]=ff(i)} if(filnet_flag=="flag") flag=X[,seq(6,28,2)] for (i in (1:M)[what0] ){ X[,i]=as.numeric(X[,i]) } if (!(version=="urban")) X=X[,c(1:2,4,seq(5,29,2))] if(version=="old") X[,4:16]=X[,4:16]/100 temp1=(X== -99.99); X[temp1]=NA names(X)=c("id","year","code",month0,"ann") if(version=="shap") X=X[X$code=="A",] if(version=="mmts") X=X[X$code=="+",] if(units=="C") X[,4:16]= round((X[,4:16]-32)*5/9,3) ##X[,16]-apply(X[,4:15],1,mean,na.rm=TRUE)#negligible up to round temp1=(X$code==" ");temp2=(X$code=="+");temp3=(X$code=="A") # blank Raw; + TOBS; A Filnet if(version=="urban") { ushcn.area=fill.array(X[temp1,c(2,4:15)]) read.ushcn=ts(c(t(ushcn.area)),start=c(as.numeric(row.names(ushcn.area)[1]),1),freq=12) } if( (version=="new")|(version=="old")) { ushcn.area=fill.array(X[temp1,c(2,4:15)]) ushcn.tobs=fill.array(X[temp2,c(2,4:15)]) ushcn.filnet=fill.array(X[temp3,c(2,4:15)]) if(sum(!is.na(ushcn.filnet[nrow(ushcn.filnet),]))==0) ushcn.filnet=ushcn.filnet[1:(nrow(ushcn.filnet)-1),] #mary only area=ts(c(t(ushcn.area)),start=c(as.numeric(row.names(ushcn.area)[1]),1),freq=12) tobs=ts(c(t(ushcn.tobs)),start=c(as.numeric(row.names(ushcn.tobs)[1]),1),freq=12) filnet=ts(c(t(ushcn.filnet)),start=c(as.numeric(row.names(ushcn.filnet)[1]),1),freq=12) read.ushcn=ts.union(area,tobs,filnet) dimnames(read.ushcn)[[2]]=c("area","tobs","filnet") if (filnet_flag=="flag") { read.ushcn[is.na(read.ushcn[,2]),"filnet"]=NA } } if (version=="shap") { shap=fill.array(X[,c(2,4:15)]) read.ushcn=ts(c(t(shap)),start=c(as.numeric(row.names(shap)[1]),1),freq=12) } if (version=="mmts") { shap=fill.array(X[,c(2,4:15)]) read.ushcn=ts(c(t(shap)),start=c(as.numeric(row.names(shap)[1]),1),freq=12) } read.ushcn } #ushcn=readLines("d:/climate/data/station/ushcn/urban_mean_fahr.dat") # #id=X$id[11] #foor USHCN sites use USHCN id (6-digit) ##http://www1.ncdc.noaa.gov/pub/data/ghcn/daily/ghcnd-stations.txt ## id="RS000024688" #Ojmyakon read.ghcnd=function(id,method="all"){ if(nchar(id)==5) loc=file.path("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/all",paste("42500",id,".dly",sep=""))#no if(nchar(id)==11) loc=file.path("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/all",paste(id,".dly",sep=""))#no if(method=="ushcn") loc=file.path("ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/all",paste("USC00",id,".dly",sep=""))#no test=try(download.file(loc,"temp.dat")) if(class(test)=="try-error") read.ghcnid=NULL else { fred=readLines("temp.dat") year=as.numeric(substr(fred,12,15));range(year)#1895 2007 month=as.numeric(substr(fred,16,17));range(year)#1895 2007 type0=substr(fred,18,21);unique(type0) #""PRCP" "TMAX" "TMIN" temp2=(type0=="TMAX");sum(temp2) #562 temp3=(type0=="TMIN");sum(temp3) #2526 if ( (sum (temp2)==0)|(sum(temp3)==0)) read.ghcnd=NA else { time0=100*year+month year1<-min(year):max(year) g<-function(x,y){g<-100*x+y;g} chron<-sort(c(outer(year1,1:12,g))); chron=data.frame(chron);names(chron)="time" widths0=c(21,rep(c(5,3),31)) y=cumsum(widths0) widths0=cbind(1+c(0,y[1:(length(y)-1)]),y ) nrow(widths0) #100 widths0=widths0[seq(2,62,2),] #maxiumms x=fred[temp2] #TMAX X=array(NA,dim=c(sum(temp2),31) ); X=data.frame(X) f=function(i){ f=substr(x,widths0[i,1],widths0[i,2]);f} for(i in 1:31){X[,i]=f(i)} for(i in 1:31){X[,i]=as.numeric(X[,i])} temp=(X== -9999);sum(temp);X[temp]=NA count=apply(!is.na(X),1,sum) xtime=time0[temp2] test=match(chron$time,xtime);temp=!is.na(test) chron$count_max=NA;chron$tmax=NA chron[temp,"count_max"]=count;chron[temp,"tmax"]= apply(X[test[temp],],1,mean,na.rm=TRUE) #minimums x=fred[temp3] X=array(NA,dim=c(sum(temp3),31) ); X=data.frame(X) f=function(i){ f=substr(x,widths0[i,1],widths0[i,2]);f} for(i in 1:31){X[,i]=f(i)} for(i in 1:31){X[,i]=as.numeric(X[,i])} temp=(X== -9999);sum(temp);X[temp]=NA count=apply(!is.na(X),1,sum) xtime=time0[temp3] test=match(chron$time,xtime);temp=!is.na(test) chron$count_min=NA;chron$tmin=NA chron[temp,"count_min"]=count;chron[temp,"tmin"]= apply(X[test[temp],],1,mean,na.rm=TRUE) chron[,c("tmax","tmin")]=round(chron[,c("tmax","tmin")]/10,2) #AVERAGES chron$tmean=apply(chron[,c("tmax","tmin")],1,mean) read.ghcnd=ts(chron[,6],start=c(floor(chron[1,1]/100),1),freq=12) } } read.ghcnd } #NORMAIZE anom=function(x,index=1961:1990) { K=floor(tsp(x)[1]) if( is.null(dim(x))) { month= 1+round(12* c(time(x))%%1) x=c( rep(NA,month[1]-1),x ,rep(NA,12- month[length(month)])) test=t(array(x,dim=c(12,length(x)/12)) ) m0=apply(test[index-K+1,],2,mean,na.rm=TRUE) test=scale(test,center=m0,scale=FALSE) anom=round(ts(c(t(test)),start=c(K,1),freq=12),2) anom} else {anom=x; for(i in 1:ncol(x)) {test=t(array(x[,i],dim=c(12,nrow(x)/12)) ) m0=apply(test[index-K+1,],2,mean,na.rm=TRUE) test=round(scale(test,center=m0,scale=FALSE),2) anom[,i]=c(t(test))}} anom } ff=function(x){ff=filter(x,rep(1/24,24));ff} collectmonthly=function(wmo,id0,usid,ccode=425,method="url"){ if(method=="url") {giss.raw=download_giss_data(wmo,dset=0); giss.adj=download_giss_data(wmo,dset=2) } if(method=="collated") {giss.raw=giss.dset0[[paste(wmo)]] ; giss.adj=giss.dset2[[paste(wmo)]]} giss.raw= giss.raw$raw[, paste(ccode,id0,"0",sep="")] fred= try(giss.adj$raw[, paste(ccode,id0,"0",sep="")]) if(class(fred)=="try-error") giss.adj=NA else giss.adj=fred ghcn.raw=read.ghcn(wmo,name0="v2d") ghcn.adj=read.ghcn(wmo,name0="v2adj") ghcn.raw=ghcn.raw[[1]][,paste(id0,0,sep="")] ghcn.adj=ghcn.adj[[1]][,paste(ccode,id0,0,sep="")] #no GHCN ADj if(usid== -999) collectmonthly=ts.union(giss.raw,giss.adj,ghcn.raw,ghcn.adj) else{ test=try(length(ushcn.collation)); if(class(test)=="try-error") load ("d:/climate/data/station/ushcn/ushcn.collation.tab") k=as.character(usid) collectmonthly=ts.union(giss.raw,giss.adj,ghcn.raw,ghcn.adj, ushcn.collation[[1]][,k],ushcn.collation[[2]][,k],ushcn.collation[[3]][,k]) dimnames(collectmonthly)[[2]][5:7]=c("ushcn.area","ushcn.tobs","ushcn.filnet") } collectmonthly } #k is row number in details USHCN # id0=details$ghcnid[k];wmo=floor(id0/1000);title0=as.character(details$name[k]);usid=details$id[k] collect=function(wmo,id0,usid,ccode=425,method0="url"){ X=collectmonthly(wmo,id0,usid,ccode=425,method=method0) test=anom(X) collect=list(X,test,ts(apply(test,2,annavg),start=floor(tsp(test)[1])),details[k,]) names(collect)=c("monthly","anom","annual","details") collect } ff=function(x){ff=filter(x,rep(1/24,24));ff} ###PLOT FUNCTIONS mollweide=function(grid,options="cru") { map("world",proj="mollweide",fill=TRUE,col="white"); lines(mapproject(list(y= -90:90,x=rep(-174.9,181) )),col=4 ) lines(mapproject(list(y= -90:90,x=rep(-174.8,181) )),col=4 ) a=seq(-60,60,30);for(i in 1:5) lines(mapproject(list(y= rep(a[i],3),x= c(-175,0,-174) )),col=4,lty=2 ) a=seq(-150,150,30); for(i in 1:length(a)) lines(mapproject(list(x=rep(a[i],181),y= c(-90:90) )),col=4,lty=2 ) index=(1:nrow(grid))[!is.na(grid$x)] for(i in index){ if(grid$long[i]== -172.5) grid$long[i]=grid$long[i]+.2 #weird little patch for Mollweide left boundary polygon(mapproject(list(x=c(grid$long[i]-2.5,grid$long[i]-2.5,grid$long[i]+2.5,grid$long[i]+2.5), y=c(grid$lat[i]-2.5,grid$lat[i]+2.5,grid$lat[i]+2.5,grid$lat[i]-2.5) )),col=grid$col[i],border=grid$col[i]) } map("world",proj="mollweide",fill=FALSE,add=TRUE); } gridf=function (year,month,method="cru",target=instr) { if (method=="cru") { breaks0=c(-10,-5,-3,-1,-.5,-.2,0,.2,.5,1,3,5,10) col1=c("#0000BB","#0014FF","#006BFF","#00C2FF","lightblue2","lightblue1", "lightgoldenrodyellow", "lightgoldenrod","#FF8700", "#FF5B00","#D70000","#800000") k= (year-1850)*12+ month long=seq(-177.5,177.5,5);lat=seq(-87.5,87.5,5) gridf= cbind( rep(long,36),as.numeric(as.character(gl(36,72,labels=lat))) ,c(target[,,k])) gridf=data.frame(gridf);names(gridf)=c("long","lat","x") #h=function(x) min((1:length(breaks0))[unlist(x["x"])