#Wczytanie danych, wstepna eksploracja, przydzielenie wspolrzednych oraz prezentacja danych
setwd('C:/Users/schul_000/Desktop/GEOSTATA KOŁO')
dane <- read.csv("data/dane.csv")
surowe_dane <- dane
summary(dane)
library('sp')
coordinates(dane) <- ~X+Y
summary(dane)
#Zbadanie typu probkowania
spplot(dane, "Fish")
library('ggplot2')
#Zbadanie rozkladu
ggplot(dane@data, aes(Fish)) + geom_histogram() #prawoskosny, duzo zer, wystepuja bardzo wysokie wartosci
boxplot(dane$Fish)
summary(dane$Fish)
sd(dane$Fish)
#ze wzgledu na typ danych (ryby sa lub ich nie ma) nie usuwamy wartosci globalnie i lokalnie
library('gstat')
#tworzymy mape semiwariogramu
vario_map <- variogram(Fish ~ 1, dane, cutoff = 1000, width = 30, map = TRUE)
plot(vario_map, threshold=10)#nie wystepuje anizotropia
#Semiwariancja
semivario <- variogram(Fish~1, dane)
plot(semivario)
#Model
model <- vgm(psill=300, model = 'Sph', range = 20, add.to = vgm(850, "Nug"))
model2 <- fit.variogram(semivario,vgm(psill=300, model = 'Sph', range = 20, add.to = vgm(850, "Nug")))
plot(semivario, model)
plot(semivario, model2)
#Siatka
siatka <- expand.grid(X=seq(from=min(surowe_dane$X)-2, to=max(surowe_dane$X)+2, by=2),
                      Y=seq(from=min(surowe_dane$Y)-2, to=max(surowe_dane$Y)+2, by=2))
coordinates(siatka)<- ~X+Y
gridded(siatka)<-TRUE
plot(siatka)
plot(dane, add = TRUE)
#kriging zwykły, bo prosty wymaga lokalnej sredniej
ok <- krige(Fish~1, dane, siatka, model=model, maxdist=30)
summary(ok)
#wynik estymacji i błąd
spplot(ok, 'var1.pred')
spplot(ok, 'var1.var')
#Ocena
#sredni blad estymacji
library('caret')
set.seed(124)
indeks <- as.vector(createDataPartition(dane$Fish, p=0.75, list=FALSE))
indeks
train <- dane[indeks, ]
test <- dane[-indeks, ]
vario_train <- variogram(Fish~1, data=train)
ok_train <- krige(Fish~1, train, test, model=model, maxdist=30) #nmax
reszta_ok <-  ok_train$var1.pred - test$Fish
summary(reszta_ok)
ME <- sum(reszta_ok)/length(reszta_ok)
ME
#pierwiastek sredniego bledu kwadratowego
RMSE <- sqrt(sum(reszta_ok^2)/length(reszta_ok))
RMSE
library('raster')
estymacja <- raster(ok['var1.pred'])
writeRaster(estymacja, filename="Nazwisko_Imie.tif", overwrite=TRUE)