空间插值 autoKrige

 library(automap)

library(tidyverse)

library(openxlsx)

library(magrittr)

library(raster)

library(sf)

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'

gansu <- st_read("~/wang_map/甘肃站点.shp")

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

gansu %<>% 

  inner_join(jiangshui,by=c('stationid'='stationid')) %>%

  st_transform(crs=mycrs) %>% 

  as_Spatial() 

gsmap <- st_read("~/wang_map/town_62.shp") %>% 

  st_transform(crs=mycrs) %>% 

  as_Spatial() 

#甘肃轮廓融合

# gsmapline <- rgeos::gUnaryUnion(gsmap)

grd <- makegrid(gsmap, n = 10000)

names(grd) <- c("X", "Y")

gridded(grd)  <- ~ X+Y # Create SpatialPixel object

plot(grd)

# coordinates(grd) <- c("X", "Y")

# gridded(grd) <- TRUE  # Create SpatialPixel object

# fullgrid(grd) 

######### 泛克里金法grd需要指定

grd <- as.data.frame(spsample(gsmap, "regular", n=10000))

names(grd)       <- c("Lon", "Lat")

coordinates(grd) <- c("Lon", "Lat")

gridded(grd)     <- TRUE  # Create SpatialPixel object

fullgrid(grd)    <- TRUE  # Create SpatialGrid object

#普通克里金

kriging_result = automap::autoKrige(X2005.1~1, gansu)

kriging_result = automap::autoKrige(X2005.1~1, gansu,grd)

#泛克里金法

kriging_result = automap::autoKrige(X2005.1~Lon+Lat,gansu,grd)

plot(kriging_result)

# str(kriging_result)

# kriging_result$krige_output$var1.pred

# Extracting parts from the autoKrige object

prediction_spdf = kriging_result$krige_output

sample_variogram = kriging_result$exp_var

variogram_model = kriging_result$var_model

raster_lay <- mask(raster(kriging_result$krige_output),gsmap)

plot(raster_lay)