空间插值 sf对象

 library(sf)

library(gstat)

library(openxlsx)

library(magrittr)

library(tidyverse)

mycrs <- '+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +datum=WGS84 +no_defs'

data_jiangshui <- read.xlsx('./Documents/wang_map/降水2005-2009修订后.xlsx')

gs_station<- st_read("./Documents/wang_map/甘肃站点.shp")  %>% 

  st_transform(mycrs)

gs_station %<>% 

  inner_join(data_jiangshui,by=c('stationid'='stationid')) 

gs<- st_read("./Documents/wang_map/town_62.shp") %>% 

  st_transform(mycrs)

gs.border <- st_boundary(gs)

gs.bou <- st_union(gs) %>% 

  st_boundary()

plot(gs.bou)

grid <- st_make_grid(gs.bou, n = c(50, 50)) %>% 

  st_transform(st_crs(gs_station))

plot(grid)

st_crs(gs_station)$proj4string

st_crs(grid)$proj4string

vg.gs <- variogram(gs_station$`2005-7` ~ 1, gs_station)

# vg.gs

# plot(vg.gs$dist, vg.gs$gamma)

# vgm.gs <- vgm(80, "Gau", 500,nugget=200)

fit.gs <-  fit.variogram(vg.gs, vgm(c("Gau", "Sph", "Mat", "Exp")), fit.kappa = TRUE)

plot(vg.gs, fit.gs)

krige.gs <- krige(gs_station$`2005-7` ~ 1, gs_station, grid,model = fit.gs)

krige_gs <- st_intersection(krige.gs, st_union(gs))

ggplot() +

  geom_sf(data = krige_gs, aes(fill = var1.pred), col = NA) +

  geom_sf(data = gs.bou) +

  theme_bw() +

  theme(text = element_text(family = "mono")) +

  scale_fill_continuous(low = "white", high = "red", name = "降水")

ggplot() +

  geom_sf(data = krige_gs, aes(fill = var1.pred), col = NA) +

  geom_sf(data = gs,color='black',size=0.1,alpha=0.09)+

  theme_bw() +

  theme(text = element_text(family = "mono")) +

  scale_fill_continuous(low = "white", high = "red", name = "降水")

# 还可以将插值的结果重新聚合到乡镇街道层次上，这里使用按面积加权的方法进行聚合操作。

street <- st_interpolate_aw(krige_gs, gs, extensive = F)

ggplot()+

  geom_sf(data = street, aes(fill = var1.pred), col = NA) +

  geom_sf(data = gs,color='black',size=0.1,alpha=0.09) +

  scale_fill_continuous(low = "white", high = "red",

                        space = "rgb", name = "lnprice")

######### 泛克里金法

coord_data <- st_coordinates(st_geometry(gs_station))

grid <- st_make_grid(gs, n = c(50, 50)) %>% 

  st_transform(st_crs(gs_station)) %>% 

  st_as_sf()

coord_grid <- st_coordinates(st_geometry(grid)) %>%

  data.frame() %>% 

  unique() %>%

  group_by(L1, L2) %>%

  summarise(X = mean(X), Y = mean(Y))

gs_station$X = coord_data[,1]

gs_station$Y = coord_data[,2]

grid$X = coord_grid$X

grid$Y = coord_grid$Y

vg.gs <- variogram(gs_station$`2005-7` ~ X+Y, gs_station)

fit.gs <-  fit.variogram(vg.gs, vgm(c("Gau", "Sph", "Mat", "Exp")), fit.kappa = TRUE)

plot(vg.gs, fit.gs)

krige.gs <- krige(gs_station$`2005-7` ~ X+Y, gs_station, grid,model = fit.gs)

krige_gs <- st_intersection(krige.gs, st_union(gs))

ggplot() +

  geom_sf(data = krige_gs, aes(fill = var1.pred), col = NA) +

  geom_sf(data = gs.bou) +

  theme_bw() +

  theme(text = element_text(family = "mono")) +

  scale_fill_continuous(low = "white", high = "red", name = "降水")