##SCRAPING KNMI DATA ### ### knmi.info - information on KNMI moels scraped from webpage in early 2009 ### ### download_html - utility function to read webpages ### ### read.knmi.models - function to read KNMI model data. Requires ### field - "tas" is the field name for temperature. Only one tested so far ### model - use alias nomenclature from knmi.info. ### scenario- use nomenclature from knmi.info. "20c3m" and "sresa1b" are examples. Radio "SCenario Runs" for nomenclature. ### landtype - default is all. Can be "land" or "sea" ### lat - south to north. default c(-90,90) ### long - 0 to 360. default c(0,360) ### version - default is anomaly. "centigrade" gets Centigrade ### Info - requires knmi.info ### ### get.knmi - function to read observation data. Just beginning on this. knmi.info=read.csv("http://www.climateaudit.org/data/models/knmi/knmi.info.csv",sep="\t",header=TRUE) 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); } read.knmi.models=function(field, model,scenario,landtype="",lat=c(-90,90), long=c(0,360), suffix= "&standardunits=true",email=Email,version="anomaly", method="default",Info=knmi.info) { #construct cgi call prefix=paste("http://climexp.knmi.nl/get_index.cgi?email=",email,sep="") suffix=paste(suffix,"&lat1=",lat[1],"&lat2=",lat[2],"&lon1=",long[1],"&lon2=",long[2],sep="") nruns= Info[Info$scenario==scenario&Info$alias==model,"runs"]; #2 url= paste(prefix,"&field=",field,"_",model,"_",scenario,"&",suffix,sep="") if( !(landtype=="")) url=paste(url,"&masktype=",landtype,sep="");url #ping KNMI to retrieve data using above call my_info=download_html(url) #my_info[grep("raw data",my_info)] #option for anomaly data #paste the location of the raw data file (anomaly) if (version=="anomaly") { rawdata_call= my_info[grep("raw data",my_info)][3]; #y n=nchar(rawdata_call) #suffix2=paste(long[1],"-",long[2],"E_",lat[1],"-",lat[2],"N",sep="") if(nruns==1) nruncall="a" else nruncall="__++a" loc=paste("http://climexp.knmi.nl/",substr(rawdata_call,10,n-16) ,sep="") ; #loc Sys.sleep(2) #retrieve the data test=readLines(loc);test[1:5] #extract the model runs in time series form from the data: ncol=1 and ncol>1 index=grep("ensemble",test) # 1 1817 3633 5449 7265 K=length(index); if(K>0) { index=c(index,length(test)+1) writeLines(substr(test,1,24),"temp.dat") test=NULL for (i in 1:K) { working=read.table("temp.dat",skip=index[i]+1,nrow= (index[i+1]-2)-(index[i]+1) ) ; test=ts.union(test,ts(working[,2],start=c(floor(working[1,1]),round( 12* (working[1,1]%%1),0)+1),freq=12)) } dimnames(test)[[2]]=1:K read.knmi.models=test} else { writeLines(substr(test,1,24),"temp.dat") working=read.table("temp.dat",nrow= length(test)-2) ; read.knmi.models= ts(as.matrix(working[,2]),start=c(floor(working[1,1]),round( 12* (working[1,1]%%1),0)+1),freq=12) } if(method=="verbose") read.knmi.models=list(data=read.knmi.models,info=list(field=field,model=as.character(model),scenario=scenario,region=region,loc=loc) ) } #end anomaly version if (version=="centigrade") { latcall=paste(long[1],"-",long[2],"E_",lat[1],"-",lat[2],"N",sep="") working=NULL for (j in 1:nruns) { loc=paste("http://climexp.knmi.nl/data/i",field,"_",model,"_",scenario,"_",latcall,"_","n","_",numbertochar(j-1),".dat",sep="" ) #test=readLines(loc) test= read.fwf(loc,widths=c(5,10,rep(14,11) ),skip=5 ) working=ts.union(working, ts( c(t( test[,2:13])),start=test[1,1],freq=12) ) } dimnames(working)[[2]]=1:nruns read.knmi.models=working } #end non-anomaly version read.knmi.models } numbertochar=function(x) {if(x<10) paste ("0",x,sep="") else paste(x) } ##fields: crutem3, hadcru #long=c(-158,-157);lat=c(21,22);field="crutem3_hadsst2" get.knmi=function(field,long=c(0,360),lat=c(-90,90),email=Email) { prefix= paste("http://climexp.knmi.nl/get_index.cgi?email=",email,sep="") fieldfix=paste("field",field,sep="=") suffix="standardunits=true" latfix= paste("lat1=",lat[1],"&lat2=",lat[2],"&lon1=",long[1],"&lon2=",long[2],sep="") # 0 to 360 longitude #logon logon=download_html( paste("http://climexp.knmi.nl/start.cgi?",email,sep="")) #ping data url=paste(prefix,fieldfix,suffix,latfix,sep="&") my_info=download_html(url) rawdata_call=my_info[grep("raw data",my_info)][3] # my_info[grep("raw data",my_info)][3]; #retrieve data #latcall= paste(long[1],"-",long[2],"E_",lat[1],"-",lat[2],"N",sep="") #loc=paste( paste("http://climexp.knmi.nl/data/i",field,sep=""), latcall,"na.txt",sep="_");loc n=nchar(rawdata_call) loc=paste("http://climexp.knmi.nl/",substr(rawdata_call,10,n-16) ,sep="") ; #loc Sys.sleep(2) #retrieve the data fred=readLines(loc) X=read.fwf(loc,widths=c(10,14)) X[grep("repeat",fred),2]=NA x=ts(rep(NA,10000),start=X[1,1],freq=12) x=window(x,end=max(X[,1],na.rm=T) ) test=match(round(time(x),2),round(X[,1],2));temp=!is.na(test); sum(temp) x[temp]=X[test[temp],2] temp=!is.na(x) get.knmi=round(window(x, end=max(time(x)[temp])),3) #if(method="verbose") get.knmi=list(data=ts( test[,2],start=test[1,1],freq=12),loc=loc) get.knmi } #long=c(-158,-157);lat=c(21,22);field="crutem3_hadsst2" get.knmi.old=function(field,long=c(0,360),lat=c(-90,90),email=Email) { prefix= paste("http://climexp.knmi.nl/get_index.cgi?email=",email,sep="") fieldfix=paste("field",field,sep="=") suffix="standardunits=true" latfix= paste("lat1=",lat[1],"&lat2=",lat[2],"&lon1=",long[1],"&lon2=",long[2],sep="") # 0 to 360 longitude #logon logon=download_html( paste("http://climexp.knmi.nl/start.cgi?",email,sep="")) #ping data url=paste(prefix,fieldfix,suffix,latfix,sep="&") my_info=download_html(url) rawdata_call=my_info[grep("raw data",my_info)][1] # my_info[grep("raw data",my_info)][3]; rawdata_call=strsplit(rawdata_call,",")[[1]][2] #retrieve data #latcall= paste(long[1],"-",long[2],"E_",lat[1],"-",lat[2],"N",sep="") #loc=paste( paste("http://climexp.knmi.nl/data/i",field,sep=""), latcall,"na.txt",sep="_");loc n=nchar(rawdata_call) loc=paste("http://climexp.knmi.nl/",substr(rawdata_call,11,n-14) ,sep="") ; #loc Sys.sleep(2) #retrieve the data fred=readLines(loc) X=read.fwf(loc,widths=c(5,rep(14,12))) temp= (X< -900);sum(temp) X[temp]=NA target=range(X[,1],na.rm=T) target= cbind(target[1]:target[2],array(NA,dim=c(target[2]-target[1]+1,12) )) test=match(target[,1],X[,1]);temp=!is.na(test); sum(temp) target[temp,2:13]=as.matrix(X)[test[temp],2:13] x= ts( c(t(target[,2:13])),start=target[1,1],freq=12) temp=!is.na(x) get.knmi=round(window(x, end=max(time(x)[temp])),3) #if(method="verbose") get.knmi=list(data=ts( test[,2],start=test[1,1],freq=12),loc=loc) get.knmi }