# Pro geostatistickou analyzu existuje v R nekolik knihoven. ###### # geoR ###### library(geoR) # ### Data: # ulozena jako objekt tridy "geodata", obsahuji souradnice poloh ($coords) # a pozorovane hodnoty ($data) v danych polohach, mohou obsahovat dalsi # udaje. Nacitat data ze souboru lze pomoci 'read.geodata'. data(package="geoR") # realna data: data(gambia) # rozsireni malarie u deti ve vesnicich v Gambii plot(gambia.borders,type="l",asp=1) points(gambia[,1:2],pch=19) gambia.map() # zvetseni mapy # simulovana data: data(s100) summary(s100) # zakladni informace o datech plot(s100) # vykresleni dat points(s100) # pouze polohy, kolecka umerna velikosti dat - viz 'points.geodata' points(s100,pt.divide="quintile") # pet ruznych velikosti # ### Variogram: cloud1 <- variog(s100,option="cloud",max.dist=1) # (||x_i-x_j||,(Z(x_i)-Z(x_j))^2/2) cloud2 <- variog(s100,option="cloud",estimator.type = "modulus",max.dist=1) bin1 <- variog(s100,uvec=seq(0,1,l=11)) bin2 <- variog(s100,uvec=seq(0,1,l=11),estimator.type="modulus") par(mfrow=c(2,2)) plot(cloud1,main="klasicky odhad") plot(cloud2,main="robustni odhad") plot(bin1,main="klasicky odhad") plot(bin2,main="robustni odhad") # bin1cl <- variog(s100,uvec=seq(0,1,l=11),bin.cloud=T) bin2cl <- variog(s100,uvec=seq(0,1,l=11),estimator.type="modulus",bin.cloud=T) par(mfrow=c(2,2)) plot(bin1cl,bin.cloud=T,main="klasicky odhad") plot(bin2cl,bin.cloud=T,main="robustni odhad") plot(bin1,main="klasicky odhad") plot(bin2,main="robustni odhad") # par(mfrow=c(1,1)) plot(bin1,type="b") # neparametricky odhad smooth <- variog(s100,option="smooth",max.dist=1,n.points=100,kernel="normal",band=0.2) lines(smooth,type="l",lty=2) # jadrovy odhad lines.variomodel(cov.model="exp",cov.pars=c(1,.3),nugget=0,max.dist=1,lwd=3) # exponencialni model s pevnymi parametry (skutecny model) legend(0.4,0.3,c("empiricky","exponencialni model","vyhlazeny"),lty=c(1,1,2),lwd=c(1,3,1)) # ### Odhady parametru: # 1. "od oka": zkouseni ruznych modelu pomoci 'lines.variomodel' a vizualni # porovnani s empirickym variogramem - 'eyefit' # 2. shoda s empirickym variogramem pomoci metody nejmensich ctvercu: # obycejnych (OLS) nebo vazenych (WLS) - 'variofit' # 3. metody zalozene na verohodnosti: maximalni verohodnost (ML) nebo # rezidualni maximalni verohodnost (REML) pro gaussovska pole - 'likfit' # 4. bayesovske metody - 'krige.bayes' # plot(bin1) lines.variomodel(cov.model="exp",cov.pars=c(1,.3),nug=0,max.dist=1,lwd=3) lines.variomodel(cov.model="mat",cov.pars=c(.85,.2),nug=0.1,kappa=1,max.dist=1,lty=2) lines.variomodel(cov.model="sph",cov.pars=c(.8,.8),nug=0.1,max.dist=1,lwd=2) # eyefit(bin1) # # fitovani modelu s nulovym nuggetem: ml <- likfit(s100,ini=c(1,0.5),fix.nugget=T) # maximalni verohodnost reml <- likfit(s100,ini=c(1,0.5),fix.nugget=T,lik.method="RML") # rezidualni maximalni verohodnost ols <- variofit(bin1,ini=c(1,0.5),fix.nugget=T,weights="equal") # obycejne nejmensi ctverce wls <- variofit(bin1,ini=c(1,0.5),fix.nugget=T) # vazene nejmensi ctverce # fitovani modelu s pevnou hodnotou nuggetu: ml.fn <- likfit(s100,ini=c(1,0.5),fix.nugget=T,nugget=0.15) reml.fn <- likfit(s100,ini=c(1,0.5),fix.nugget=T,nugget=0.15,lik.method="RML") ols.fn <- variofit(bin1,ini=c(1,0.5),fix.nugget=T,nugget=0.15,weights="equal") wls.fn <- variofit(bin1,ini=c(1,0.5),fix.nugget=T,nugget=0.15) # fitovani modelu s odhadnutym nuggetem: ml.n <- likfit(s100,ini=c(1,0.5),nugget=0.5) reml.n <- likfit(s100,ini=c(1,0.5),nugget=0.5,lik.method="RML") ols.n <- variofit(bin1,ini=c(1,0.5),nugget=0.5,weights="equal") wls.n <- variofit(bin1,ini=c(1,0.5),nugget=0.5) # obrazky (nafitovane modely a empiricky variogram): par(mfrow=c(1,3)) plot(bin1,main=expression(paste(tau^2==0))) lines(ml,max.dist=1) lines(reml,lwd=2,max.dist=1) lines(ols,lty=2,max.dist=1) lines(wls,lty=2,lwd=2,max.dist=1) legend(0.55,0.3,legend=c("ML","REML","OLS","WLS"),lty=c(1,1,2,2),lwd=c(1,2,1,2),cex=0.7) plot(bin1,main=expression(paste(tau^2==0.15))) lines(ml.fn,max.dist=1) lines(reml.fn,lwd=2,max.dist=1) lines(ols.fn,lty=2,max.dist=1) lines(wls.fn,lty=2,lwd=2,max.dist=1) legend(0.55,0.3,legend=c("ML","REML","OLS","WLS"),lty=c(1,1,2,2),lwd=c(1,2,1,2),cex=0.7) plot(bin1,main=expression(paste("odhadnute ", tau^2))) lines(ml.n,max.dist=1) lines(reml.n,lwd=2,max.dist=1) lines(ols.n,lty=2,max.dist=1) lines(wls.n,lty=2,lwd=2,max.dist=1) legend(0.55,0.3,legend=c("ML","REML","OLS","WLS"),lty=c(1,1,2,2),lwd=c(1,2,1,2),cex=0.7) # summary: ml.n # beta ... parametr stredni hodnoty # tausq ... nugget (c_0) # sigmasq ... prah-nugget # phi ... parametr rozsahu # dalsi mozne parametry: kappa ... dalsi parametr v modelu variogramu/kovariogramu # lambda ... parametr Boxovy-Coxovy transformace # psiA, psiR ... parametry neizotropnich modelu summary(ml.n) # ### Kriging: par(mfrow=c(1,1)) plot(s100$coords,xlim=c(0,1.2),ylim=c(0,1.2),xlab="X",ylab="Y") loci <- matrix(c(0.2,0.6,0.2,1.1,0.2,0.3,1.0,1.1),ncol=2) text(loci, as.character(1:4), cex=1.3, col="red") polygon(x=c(0,1,1,0),y=c(0,0,1,1),lty=2) # predikce ve ctyrech zvolenych mistech: kc4 <- krige.conv(s100,locations=loci,krige=krige.control(obj.m=wls)) # ordinary kriging, parametry odhadnute pomoci wls, nulovy nugget kc4$predict # predikovane hodnoty kc4$krige.var # chyby predikce # predikce v mistech pozorovani: kc3 <- krige.conv(s100,locations=s100$coords[1:3,],krige=krige.control(obj.m=wls)) s100$data[1:3] kc3$predict # predikce interpoluje data # definovani mrizky: pred.grid <- expand.grid(seq(0,1,l=51),seq(0,1,l=51)) # vypocet predikce na mrizce: kc <- krige.conv(s100,loc=pred.grid,krige=krige.control(obj.m=ml)) # ordinary kriging, odhad pomoci maximalni verohodnosti image(kc,loc=pred.grid,col=gray(seq(1,0.1,l=30)),xlab="X",ylab="Y") points(s100,add=T) # ### Bayesovsky kriging (pro gaussovske modely): pr <- prior.control(phi.discrete=seq(0,2,l=51),phi.prior="reciprocal") # specifikace apriorniho rozdeleni: # phi - diskretni rozdeleni na mrizce 0,0.04, ..., 2 bsp4 <- krige.bayes(s100,loc=loci,prior=pr,output=output.control(n.post=10000)) par(mfrow=c(1,3)) hist(bsp4$posterior$sample$beta,main="",xlab=expression(beta),prob=T) hist(bsp4$posterior$sample$sigmasq,main="",xlab=expression(sigma^2),prob=T) hist(bsp4$posterior$sample$phi,main="",xlab=expression(phi),prob=T) # histogramy s aposteriornim rozdelenim par(mfrow=c(1,1)) plot(bin1,ylim=c(0,1.5)) lines(bsp4,max.dist=1.2,summary=mean) lines(bsp4,max.dist=1.2,summary=median,lty=2) lines(bsp4,max.dist=1.2,summary="mode",post="parameters",lwd=2,lty=2) legend(0.25,0.4,legend=c("variogram posterior mean","variogram posterior median", "parameters posterior mode"),lty=c(1,2,2),lwd=c(1,1,2),cex=0.8) # porovnani empirickeho variogramu s variogramem odhadnutym bayesovsky na zaklade # sumarnich statistik (stredni hodnota, median, modus) aposteriorniho rozdeleni par(mfrow=c(2,2)) for(i in 1:4){ curve(dnorm(x,mean=kc4$pred[i],sd=sqrt(kc4$krige.var[i])),-4,4,lty=2,xlab=paste("Location",i),ylab="density",ylim=c(0,1.1)) lines(density(bsp4$predic$simul[i,])) } # prediktivni hustoty: plna cara ... bayesovsky, # prerusovana cara ... hustota normalniho rozdeleni pred.grid <- expand.grid(seq(0,1,l=51),seq(0,1,l=51)) # bayesovska predikce bsp <- krige.bayes(s100,loc=pred.grid,prior=pr,output=output.control(n.predictive=2)) par(mfrow=c(2,2)) image(bsp,loc=pred.grid,main="predikce",col=gray(seq(1,0.1,l=30))) image(bsp,val="variance",loc=pred.grid,main="predikcni rozptyl",col=gray(seq(1,0.1,l=30))) image(bsp,val="simulation",number.col=1,loc=pred.grid,main="simulace \n z predikcniho rozdeleni",col=gray(seq(1,0.1,l=30))) image(bsp,val="simulation",number.col=2,loc=pred.grid,main="jina simulace \n z predikcniho rozdeleni",col=gray(seq(1,0.1,l=30))) # ### Simulace gaussovskych nahodnych poli sim1 <- grf(100,cov.pars=c(1,.25)) par(mfrow=c(1,2)) points.geodata(sim1,main="simulovana data") plot(sim1,max.dist=1,main="skutecny a empiricky variogram") # sim2 <- grf(441,grid="reg",cov.pars=c(1, .25)) par(mfrow=c(1,2)) image(sim2,main="simulace na hrube mrizi",col=gray(seq(1,.1,l=30))) sim3 <- grf(40401, grid="reg", cov.pars=c(10, .2), met="circ") image(sim3,main="simulace na jemnejsi mrizi",col=gray(seq(1,.1,l=30))) # # Simulovani nahodnych poli je o neco lepsi v balicku 'RandomFields' # ############## # RandomFields ############## library(RandomFields) # PrintModelList() # vsechny implementovane modely model <- "spherical" mean <- 0 variance <- 4 nugget <- 1 scale <- 10 # vyznam parametru - viz help(CovarianceFct) nebo help(Variogram) x <- 0:50 y <- 0:50 f <- GaussRF(x,y,model=model,grid=T,param=c(mean,variance,nugget,scale)) image(x,y,f) # realizace gaussovskeho nahodne pole se sferickou kovariancni funkci cov <- function(h) { CovarianceFct(h,model="spherical",param=c(mean,variance,nugget,scale)) } curve(cov,0,5,n=501,main="Sfericky kovariogram s nuggetem") # # Maternuv-Whittleuv model pro ruzne hodnoty parametru hladkosti: f1 <- GaussRF(x,y,model="whittlematern",grid=T,param=c(0,1,0,10,0.5)) f2 <- GaussRF(x,y,model="whittlematern",grid=T,param=c(0,1,0,10,1.5)) f3 <- GaussRF(x,y,model="whittlematern",grid=T,param=c(0,1,0,10,2.5)) f4 <- GaussRF(x,y,model="whittlematern",grid=T,param=c(0,1,0,10,10)) par(mfrow=c(2,2)) image(x,y,f1) image(x,y,f2) image(x,y,f3) image(x,y,f4) # # interaktivni volba modelu: par(mfrow=c(1,1)) dx <- runif(300,0,8) dy <- runif(300,0,8) dz <- GaussRF(x=dx,y=dy,grid=F,model="gaus",param=c(1,2,1,2)) ev <- EmpiricalVariogram(x=dx,y=dy,data=dz,grid=F,bin=(-1:20)/4) x <- seq(1,5,0.1) ShowModels(x=x,y=x,empirical=ev) # # realna data: data(soil) names(soil) pts <- soil[,1:2] # polohy d <- soil$moisture # data plot(pts) colour <- rainbow(100) plot(pts,col=colour[1+99*(d-min(d))/(max(d)-min(d))],pch=16) ev <- EmpiricalVariogram(pts,data=d,grid=F,bin=c(-1,seq(0,175,l=20))) x <- seq(-150,150,l=121) y <- x ShowModels(0:175,x=x,y=y,empirical=ev) # # Podobne jako v 'geoR' lze odhadovat parametry metodou nejmensich ctvercu # nebo maximalni verohodnosti, slouzi k tomu 'fitvario'. # Prostorova predikce se pak provadi pomoci 'Kriging'. # ####### # gstat ####### library(gstat) # data(package="gstat") data(package="sp") # reka Masa data(meuse.riv) plot(meuse.riv,type="l",asp=1) data(meuse.grid) coordinates(meuse.grid) = c("x", "y") gridded(meuse.grid) = TRUE image(meuse.grid, "dist", add = TRUE) # koncentrace tezkych kovu v pude data(meuse) points(meuse) help(meuse) summary(meuse) plot(y~x,meuse) spplot(meuse[,1:4]) # vice v demo(examples) # ### Existuji i dalsi knihovny (napr. 'spatial' nebo 'sgeostat'), ktere lze pouzit ### pro geostatisticke problemy. # ################## # Besselovy funkce ################## par(mfrow=c(1,1)) nus <- c(0, 1, 3, 5, 10, 20) # prvniho druhu: x <- seq(0, 40, len=801) plot(x,x,ylim=c(-.8,.8),ylab="",type="n",main="Besselovy funkce J_nu(x)") for(nu in c(nus,nus+.5)) lines(x,besselJ(x,nu=nu),col=nu+2) legend(30,-.18,legend=paste("nu=",paste(nus,nus+.5,sep=",")),col=nus+2,lwd=1) # modifikovana druheho druhu: x <- seq(0, 10, len=501) plot(x,x,ylim=c(0,3),ylab="",type="n",main="Besselovy funkce K_nu(x)") for(nu in c(nus,nus+.5)) lines(x,besselK(x,nu=nu),col=nu+2) legend(7.5,1,legend=paste("nu=",paste(nus,nus+.5,sep=",")),col=nus+2,lwd=1) # # chovani u nuly: x0 <- 2^(-20:10) plot(x0,x0^-8,log="xy",ylab="",type="n", main="Besselova funkce J_nu blizko nuly\n log - log meritko") for(nu in c(nus,nus+.5)) lines(x0,besselJ(x0,nu=nu),col=nu+2) legend(3,1e50,legend=paste("nu=",paste(nus,nus+.5,sep=",")),col=nus+2,lwd=1) # plot(x0,x0^-8,log="xy",ylab="",type="n", main="Besselova funkce K_nu blizko 0\n log - log meritko") for(nu in c(nus, nus+.5)) lines(x0,besselK(x0,nu=nu),col=nu+2) legend(3,1e50,legend=paste("nu=",paste(nus,nus+.5,sep=",")),col=nus+2,lwd=1) # # zakladni funkce ze spektralniho rozkladu izotropni kovariancni funkce: omega <- function(d,t) { nu <- d/2-1 (2/t)^nu*gamma(d/2)*besselJ(t,nu) } omega1 <- function(x) sapply(x,omega,d=1) # priprava pro vykresleni omega2 <- function(x) sapply(x,omega,d=2) omega3 <- function(x) sapply(x,omega,d=3) omegainfty <- function(x) exp(-x^2) plot(c(0,35),c(-1,1),type="n") curve(omega1,0,30,col="green",add=T) curve(omega2,0,30,col="blue",add=T) curve(omega3,0,30,col="grey",add=T) curve(omegainfty,0,30,col="brown",add=T) legend(28,1,legend=paste("d=", paste(c(1,2,3,"infty"))),col=c("green","blue","grey","brown"),lwd=1) # q() # koncime