/kriging_test.R

## kriging_test.R
#Edzer Pebesma
#> Institute for Geoinformatics (ifgi), University of Münster
#> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
#> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
#> http://www.52north.org/geostatistics      e.pebesma at wwu.de


#################################################
#krige wil alleen als er geen NA waarden in de data.frame zitten. Heel onhandig
#data(meuse.grid) # only the non-missing valued cells
#coordinates(meuse.grid) = c("x", "y") # promote to SpatialPointsDataFrame
#gridded(meuse.grid) <- TRUE # promote to SpatialPixelsDataFrame
#x = as(meuse.grid, "SpatialGridDataFrame") # creates the full grid  (with NA)
##################################################


require(foreign)
require(RColorBrewer)
require(rgdal)
require(sp)
require(maptools)
require(gstat)
require(raster)

#set projections

#prj <- CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel
#+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs")

prj <- CRS(as.character(NA))      # alles op geen projectie zetten is makkelijker

disc= 'F'
basis.dir = paste (disc, ':/BA8186 NL-D radarcomposiet/04 Projectgegevens/Kriging/' ,sep='')
data.dir = paste (basis.dir, '/data', sep = '')
fig.dir =  paste (basis.dir, '/plaatjes', sep = '')
setwd(data.dir)

grondstations = read.csv('grondstations.csv',  header =T, colClasses = "character")
data.2011 = read.csv ('data_2011.csv', header = T, colClasses = "character")
names(data.2011)[1] = 'date'

# er zijn 428 grondstations maar data.2011 heeft 756 waarden per datum!! ?????
#> dim(grondstations)
#[1] 428  26
#> dim(data.2011)
#[1] 368 756
## [1] "ID"         "EXT_ID"     "NAAM"       "PARENTID"   "OMSCHRIJVI" "REGIO"      "AFKORTING"  "INSTANTIE"  "KWALITEIT"  "RESOLUTIE"  "FREQUENTIE"
#[12] "TIJDZONE"   "X"          "Y"          "COMPLEET"   "TYPE"       "TYPE_DATA"  "KWALITEIT_" "VERBINDING" "WNS1400_1M" "WNS1400_5M" "WNS1400_10"
#[23] "WNS1400_15" "WNS1400_20" "WNS1400_1H" "WNS1400_1D"


year = 2011
month = 1
#for(day in 1:31)
#{
day=13
date = paste (year, '-', formatC(month,width=2,flag='0'), '-', formatC(day,width=2,flag='0'),' ', '08:00:00', sep= '')

data.2011.date = as.numeric(data.2011[ which (data.2011$date == date), ])[-1]
data.2011.date[data.2011.date==-999]=NA

x = as.numeric( grondstations$X )
y = as.numeric( grondstations$Y )
raingauge.data = data.2011.date [328:755]  #!!!!! Dit is gedaan om toch tot het aantal ID's te komen!!
raingauge = na.omit(data.frame( cbind ( x,  y , raingauge.data)) )     #heel belangrijk, anders gaat het straks bij kriging fout!!! (error in dimensions)

raingaugesp = raingauge
coordinates(raingaugesp)=~x+y
projection (raingaugesp) = prj   #zelfde projectie als radarbeeld
pts = list("sp.points", raingaugesp, pch = 3, col = "grey")  #de locatie van de regenmeters


#inlezen radarbeeld
setwd(data.dir)
radar.file = paste ('aggregate.m1_',year, formatC(month,width=2,flag='0'), formatC(day,width=2,flag='0'), '080000.tif',sep='')
radarbeeld <- raster(radar.file)
radarbeeldsp <- as(sampleRegular(radarbeeld, 5e5, asRaster=TRUE), 'SpatialGridDataFrame')
projection (radarbeeldsp) = prj
names (radarbeeldsp) = 'radar.data'
rb = as.data.frame(radarbeeldsp)   #zonder NA waarden, anders gaat het fout bij krigen
coordinates(rb)=c('s1','s2')
gridded(rb) = TRUE

#inverse distance raingauges
raing.id = krige (raingauge.data~1, raingaugesp, radarbeeldsp)
#spplot(raing.id, "var1.pred",  sp.layout = list(pts), 	main = "inverse distance raingauges")


#radarwaarde op plek van raingauge
raing.radar=krige( radar.data~1,rb,raingaugesp,nmax=1)       #zowel raingauge als radar hebben nu geen projectie
raingaugesp$radar.data = raing.radar$var1.pred

#setwd(fig.dir)
#fig=paste('raing_radar_sp',year,formatC(month,width=2,flag='0'),formatC(day,width=2,flag='0'),'.png',sep='')
#png(fig,bg='white')
#sp.theme(TRUE)
#print(spplot( raingaugesp, main = date)  )
#dev.off()

lm.raingauge=lm(raingaugesp$raingauge.data~raingaugesp$radar.data)

setwd(fig.dir)
fig = paste('raing_radar_scatter2',year,formatC(month,width=2,flag='0'),formatC(day,width=2,flag='0'),'.png',sep='')
png(fig,bg='white')
print(xyplot(raingaugesp$raingauge.data~raingaugesp$radar.data,panel=
               function(x,y,...) {
                 panel.xyplot(x,y,...)
                 panel.abline(c(0,1))
                 panel.lmline(x,y,...,lty=2)
               },xlab='radar',ylab='raingauge',
#             xlim=c(min(raingaugesp$raingauge.data,raingaugesp$radar.data),max(raingaugesp$raingauge.data,raingaugesp$radar.data)),
#             ylim=c(min(raingaugesp$raingauge.data,raingaugesp$radar.data),max(raingaugesp$raingauge.data,raingaugesp$radar.data)),
              xlim=c(0,40),
              ylim=c(0,40),
             main=paste('correlation for',date,':r =',round(summary(lm.raingauge)$r.squared^0.5,digits=4))))
dev.off()
#}
#deze grafieken saven

#################
#KRIGING
#################
#standardized pooled variogram of small extent according to paper in Journal of Hydrometeorolgy; Schuurmans et al.

            #small extent
            psill.small=1.270
            model.small='Sph'
            range.small=10000 #in meters
            nugget.small=0.172
            vgm.small.pooled=vgm(psill.small,model.small,range.small,nugget.small)

            #medium extent:
            model.med='Sph'
            nugget.med=0.035
            psill.med.1=0.473
            psill.med.2=8.994
            range.med.1=20000 #in meters
            range.med.2=1040000 #in meters
            vgm.medium.pooled=vgm(psill.med.2,model.med,range.med.2,add.to=vgm(psill.med.1,model.med,range.med.1,nugget.med))

            #large extent:
            model.large='Sph'
            nugget.large=0.053
            psill.large.1=0.281
            psill.large.2=0.795
            range.large.1=23000 #in meters
            range.large.2=156000 #in meters
            vgm.large.pooled=vgm(psill.large.2,model.large,range.large.2,add.to=vgm(psill.large.1,model.large,range.large.1,
            nugget.large))

#variogram van regenmeters: pooled variogram is gestandariseerd dus corrigeren voor variantie van regenmeter data
        vgm.raing=vgm.small.pooled
        vgm.raing$psill=vgm.small.pooled$psill*var(raingaugesp$raingauge.data,na.rm=TRUE)

        v.ok = variogram(raingauge.data~1, raingaugesp)
        #ok.model = fit.variogram(v.ok, vgm(25, "Exp", 120000, 100))
        plot(v.ok,vgm.raing,plot.numbers=T,main=paste('medium extent variogram',date))

    # gstat object voor ordinary kriging
        g.ok=NULL
        g.ok=gstat(g.ok,id='gauge',formula=raingauge.data~1,data=raingaugesp,model=vgm.raing,nmax = 40)

    #predictie naar modelknooppunten (mn) met ordinairy kriging.
        ok.out=predict(g.ok,radarbeeldsp,nsim=0)
        ok.out$gauge.pred[ok.out$gauge.pred<0]=0

#variogram van de residuen
#fit.variogram(object, model, fit.sills = TRUE, fit.ranges = TRUE,
#	fit.method = 7, debug.level = 1, warn.if.neg = FALSE )

v.res=variogram(raingauge.data~radar.data,raingaugesp,cutoff=50000,width=5000) #variogram van de residuen
m.res=fit.variogram(v.res,vgm(0,'Exp',25000,30),fit.ranges=F)
print(plot(v.res,m.res,plot.numbers=T,main=paste('medium extent residual variogram',date)))

#universal kriging
g.ked = NULL
g.ked=gstat(g.ked, id='gauge',raingauge.data~radar.data,
            data=raingaugesp,model=m.res,nmax = 40)
ked.out=predict(g.ked,radarbeeldsp,nsim=0)
ked.out$gauge.pred[ked.out$gauge.pred<0]=0


#plot resultaten bij elkaar
radar = radarbeeldsp
radar[['a']] = radarbeeldsp [["radar.data"]]
radar[['b']] = raing.id [["var1.pred"]]
radar[['c']] = ok.out [['gauge.pred']]
radar[['d']] = ked.out [['gauge.pred']]


png(file = paste("sp_all",date,'.png",sep=''),  bg = "white")
sp.theme(TRUE)
spplot(radar, c("a", "b", 'c', 'd'),
	names.attr = c("radarbeeld", "inverse distance raingauge", 'ordinary kriging', 'KED'),
	as.table = TRUE, main = "radar-raingauge interpolation",
	 at=seq(0,40,by=2),
	sp.layout = list(pts)
)
dev.off()

png(file = "myplot.png", radar.all, bg = "transparent")

png(file="myplot.png", bg="transparent")
plot(1:10)
rect(1, 5, 3, 7, col="white")
dev.off()


##############
#Universal kriging (kriging with external drift)
#############
lm.raindata=lm(raindata$raingauge~raindata$radar)
#summary(lm.raindata)
i='2011 - 01 - 13'
print(xyplot(raindata$raingauge~raindata$radar,panel=
               function(x,y,...) {
                 panel.xyplot(x,y,...)
                 panel.abline(c(0,1))
                 panel.lmline(x,y,...,lty=2)
               },xlab='radar',ylab='raingauge',
             xlim=c(min(raindata$raingauge,raindata$radar),max(raindata$raingauge,raindata$radar)),
             ylim=c(min(raindata$raingauge,raindata$radar),max(raindata$raingauge,raindata$radar)),
             main=paste('correlation for',i,':r =',round(summary(lm.raindata)$r.squared^0.5,digits=4))))


#variogram van de residuen
v.res=variogram(raingauge~radar,raindata,cutoff=50000,width=5000) #variogram van de residuen
m.res=fit.variogram(v.res,vgm(1,'Sph',25000,1))
if (m.st$range[2]==25000) #indien het variogram model van de regenmeters niet eenduidig was, is een range van 25000 meter opgelegd. Dit ook dan doen voor residue variogram
{
  m.res=fit.variogram(v.res,vgm(1,'Sph',150000,1),fit.ranges=FALSE)
}
print(plot(v.res,m.res,plot.numbers=T,main=paste('medium extent residual variogram',i)))


#universal kriging
g.ked = NULL
g.ked=gstat(g.ked, id='raingauge',raingauge~radar,
            data=raindata,model=m.res,nmax = 40)
ked.out=predict(g.ked,radarbeeld,nsim=0)

ked.trend=predict(g.ked,radarbeeld,nsim=0,BLUE=T)
trend.coef=sd(ked.trend$raingauge.pred)/sd(radarbeeld$radar)

gridded(ked.out)=T
print(spplot(ked.out,c('raingauge.pred'),sp.layout = list(pts),cuts=51,main=paste('ked pred sqrt rainfall',i)))
print(spplot(ked.out,c('STNmRaing.var'),sp.layout = list(pts),cuts=51,main=paste('ked var sqrt rainfall',i)))

#cross validation
crval.ked=gstat.cv(g.ked)
filename=paste('crvalSTNl_ked_',i,'.txt',sep='')
write.table(crval.ked,filename,row.names=FALSE)
	#Edzer Pebesma
	#> Institute for Geoinformatics (ifgi), University of Münster
	#> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
	#> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
	#> http://www.52north.org/geostatistics e.pebesma at wwu.de


	#################################################
	#krige wil alleen als er geen NA waarden in de data.frame zitten. Heel onhandig
	#data(meuse.grid) # only the non-missing valued cells
	#coordinates(meuse.grid) = c("x", "y") # promote to SpatialPointsDataFrame
	#gridded(meuse.grid) <- TRUE # promote to SpatialPixelsDataFrame
	#x = as(meuse.grid, "SpatialGridDataFrame") # creates the full grid (with NA)
	##################################################


	require(foreign)
	require(RColorBrewer)
	require(rgdal)
	require(sp)
	require(maptools)
	require(gstat)
	require(raster)

	#set projections

	#prj <- CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel
	#+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs")

	prj <- CRS(as.character(NA)) # alles op geen projectie zetten is makkelijker

	disc= 'F'
	basis.dir = paste (disc, ':/BA8186 NL-D radarcomposiet/04 Projectgegevens/Kriging/' ,sep='')
	data.dir = paste (basis.dir, '/data', sep = '')
	fig.dir = paste (basis.dir, '/plaatjes', sep = '')
	setwd(data.dir)

	grondstations = read.csv('grondstations.csv', header =T, colClasses = "character")
	data.2011 = read.csv ('data_2011.csv', header = T, colClasses = "character")
	names(data.2011)[1] = 'date'

	# er zijn 428 grondstations maar data.2011 heeft 756 waarden per datum!! ?????
	#> dim(grondstations)
	#[1] 428 26
	#> dim(data.2011)
	#[1] 368 756
	## [1] "ID" "EXT_ID" "NAAM" "PARENTID" "OMSCHRIJVI" "REGIO" "AFKORTING" "INSTANTIE" "KWALITEIT" "RESOLUTIE" "FREQUENTIE"
	#[12] "TIJDZONE" "X" "Y" "COMPLEET" "TYPE" "TYPE_DATA" "KWALITEIT_" "VERBINDING" "WNS1400_1M" "WNS1400_5M" "WNS1400_10"
	#[23] "WNS1400_15" "WNS1400_20" "WNS1400_1H" "WNS1400_1D"



	year = 2011
	month = 1
	#for(day in 1:31)
	#{
	day=13
	date = paste (year, '-', formatC(month,width=2,flag='0'), '-', formatC(day,width=2,flag='0'),' ', '08:00:00', sep= '')

	data.2011.date = as.numeric(data.2011[ which (data.2011$date == date), ])[-1]
	data.2011.date[data.2011.date==-999]=NA

	x = as.numeric( grondstations$X )
	y = as.numeric( grondstations$Y )
	raingauge.data = data.2011.date [328:755] #!!!!! Dit is gedaan om toch tot het aantal ID's te komen!!
	raingauge = na.omit(data.frame( cbind ( x, y , raingauge.data)) ) #heel belangrijk, anders gaat het straks bij kriging fout!!! (error in dimensions)

	raingaugesp = raingauge
	coordinates(raingaugesp)=~x+y
	projection (raingaugesp) = prj #zelfde projectie als radarbeeld
	pts = list("sp.points", raingaugesp, pch = 3, col = "grey") #de locatie van de regenmeters


	#inlezen radarbeeld
	setwd(data.dir)
	radar.file = paste ('aggregate.m1_',year, formatC(month,width=2,flag='0'), formatC(day,width=2,flag='0'), '080000.tif',sep='')
	radarbeeld <- raster(radar.file)
	radarbeeldsp <- as(sampleRegular(radarbeeld, 5e5, asRaster=TRUE), 'SpatialGridDataFrame')
	projection (radarbeeldsp) = prj
	names (radarbeeldsp) = 'radar.data'
	rb = as.data.frame(radarbeeldsp) #zonder NA waarden, anders gaat het fout bij krigen
	coordinates(rb)=c('s1','s2')
	gridded(rb) = TRUE

	#inverse distance raingauges
	raing.id = krige (raingauge.data~1, raingaugesp, radarbeeldsp)
	#spplot(raing.id, "var1.pred", sp.layout = list(pts), main = "inverse distance raingauges")


	#radarwaarde op plek van raingauge
	raing.radar=krige( radar.data~1,rb,raingaugesp,nmax=1) #zowel raingauge als radar hebben nu geen projectie
	raingaugesp$radar.data = raing.radar$var1.pred

	#setwd(fig.dir)
	#fig=paste('raing_radar_sp',year,formatC(month,width=2,flag='0'),formatC(day,width=2,flag='0'),'.png',sep='')
	#png(fig,bg='white')
	#sp.theme(TRUE)
	#print(spplot( raingaugesp, main = date) )
	#dev.off()

	lm.raingauge=lm(raingaugesp$raingauge.data~raingaugesp$radar.data)

	setwd(fig.dir)
	fig = paste('raing_radar_scatter2',year,formatC(month,width=2,flag='0'),formatC(day,width=2,flag='0'),'.png',sep='')
	png(fig,bg='white')
	print(xyplot(raingaugesp$raingauge.data~raingaugesp$radar.data,panel=
	function(x,y,...) {
	panel.xyplot(x,y,...)
	panel.abline(c(0,1))
	panel.lmline(x,y,...,lty=2)
	},xlab='radar',ylab='raingauge',
	# xlim=c(min(raingaugesp$raingauge.data,raingaugesp$radar.data),max(raingaugesp$raingauge.data,raingaugesp$radar.data)),
	# ylim=c(min(raingaugesp$raingauge.data,raingaugesp$radar.data),max(raingaugesp$raingauge.data,raingaugesp$radar.data)),
	xlim=c(0,40),
	ylim=c(0,40),
	main=paste('correlation for',date,':r =',round(summary(lm.raingauge)$r.squared^0.5,digits=4))))
	dev.off()
	#}
	#deze grafieken saven

	#################
	#KRIGING
	#################
	#standardized pooled variogram of small extent according to paper in Journal of Hydrometeorolgy; Schuurmans et al.

	#small extent
	psill.small=1.270
	model.small='Sph'
	range.small=10000 #in meters
	nugget.small=0.172
	vgm.small.pooled=vgm(psill.small,model.small,range.small,nugget.small)

	#medium extent:
	model.med='Sph'
	nugget.med=0.035
	psill.med.1=0.473
	psill.med.2=8.994
	range.med.1=20000 #in meters
	range.med.2=1040000 #in meters
	vgm.medium.pooled=vgm(psill.med.2,model.med,range.med.2,add.to=vgm(psill.med.1,model.med,range.med.1,nugget.med))

	#large extent:
	model.large='Sph'
	nugget.large=0.053
	psill.large.1=0.281
	psill.large.2=0.795
	range.large.1=23000 #in meters
	range.large.2=156000 #in meters
	vgm.large.pooled=vgm(psill.large.2,model.large,range.large.2,add.to=vgm(psill.large.1,model.large,range.large.1,
	nugget.large))

	#variogram van regenmeters: pooled variogram is gestandariseerd dus corrigeren voor variantie van regenmeter data
	vgm.raing=vgm.small.pooled
	vgm.raing$psill=vgm.small.pooled$psill*var(raingaugesp$raingauge.data,na.rm=TRUE)

	v.ok = variogram(raingauge.data~1, raingaugesp)
	#ok.model = fit.variogram(v.ok, vgm(25, "Exp", 120000, 100))
	plot(v.ok,vgm.raing,plot.numbers=T,main=paste('medium extent variogram',date))

	# gstat object voor ordinary kriging
	g.ok=NULL
	g.ok=gstat(g.ok,id='gauge',formula=raingauge.data~1,data=raingaugesp,model=vgm.raing,nmax = 40)

	#predictie naar modelknooppunten (mn) met ordinairy kriging.
	ok.out=predict(g.ok,radarbeeldsp,nsim=0)
	ok.out$gauge.pred[ok.out$gauge.pred<0]=0

	#variogram van de residuen
	#fit.variogram(object, model, fit.sills = TRUE, fit.ranges = TRUE,
	# fit.method = 7, debug.level = 1, warn.if.neg = FALSE )

	v.res=variogram(raingauge.data~radar.data,raingaugesp,cutoff=50000,width=5000) #variogram van de residuen
	m.res=fit.variogram(v.res,vgm(0,'Exp',25000,30),fit.ranges=F)
	print(plot(v.res,m.res,plot.numbers=T,main=paste('medium extent residual variogram',date)))

	#universal kriging
	g.ked = NULL
	g.ked=gstat(g.ked, id='gauge',raingauge.data~radar.data,
	data=raingaugesp,model=m.res,nmax = 40)
	ked.out=predict(g.ked,radarbeeldsp,nsim=0)
	ked.out$gauge.pred[ked.out$gauge.pred<0]=0



	#plot resultaten bij elkaar
	radar = radarbeeldsp
	radar[['a']] = radarbeeldsp [["radar.data"]]
	radar[['b']] = raing.id [["var1.pred"]]
	radar[['c']] = ok.out [['gauge.pred']]
	radar[['d']] = ked.out [['gauge.pred']]


	png(file = paste("sp_all",date,'.png",sep=''), bg = "white")
	sp.theme(TRUE)
	spplot(radar, c("a", "b", 'c', 'd'),
	names.attr = c("radarbeeld", "inverse distance raingauge", 'ordinary kriging', 'KED'),
	as.table = TRUE, main = "radar-raingauge interpolation",
	at=seq(0,40,by=2),
	sp.layout = list(pts)
	)
	dev.off()

	png(file = "myplot.png", radar.all, bg = "transparent")

	png(file="myplot.png", bg="transparent")
	plot(1:10)
	rect(1, 5, 3, 7, col="white")
	dev.off()


	##############
	#Universal kriging (kriging with external drift)
	#############
	lm.raindata=lm(raindata$raingauge~raindata$radar)
	#summary(lm.raindata)
	i='2011 - 01 - 13'
	print(xyplot(raindata$raingauge~raindata$radar,panel=
	function(x,y,...) {
	panel.xyplot(x,y,...)
	panel.abline(c(0,1))
	panel.lmline(x,y,...,lty=2)
	},xlab='radar',ylab='raingauge',
	xlim=c(min(raindata$raingauge,raindata$radar),max(raindata$raingauge,raindata$radar)),
	ylim=c(min(raindata$raingauge,raindata$radar),max(raindata$raingauge,raindata$radar)),
	main=paste('correlation for',i,':r =',round(summary(lm.raindata)$r.squared^0.5,digits=4))))


	#variogram van de residuen
	v.res=variogram(raingauge~radar,raindata,cutoff=50000,width=5000) #variogram van de residuen
	m.res=fit.variogram(v.res,vgm(1,'Sph',25000,1))
	if (m.st$range[2]==25000) #indien het variogram model van de regenmeters niet eenduidig was, is een range van 25000 meter opgelegd. Dit ook dan doen voor residue variogram
	{
	m.res=fit.variogram(v.res,vgm(1,'Sph',150000,1),fit.ranges=FALSE)
	}
	print(plot(v.res,m.res,plot.numbers=T,main=paste('medium extent residual variogram',i)))


	#universal kriging
	g.ked = NULL
	g.ked=gstat(g.ked, id='raingauge',raingauge~radar,
	data=raindata,model=m.res,nmax = 40)
	ked.out=predict(g.ked,radarbeeld,nsim=0)

	ked.trend=predict(g.ked,radarbeeld,nsim=0,BLUE=T)
	trend.coef=sd(ked.trend$raingauge.pred)/sd(radarbeeld$radar)

	gridded(ked.out)=T
	print(spplot(ked.out,c('raingauge.pred'),sp.layout = list(pts),cuts=51,main=paste('ked pred sqrt rainfall',i)))
	print(spplot(ked.out,c('STNmRaing.var'),sp.layout = list(pts),cuts=51,main=paste('ked var sqrt rainfall',i)))

	#cross validation
	crval.ked=gstat.cv(g.ked)
	filename=paste('crvalSTNl_ked_',i,'.txt',sep='')
	write.table(crval.ked,filename,row.names=FALSE)