##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) } 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) #my_info[grep("raw data",my_info)] # my_info[grep("raw data",my_info)][3]; Sys.sleep(4) #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 test=read.fwf(loc,widths=c(10,14)) test=test[!is.na(test[,1]),] get.knmi=ts( test[,2],start=test[1,1],freq=12) #if(method="verbose") get.knmi=list(data=ts( test[,2],start=test[1,1],freq=12),loc=loc) get.knmi }