如何为 R 中的预测创建坐标网格
How to create grid for coordinates for Prediction in R
我正在使用高斯过程模型进行预测，现在我需要根据数据中的坐标使用网格文件，但我没有，我也没有'不知道如何创建它。
我在此 上遵循了 post，但它显示的网格位于宾夕法尼亚州，而不是我的数据坐标所在的芝加哥！
所以我很困惑哪种是创建包含数据中其他列的网格文件的理想方式。
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
数据包括更多的列，例如，小时、天、飞蛾、年、速度，这些观测值是针对 2 个月期间的 10 个位置的，但我只将坐标放在这里以了解如何创建网格.
以下是我按照上述 link 创建网格的步骤：
# 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)
但是当我查看地图时，网格位于宾夕法尼亚地区而不是芝加哥！你知道为什么吗？我选择了我的纬度：41.99795 和我的长：-87.67459，当我把它们放在 Google 地图上时，它显示了芝加哥区域，但在 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")
我不确定您的代码发生了什么，但您的来源似乎设置不正确。我已经更新了上面的代码以在芝加哥上生成一个网格。我从 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)。检查您正在使用的函数和包是否也是这种情况。
编辑：
如何挑选产地？这个特定代码的原点是网格的左下角，所以我去了 Google 地图并选择了一个略微超出城市范围的随机点（基于 Google 的数据），所以在极限的下方和左侧。
如果这些是你的积分
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)
我正在使用高斯过程模型进行预测，现在我需要根据数据中的坐标使用网格文件，但我没有，我也没有'不知道如何创建它。
我在此
所以我很困惑哪种是创建包含数据中其他列的网格文件的理想方式。
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
数据包括更多的列，例如，小时、天、飞蛾、年、速度，这些观测值是针对 2 个月期间的 10 个位置的，但我只将坐标放在这里以了解如何创建网格.
以下是我按照上述 link 创建网格的步骤：
# 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)
但是当我查看地图时，网格位于宾夕法尼亚地区而不是芝加哥！你知道为什么吗？我选择了我的纬度：41.99795 和我的长：-87.67459，当我把它们放在 Google 地图上时，它显示了芝加哥区域，但在 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")
我不确定您的代码发生了什么，但您的来源似乎设置不正确。我已经更新了上面的代码以在芝加哥上生成一个网格。我从 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)。检查您正在使用的函数和包是否也是这种情况。
编辑： 如何挑选产地？这个特定代码的原点是网格的左下角，所以我去了 Google 地图并选择了一个略微超出城市范围的随机点（基于 Google 的数据），所以在极限的下方和左侧。
如果这些是你的积分
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)