Simple kriging playground.

Mój playground z Małopolską.

library(tidyverse) 
library(gstat)
library(sf)

# ---02. Przerabiam dane wysokościowe

dane_plik <- "dane/malopolskie.txt"

dane_nmt <- read_delim(dane_plik, " ", col_names = c("x", "y", "z"), col_types = "ddd") %>% 
  sample_n(size = 1000, replace=F) %>% 
  st_as_sf( coords = c("x", "y"), crs = 2180, agr = "constant", precision = 0.1) %>% 
  st_transform(crs = 4326) %>%
  as_Spatial(.)


#---03. Przerabiam shp na grid

woj_TERYT <- "12"
malopolska_shp <- st_read("dane/województwa.shp") %>% 
  st_transform(crs = 4326) %>%
  filter(jpt_kod_je==woj_TERYT) 


# grid
bbox<- st_bbox(malopolska_shp)

lon <- seq(bbox$xmin, bbox$xmax, length.out = 500)
lat <- seq(bbox$ymin, bbox$ymax, length.out = 500)
grd <- expand.grid(lon = lon, lat = lat) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, agr = "constant") %>%   
  st_join(malopolska_shp, left=FALSE) %>%  
  as_Spatial(.)

malopolska_shp<-malopolska_shp %>% as_Spatial()


#---04. Variogram

#Jak w lm  z~(jakaś zmienna np stężenie w glebie ~ odległość od rzeki)


dt.vgm <- variogram(z~1,dane_nmt)

class(dt.vgm)

dt.fit <-
  fit.variogram(dt.vgm, model = vgm(1,"Lin",600,1)) # fit model

# vgm() list of models

plot(dt.vgm, dt.fit)


### 05. Kriging


lzn.kriged <- krige((z) ~ 1, dane_nmt , grd, model=dt.fit)

lzn.kriged %>% as.data.frame %>% rename(lon=coords.x1, lat=coords.x2) %>% 
  ggplot(aes(x=lon, y=lat)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient2(low="green", mid = "yellow",  high="red",midpoint = 0) +
  theme_bw()