如何为 R 中的预测创建坐标网格

Question

我正在使用高斯过程模型进行预测，现在我需要根据数据中的坐标使用网格文件，但我没有，我也不知道如何创建它。

我关注了link 上的帖子，但它显示的是宾夕法尼亚州的网格，而不是我的数据坐标所在的芝加哥！

所以我很困惑，哪种方法是创建包含数据中其他列的网格文件的理想方式。

station <- data.frame(lat = c(41.997946, 41.960669, 41.960669, 41.960669,41.909269,41.931841,41.909269,41.910561,41.866129,41.866129), long = c(-87.654561, -87.747456, -87.67459, -87.646438,-87.747456,-87.67459,-87.67459,-87.619112,-87.747456,-87.691617),station = 1:10)

 station
            lat      long station
    1  41.99795 -87.65456       1
    2  41.96067 -87.74746       2
    3  41.96067 -87.67459       3
    4  41.96067 -87.64644       4
    5  41.90927 -87.74746       5
    6  41.93184 -87.67459       6
    7  41.90927 -87.67459       7
    8  41.91056 -87.61911       8
    9  41.86613 -87.74746       9
    10 41.86613 -87.69162      10

数据包含更多列，例如 Hour、Day、Moths、Year、Speed，这些观察是针对 2 个月内的 10 个位置的，但我只是将坐标放在这里以了解如何创建网格。

这是我按照上面的链接创建网格的步骤：

# Set the projection. They were latitude and longitude, so use WGS84 long-lat projection
proj4string(station) <- CRS("+init=epsg:4326")

# View the station location using the mapview function
mapview(station)

#3. Determine the origin
# Set the origin
ori_t <- SpatialPoints(cbind(-87.67459, 41.99795), proj4string =  CRS("+init=epsg:4326")) 
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))
coordinates(ori_t)

#ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))
#coordinates(ori_t)


# The origin has been rounded to the nearest 100
x_ori <- round(coordinates(ori_t)[1, 1]/100) * 100
y_ori <- round(coordinates(ori_t)[1, 2]/100) * 100

# Define how many cells for x and y axis
x_cell <- 250
y_cell <- 200

# Define the resolution to be 1000 meters
cell_size <- 1000

# Create the extent
ext <- extent(x_ori, x_ori + (x_cell * cell_size), y_ori, y_ori + (y_cell * cell_size)) 


# Initialize a raster layer
ras <- raster(ext)

# Set the resolution to be
res(ras) <- c(cell_size, cell_size)
ras[] <- 0

# Project the raster
projection(ras) <- CRS("+init=epsg:3857")

# Create interactive map
mapview(station) + mapview(ras)

但是当我查看地图时，网格位于宾夕法尼亚地区而不是芝加哥！你知道为什么吗？我选择了我的 lat : 41.99795 和我的 long:-87.67459 ，当我把它们放在谷歌地图上时，它显示了 Chicago area ，但在 R 上没有显示相同的区域！！

# Convert to spatial pixel
st_grid <- rasterToPoints(ras, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")

另外，当我保存网格文件时，如何使用网格坐标提示其他列？因为它只显示新的 long 和 lat 列

write.csv(st_grid, file = "st_grid.csv")

Answer 1

如果这些是你的观点

library(sp)
station <- data.frame(lat = c(41.997946, 41.960669, 41.960669, 41.960669,41.909269,41.931841,41.909269,41.910561,41.866129,41.866129), long = c(-87.654561, -87.747456, -87.67459, -87.646438,-87.747456,-87.67459,-87.67459,-87.619112,-87.747456,-87.691617),station = 1:10)
coordinates(station) = ~ long+lat
proj4string(station) <- CRS("+proj=longlat +datum=WGS84")
stp <- spTransform(station, CRSobj = CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m"))

你可以的

library(raster)
r <- raster(stp, res=250)

您可以使用extend 或像这样进一步操作它（扩展10公里，然后舍入，不改变分辨率）

rr <- setExtent(r, round(extent(r)+10000,-3), keepres=TRUE)

Answer 2

我不确定您的代码发生了什么，但您的来源似乎设置不正确。我已经更新了上面的代码以在芝加哥生成一个网格。我从 Google 地图中随机选择了一个起点，并修改了 x_cell 和 y_cell 以生成一张大小合理的城市地图。

library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)

station <- data.frame(lat = c(41.997946, 41.960669, 41.960669, 41.960669,41.909269,41.931841,41.909269,41.910561,41.866129,41.866129),
                      long = c(-87.654561, -87.747456, -87.67459, -87.646438,-87.747456,-87.67459,-87.67459,-87.619112,-87.747456,-87.691617),
                      station = 1:10)
coordinates(station) <- ~long + lat
# Set the projection. They were latitude and longitude, so use WGS84 long-lat projection
proj4string(station) <- CRS("+init=epsg:4326")

# View the station location using the mapview function
mapview(station)

#3. Determine the origin
# Set the origin
ori <- SpatialPoints(cbind(-87.872660, 41.619136), proj4string =  CRS("+init=epsg:4326")) 
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))

# The origin has been rounded to the nearest 100
x_ori <- round(coordinates(ori_t)[1, 1]/100) * 100
y_ori <- round(coordinates(ori_t)[1, 2]/100) * 100

# Define how many cells for x and y axis
x_cell <- 60
y_cell <- 80

# Define the resolution to be 1000 meters
cell_size <- 1000

# Create the extent
ext <- extent(x_ori, x_ori + (x_cell * cell_size), y_ori, y_ori + (y_cell * cell_size)) 


# Initialize a raster layer
ras <- raster(ext)

# Set the resolution to be
res(ras) <- c(cell_size, cell_size)
ras[] <- 0

# Project the raster
projection(ras) <- CRS("+init=epsg:3857")

# Create interactive map
mapview(station) + mapview(ras)

这是我最后得到的图像：

至于您的其他问题，我不确定您是否应该将网格与您的数据结合起来。根据您提到的问题中链接的tutorial，例如krige 使用数据meuse 和网格meuse.grid 作为参数：lzn.kriged <- krige(log(zinc) ~ 1, meuse, meuse.grid, model=lzn.fit)。检查您正在使用的函数和包是否也是这种情况。

编辑： 如何选择产地？这个特定代码的原点是网格的左下角，所以我去谷歌地图并选择了一个稍微超出城市范围的随机点（基于谷歌的数据），所以有点低于和有点在左边从极限。