Home › Zpracování dat v GIS › CV11 – Úvod do geostatistiky

CV11 – Úvod do geostatistiky

Zadání cvičení:

  • Vykreslete plochu studované oblasti s vyznačením hodnot sledované veličiny. Interpretujte vývoj hodnot.
  • Proveďte základní statistické posouzení pro normální i logaritmované hodnoty.
  • Vypočítejte základní charakteristiky za předpokladu normálního a lognormálního rozdělení (aritmetický průměr, rozptyl, koef. šikmosti a koef. špičatosti, kvantilové charakteristiky). Popište typ distribuce (jednoduchá nebo smíšená, v případě smíšené odhadnout rozpětí a střední hodnotu dílčích distribucí).
  • Pro normální i logaritmované hodnoty vykreslete histogram a rankitový graf (probability plot). Interpretujte.
  • Vypočtěte hodnoty izotropního semivariogramu (v případě vhodnější LN distribuce použijte logaritmováni hodnot při výpočtu semivariogramu), vykreslete graf semivariogramu, vyberte vhodný teoretický typ, vypočtěte teoretické hodnoty semivariogramu vybraného typu pro příslušné vzdálenosti a vyneste do společného grafu.
  • Analyzujte 4 základní směry (opět celá analýza semivariogramu) za předpokladu geometrické anizotropie. Vyhodnoťte jednotlivé směry.
  • Podle zjištěných hodnot proveďte krigování v oblasti (bodové, základní krigování).

Postup práce v prostředí R

Data ke cvičení: Doubrava.xls
Data v souboru Doubrava.xls, v tomto případě budete pracovat v proměnnou kal.
Data jsou z vrtů realizovaných v odkališti Pilňok u dolu Doubrava jako střední hodnoty příslušných vlastností kalu.
X a Y relativní souřadnice X a Y v metrech,
Pisek obsah písčité frakce (průměr zrn > 0.063 mm) v %,
ALPE obsah prachové frakce (průměr zrn od 0.004 do 0.063 mm) v %,
kal obsah aleuropelitické frakce (průměr zrn pod 0.004 mm) v %,
Ad popelnatost v bezvodém vzorku v %,
W vlhkost v %,
S obsah síry v %.

#Nastavení pracovního adresáře na disk, kde máte uložená data
setwd(“D:/PAD/Geostatistika”)

#Načtení knihovny pro excel a import dat
File – Import Dataset – From Excel

#Zapnutí knihoven pro práci s prostorovými daty (nutné předem instalovat jako Pluginy)
library(ggplot2)
library(scales)
library(moments)

#Vykreslení histogramu hodnot a výpočet základní popisné statistiky.
hist(Doubrava$Kal)
summary(Doubrava$Kal)

#Výpočet koeficientu šikmosti a špičatosti. Koeficient šikmosti popisuje nesymetrii souboru. Pokud se pohybuje okolo 0, jsou hodnoty rozděleny rovnoměrně vpravo i vlevo. Kladná šikmost znamená, že vpravo od průměru se vyskytují odlehlejší hodnoty, většina hodnot je vlevo. Záporná šikmost znamená, že se hodnoty vyskytují většinou vpravo. Koeficient špičatosti pak porovnává dané rozdělení s normálním rozdělením náhodné veličiny. Objektivní koeficient špičatosti je čtvrtým normovaným momentem zmenšeným o hodnotu 3 (nutné odečíst tedy hodnotu 3)

skewness(Doubrava$Kal)
kurtosis(Doubrava$Kal)

#Pro testování normality použijeme Shapiro-Wilkův test normality. Testuje se hypotéza, že náhodný výběr pochází z normálního rozdělení. Pokud je hodnota p-value menší než 0,05, pak zamítáme nulovou hypotézu. Dle výsledku koeficientu šikmosti a testu normality se zamítá nulová hypotéza o normálním rozdělení. Dalším krokem bude testování logaritmovaných hodnot.

shapiro.test(Doubrava$Kal)

#Pro vykreslení kvantil-kvantilového grafu je nutné mít nainstalovanou knihovnu car. Pravděpodobnostní křivka pro dané rozdělení není zcela přímková. Křivka mění svůj směr, což znamená změnu distribuce, z toho lze usoudit, že distribuce je smíšená.

library(car)

qqPlot(Doubrava$Kal)

#Vykreslení plochy studované oblasti. Z výsledného obrázku je patrné, že rozložení hodnot kalové frakce v odkališti není pravidelné. I když se některé hodnoty prolínají je vidět náznak růstu hodnot od severozápadu k jihovýchodu. V blízkosti středu a na severovýchodě je oblast bez hodnot.

Doubrava$cat <- cut(Doubrava$Kal, breaks = c(13,53,60,74,92))hodnota.plot <- ggplot(aes(x = X, y = Y), data = Doubrava)hodnota.plot <- hodnota.plot + geom_point(aes(color = cat))hodnota.plot <- hodnota.plot + coord_equal()hodnota.plot <- hodnota.plot + scale_color_brewer(palette = “YlGnBu”)hodnota.plot  

#Převod hodnot do logaritmické škály a testovány normality.

Doubrava$log <- log(Doubrava$Kal)hist(Doubrava$log)summary(Doubrava$log)shapiro.test(Doubrava$log)qqPlot(Doubrava$log)  

#Vzhledem k tomu, že oba koeficienty šikmosti i špičatosti vycházejí hůře než v případě hodnot bez logaritmizace, v dalších krocích se bude pracovat s původními hodnotami.

#Pro vykreslení semivariačního mraku a samotnou interpolaci je nutné mít nainstalované a zapnuté následující knihovny.

library(sp)library(gstat)

library(lattice)

#Přiřazení prostorových souřadnic.

oblast.sp <- Doubravacoordinates(oblast.sp) <- ~X + Y 

#Vykreslení semivariančního mraku. Geostatistická interpolace vychází z vlastností prostorové veličiny. Tyto vlastnosti lze popsat pomocí variogramu a kovariance. Jelikož tyto nástroje popisují strukturu prostorové veličiny, označují se jako strukturní funkce.
vario.cloud <- variogram(Doubrava$Kal ~ 1, oblast.sp, cloud = TRUE)
plot(vario.cloud)

#Vykreslení semivariogramu v základních směrech 0, 45, 90 a 135°. Pokud je prostorová veličina izotropní, což znamená, že se její vlastnosti nemění v závislosti na směru, je možné pracovat s jednorozměrných semivariogramem.
vario.aniso <- variogram(Doubrava$Kal ~ 1, oblast.sp, alpha = c(0, 45, 90, 135))
plot(vario.aniso)

#Výpočet párů vstupujících do semivariogramu, původní nastavení.

vario <- variogram(Doubrava$Kal ~ 1, data = oblast.sp)
vario

#Výpočet párů dle zadaných parametrů. Otestujte změny v počtu párů spadajících do jednoho kroku změnou parametru width (délka kroku) a cutoff (vzdálenost, do které se počítá variogram. V každé skupině by mělo být alespoň 15 párů, ze kterých se bude počítat teoretický variogram. Se zvolenými parametry pak vykreslete variogram pro všechny směry a zhodnoťte, ve kterém směru dochází k největším změnám.

vario <- variogram(Doubrava$Kal ~ 1, data = oblast.sp, width=10, cutoff=500)
plot(vario)
vario.aniso <- variogram(Doubrava$Kal ~ 1, oblast.sp, alpha = c(0, 45, 90, 135), width=50, cutoff=800)
plot(vario.aniso)

#Zobrazte si funkci vgm a prozkoumejte její parametry. Zkuste definovat parametry semivariogramu – práh, dosah, zbytkový rozptyl. Hodnoty dosazujte tak, aby výsledná křivka, co nejlépe vystihovala tvar variogramu.
?vgm
krivka.sph <- vgm(255, “Gau”, 400, 20)
krivka.sph <- vgm(255, “Gau”, 400, 20)
plot(krivka.sph)
plot(vario, krivka.sph)

#Vytvoření gridu ve směru osy X a Y, do kterého bude probíhat výpočet krigování.

x.range <- as.integer(c(0,900))
y.range <- as.integer(c(0,900))
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1), y=seq(from=y.range[1], to=y.range[2], by=1))
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE
points(Doubrava, pch=1, col=”red”, cex=1)

#Krigování pomocí zvoleného semivariogramu. Vykreslení hodnot pred=predikované hodnoty, a var=rozptyl hodnot.

krig <- krige(formula = Doubrava$Kal ~ 1, locations = oblast.sp, newdata = grd, model = krivka.sph)
color.pal <- colorRampPalette(c(“dark red”,”orange”, “light yellow”))
color.palr <- colorRampPalette(c(“light yellow”,”orange”, “dark red”))
spplot(krig[“var1.pred”], col.regions = color.pal)
spplot(krig[“var1.var”], col.regions = color.palr)