克里金预测函数为每个值返回 NAN [r]

如何解决克里金预测函数为每个值返回 NAN [r]

我正在尝试对我的数据执行克里金法预测。

我遵循一个非常标准的流程。 首先，我准备了一个应该插入克里金法预测的网格。 我准备要插入的数据。 我确保两个对象的 CRS 相同。然后作为下一步，选择变异函数。最后，我们使用克里金法对值进行插值。

然而，由于某种原因，返回的克里金法对象包含所有 NAN 预测。

我尝试了不同的 CRS 系统，但无法想出解决此问题的方法。

因此，我想征求任何建议。

示例：

library(tidyverse)
library(sf)
library(raster)
library(stars)
library(gstat)
library(ggmap)

# get data. using your example,we'll take spatial data from the raster package
Regions <- getData("GADM",country = "CZ",level = 1)

Regions <- Regions[Regions$NAME_1 == 'Prague',]

Regions <- 
  Regions %>% 
  st_as_sf()

grid_spacing <- 0.005

# Create a grid for the border --------------------------------------------

polygony <- 
  st_make_grid(Regions,square = T,cellsize = c(grid_spacing,grid_spacing)) %>%
  st_sf()

grid <- st_intersection(polygony,Regions)

plot(grid)

grid2 <- as_Spatial(grid$geometry)
plot(grid2)

coordinates(grid2)

spdf <- grid2
grd <- makegrid(spdf,n = 3000,pretty = F)

colnames(grd) <- c('x','y')

grd %>% plot()

grd_pts <- SpatialPoints(coords = grd,proj4string=CRS(proj4string(spdf)))

# find all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf,]
grd_pts_in %>% plot()

# transform grd_pts_in back into a data frame
gdf <- as.data.frame(coordinates(grd_pts_in)) 

gdf %>% 
  mutate(value = rnorm(nrow(gdf))) %>% 
  ggplot(aes(x,y,fill =  value)) + 
  geom_tile()

grid2 %>% plot()
gdf %>% plot()  

system <- "+proj=utm +zone=19 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0"
system2 <- "+proj=longlat +datum=WGS84 +no_defs "



# Data to be interpolated -------------------------------------------------

df <- 
  tibble(
long = sample(coordinates(grid2)[,1],500),lat = sample(coordinates(grid2)[,2],val = rnorm(500,5000,10)
)

df %>% 
  ggplot(aes(long,lat,color = val)) + 
  geom_point()

# Kriging Predictions Itself ----------------------------------------------

resid2 <- df %>% 
  st_as_sf(coords = c("long","lat"),crs = system2)


cutoff = 3
lzn.vgm <- variogram(val ~ 1,data = resid2,cutoff = cutoff)
plot(lzn.vgm)

nugget <- 100 # is the value of the semivariance at zero distance.
psill <- max(lzn.vgm$gamma) - nugget # is the difference between the sill and the nugget
range <- 0.2 # distance where the model first begins to flattens out.

lzn.fit.Ste <- fit.variogram(lzn.vgm,model = vgm(
  model = "Ste",nugget = nugget,sill = psill,range = range))

plot(lzn.vgm,lzn.fit.Ste)


coordinates(grid2)

POLY2 <- as(grid2,"Spatial")
RES <- as(resid2,"Spatial")

crs(RES) = CRS(system2)
crs(grid2)  = CRS(system2)

st_crs(grid2) == st_crs(RES)

OK =
  krige(val ~ 1,RES,grid2,model = lzn.fit.Ste)

OK %>% 
  as.data.frame()