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()
