博文

目前显示的是 七月, 2021的博文

空间插值 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 &

空间插值 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(