#Wczytanie danych, wstepna eksploracja, przydzielenie wspolrzednych oraz prezentacja danych
setwd('C:/Users/schul_000/Desktop/GEOSTATA KOŁO/A')
dane <- read.csv("data/dane.csv")
surowe_dane <- read.csv("data/dane.csv")
summary(dane)
library('sp')
coordinates(dane) <- ~X+Y
summary(dane)
spplot(dane, "Pb")
spplot(dane, "Zn")
library('ggplot2')
#Zbadanie rozkladu Pb i podstawowych statystyk
ggplot(dane@data, aes(Pb)) + geom_histogram() #prawoskosny, duzo wartosci odstajacych
boxplot(dane$Pb)
summary(dane$Pb)
sd(dane$Pb)
#Zbadanie rozkladu Zn i podstawowych statystyk
ggplot(dane@data, aes(Zn)) + geom_histogram() #prawoskosny, duzo wartosci odstajacych
boxplot(dane$Zn)
summary(dane$Zn) #wartosc NA!!
sd(dane$Zn)
#usuniecie wartosci globalnie odstajacym w jednym i drugim zbiorze
q1 <- as.numeric(quantile(dane$Pb, .25))
q3 <- as.numeric(quantile(dane$Pb, .75))
summary(dane$Pb > q3 + 1.5*(q3-q1))
summary(dane$Pb < q1 - 1.5*(q3-q1))
dane@data[c(3, 92, 29, 74, 31, 102, 20, 27, 25, 23, 30, 26, 10),2] <- NA
ggplot(dane@data, aes(Pb)) + geom_histogram()
boxplot(dane$Pb)
spplot(dane, "Pb")
q1 <- as.numeric(quantile(dane$Zn, .25, na.rm = TRUE))
q3 <- as.numeric(quantile(dane$Zn, .75, na.rm = TRUE))
summary(dane$Zn > q3 + 1.5*(q3-q1))
summary(dane$Zn < q1 - 1.5*(q3-q1))
dane@data[c(92, 29, 102, 99),3] <- NA
ggplot(dane@data, aes(Pb)) + geom_histogram()
boxplot(dane$Pb)
spplot(dane, "Pb")
#stworzenie map semiwariogramów
library('gstat')
dane_pb <- dane[!is.na(dane$Pb), ]
vario_map_pb <- variogram(Pb ~ 1, dane_pb, cutoff = 100, width = 1, map = TRUE)
plot(vario_map_pb, threshold=1)
dane_zn <- dane[!is.na(dane$Zn), ]
vario_map_zn <- variogram(Zn ~ 1, dane_zn, cutoff = 100, width = 1, map = TRUE)
plot(vario_map_zn, threshold=1)
#anizotropia na kacie 45 stopni
vario_kier_pb <- variogram(Pb ~ 1, dane_pb, alpha = c(45, 90, 135, 180))
plot(vario_kier_pb)
vario_kier_zn <- variogram(Zn ~ 1, dane_zn, alpha = c(45, 90, 135, 180))
plot(vario_kier_zn)
#model pb
vario_kier_fit_pb <- vgm(psill = 0.15, model = "Sph", range = 3, nugget = 0.1, anis = c(45, 0.25))
plot(vario_kier_pb, vario_kier_fit_pb, as.table = TRUE)
#model zn
vario_kier_fit_zn <- vgm(psill = 2, model = "Sph", range = 3, nugget = 0.1, anis = c(45, 0.25))
plot(vario_kier_zn, vario_kier_fit_zn, as.table = TRUE)
#siatka
siatka <- expand.grid(X=seq(from=min(surowe_dane$X)-1, to=max(surowe_dane$X)+1, by=0.25),
                      Y=seq(from=min(surowe_dane$Y)-1, to=max(surowe_dane$Y)+1, by=0.25))
coordinates(siatka)<- ~X+Y
gridded(siatka)<-TRUE
plot(siatka)
plot(dane, add = TRUE)
#kriging pb
ok_pb <- krige(Pb~1, dane_pb, siatka, model=vario_kier_fit_pb, maxdist=20)
summary(ok_pb)
spplot(ok_pb, 'var1.pred')
spplot(ok_pb, 'var1.var')
#kriging zn
ok_zn <- krige(Zn~1, dane_zn, siatka, model=vario_kier_fit_zn, maxdist=20)
summary(ok_zn)
spplot(ok_zn, 'var1.pred')
spplot(ok_zn, 'var1.var')
#Ocena
#Pb
#walidacja podzbiorem
#sredni blad estymacji
library('caret')
set.seed(124)
indeks_pb <- as.vector(createDataPartition(dane_pb$Pb, p=0.75, list=FALSE))
train_pb <- dane_pb[indeks_pb, ]
test_pb <- dane_pb[-indeks_pb, ]
vario_train_pb <- variogram(Pb~1, data=train_pb)
ok_train_pb <- krige(Pb~1, train_pb, test_pb, model=vario_kier_fit_pb, maxdist=20)
reszta_ok_pb <-  ok_train_pb$var1.pred - test_pb$Pb
summary(reszta_ok_pb)
ME <- sum(reszta_ok_pb)/length(reszta_ok_pb)
ME
#pierwiastek sredniego bledu kwadratowego
RMSE <- sqrt(sum(reszta_ok_pb^2)/length(reszta_ok_pb))
RMSE
#Zn
#sredni blad estymacji
indeks_zn <- as.vector(createDataPartition(dane_zn$Zn, p=0.75, list=FALSE))
train_zn <- dane_zn[indeks_zn, ]
test_zn <- dane_zn[-indeks_zn, ]
vario_train_zn <- variogram(Zn~1, data=train_zn)
ok_train_zn <- krige(Zn~1, train_zn, test_zn, model=vario_kier_fit_zn, maxdist=20)
reszta_ok_zn <-  ok_train_zn$var1.pred - test_zn$Zn
summary(reszta_ok_zn)
ME <- sum(reszta_ok_zn)/length(reszta_ok_zn)
ME
#pierwiastek sredniego bledu kwadratowego
RMSE <- sqrt(sum(reszta_ok_zn^2)/length(reszta_ok_zn))
RMSE