November 27, 2012 10:06
diff --git a/kriging_test.R b/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)