### Geostatistika - krigovani ### ################################# library(geoR) # Realna data: namerene vysky v 52 mistech ctvercoveho okna o strane 6.7 jednotek # (1 jednotka = 50 stop = 15,24 metru) data(elevation) # vysky v rozmezi 690 az 960 jednotek (1 jednotka = 10 stop = 3,05 metru): points(elevation,cex.min=1,cex.max=4) plot(elevation,lowess=TRUE) # pred vlastnim krigovanim (prostorovou projekci) potrebujeme mit urceny model kovariancni struktury, # odhadneme jako na minulem cviceni pomoci maximalni verohodnosti. # Model s konstantni stredni hodnotou: ml0 <- likfit(elevation, ini=c(3000,2), cov.model="matern", kappa=1.5) ml0 # Model s linearnim trendem: ml1 <- likfit(elevation, trend="1st", ini=c(1300,2), cov.model="matern", kappa=1.5) ml1 ### Jednoduche krigovani ### ############################ # Pripousti nekonstantni stredni hodnotu, musime ji vsak znat (ci odhadnout). # pouzijeme parametry odhadnute metodou maximalni verohodnosti a nalezneme E[Z(x) | Z] (v pripade gaussovskeho pole, resp. obecne NNLO) locs <- pred_grid(c(0,6.3), c(0,6.3), by=0.1) # mrizka, na ktere budeme hledat predikci KC <- krige.control(type="sk", obj.mod=ml0) # jednoduche krigovani pro stacionarni model, nastaveni odhadu sk <- krige.conv(elevation, krige=KC, loc=locs) # vlastni odhad KCt <- krige.control(type="sk", obj.mod=ml1, trend.d = "1st", trend.l="1st") # jednoduche krigovani pro model s linearnim trendem, nastaveni odhadu skt <- krige.conv(elevation, krige=KCt, loc=locs) # vlastni odhad # vykresleni predikovanych hodnot: pred.lim <- range(c(sk$pred, skt$pred)) # jakou skalu hodnot potrebujeme zobrazovat pro predikovane hodnoty pred.lim # vsimneme si - minimum a maximum predikovanych hodnot presne odpovida minimu a maximu pozorovanych hodnot, # tj. predikci nevybocime z intervalu pozorovanych hodnot: summary(elevation) # Dnes ve vsech grafech: tmave barvy = vysoke hodnoty, svetle barvy = nizke hodnoty. windows() image(sk, col=gray(seq(1,0,l=51)), zlim=pred.lim) # model s konst. stredni hodnotou contour(sk, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) windows() image(skt, col=gray(seq(1,0,l=51)), zlim=pred.lim) # model s linearnim trendem contour(skt, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) # vykresleni standardnich chyb predikce (odmocnina z ctvercove chyby predikce) sd.lim <- range(sqrt(c(sk$krige.var, skt$krige.var))) # jakou skalu hodnot potrebujeme zobrazovat pro chybu predikce windows() image(sk, value=sqrt(sk$krige.var), col=gray(seq(1,0,l=51)), zlim=sd.lim) # model s konst. stredni hodnotou contour(sk, value=sqrt(sk$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch="+") windows() image(skt, value=sqrt(skt$krige.var), col=gray(seq(1,0,l=51)), zlim=sd.lim) # model s linearnim trendem contour(skt, value=sqrt(skt$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch="+") # Vetsi rozptyl u modelu s linearnim trendem je zpusoben prave odhadem parametru trendu - vnasi do problemu # dalsi variabilitu. ### Obycejne krigovani ### ########################## # Predpokladame konstantni stredni hodnotu, predikce ma tvar linearni kombinace pozorovanych hodnot, # kde soucet vah = 1 (pripoustime, ze nektere vahy jsou zaporne ci vetsi nez 1). # 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 barev # odhad kovariancni struktury, tentokrat metodou vazenych nejmensich ctvercu: bin1 <- variog(s100,uvec=seq(0,1,l=11)) plot(bin1) wls <- variofit(bin1,ini=c(1,0.5),fix.nugget=T) # fitovani modelu s nulovym nuggetem wls # popis nafitovaneho modelu kovariancni struktury # kde mame data a kde budeme odhadovat: 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, type="ok")) # ordinary kriging, parametry odhadnute pomoci wls, nulovy nugget kc4$predict # predikovane hodnoty kc4$krige.var # chyby predikce # predikce v mistech pozorovani (predikce interpoluje data): kc3 <- krige.conv(s100,locations=s100$coords[1:3,],krige=krige.control(obj.m=wls, type="ok")) s100$data[1:3] kc3$predict # 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=wls, type="ok")) image(kc,loc=pred.grid,col=gray(seq(1,0.1,l=30)),xlab="X",ylab="Y") points(s100,add=T)