在 R 中为 gstat 中的克里金法创建网格
Create Grid in R for kriging in gstat
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
考虑这些坐标，这些坐标对应于测量降雨数据的气象站。
R 中 gstat 包的介绍使用了 meuse 数据集。在本教程的某个时刻：https://rpubs.com/nabilabd/118172，伙计们在这行代码中使用了 "meuse.grid"：
data("meuse.grid")
我没有这样的文件，我不知道如何创建它，我可以使用这些坐标创建一个吗？或者至少将我指向 material，它讨论了如何为自定义区域创建自定义网格（即不使用来自 GADM 的行政边界）。
可能措辞不对，甚至不知道这个问题对精通R的人是否有意义。仍然，很想听到一些指导，或者至少是提示。非常感谢！
R 的新手总数和统计数据。
编辑：查看我发布的教程中的示例网格，这就是我想要制作的东西。
编辑 2：这种方法可行吗？ https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html
如果您的研究区域是一个多边形，作为 SpatialPolygons 导入，您可以使用 package raster 对其进行栅格化，或者使用 sp::spsample 使用采样类型 [=12] 对其进行采样=].
如果您没有这样的多边形，您可以使用 expand.grid 创建规则分布在矩形 long/lat 区域上的点，使用 seq 生成一系列长和纬度值。
我将分享我为克里金法创建网格的方法。可能有更有效或更优雅的方法来完成相同的任务，但我希望这将成为促进一些讨论的开始。
发帖人想的是每 10 个像素 1 公里，但这可能太多了。我将创建一个像元大小等于 1 公里 * 1 公里的网格。此外，原发布者没有指定网格的原点，所以我会花一些时间来确定一个好的起点。我还假设 Spherical Mercator 投影坐标系是适合投影的选择。这是 Google 地图或开放街道地图的常见投影。
1。加载包
我将使用以下软件包。 sp、rgdal 和 raster 是为空间分析提供许多有用功能的包。 leaflet 和 mapview 是用于快速探索空间数据可视化的包。
# Load packages
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
2。站点位置的探索性可视化
我创建了一个交互式地图来检查四个站点的位置。因为楼主提供了这四个站的经纬度，所以我可以用Latitude/Longitude做一个SpatialPointsDataFrame投影。请注意，Latitude/Longitude 投影的 EPSG 代码是 4326。要了解有关 EPSG 代码的更多信息，请参阅本教程 (https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf)。
# Create a data frame showing the **Latitude/Longitude**
station <- data.frame(lat = c(7.16, 8.6, 8.43, 8.15),
long = c(124.21, 123.35, 124.28, 125.08),
station = 1:4)
# Convert to SpatialPointsDataFrame
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)
mapview 函数将创建交互式地图。我们可以使用这张地图来确定什么是适合网格原点的。
3。确定产地
查看地图后，我确定原点可能在经度 123 和纬度 7 附近。该原点将位于网格的左下角。现在我需要在球面墨卡托投影下找到代表同一点的坐标。
# Set the origin
ori <- SpatialPoints(cbind(123, 7), proj4string = CRS("+init=epsg:4326"))
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))
我首先根据原点的经纬度创建了一个SpatialPoints对象。之后我使用spTransform进行工程改造。对象 ori_t 现在是球面墨卡托投影的原点。请注意，球面墨卡托的 EPSG 代码是 3857.
要查看坐标值，我们可以使用coordinates函数，如下所示。
coordinates(ori_t)
coords.x1 coords.x2
[1,] 13692297 781182.2
4。确定网格的范围
现在我需要确定可以覆盖所有四个点的网格范围以及克里金法所需的区域，这取决于像元大小和像元数量。以下代码根据信息设置范围。我已决定像元大小为 1 km * 1 km，但我需要试验 x 和 y 方向的良好像元数。
# 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))
根据我创建的范围，我可以创建一个数字都等于0的栅格图层。然后我可以再次使用mapview函数来查看光栅和四个站点是否匹配。
# 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)
这个过程我重复了好几次。最后，我决定 x 和 y 方向的单元格数分别为 250 和 200。
5。创建空间网格
现在我已经创建了一个具有适当范围的栅格图层。我可以先将此光栅保存为 GeoTiff 以备将来使用。
# Save the raster layer
writeRaster(ras, filename = "ras.tif", format="GTiff")
最后，要使用包 gstat 中的克里金法函数，我需要将栅格转换为 SpatialPixels。
# Convert to spatial pixel
st_grid <- rasterToPoints(ras, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")
st_grid 是可用于克里金法的 SpatialPixels。
这是一个确定合适网格的迭代过程。在整个过程中，用户可以根据分析需要更改投影、原点、单元格大小或单元格数量。
@yzw 和@Edzer 提出了创建 规则矩形网格 的优点，但有时需要 在定义的多边形，通常用于克里金.
这是一个文档稀少的主题。 One good answer can be found here。我用下面的代码扩展它：
考虑内置的 meuse 数据集。 meuse.grid 是一个不规则形状的网格。 我们如何为我们独特的研究区域制作像 meuse.grid 这样的网格？
library(sp)
data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))
想象一个形状不规则的 SpatialPolygon 或 SpatialPolygonsDataFrame，称为 spdf。您首先在其上构建一个规则的矩形网格，然后通过不规则形状的多边形对规则网格中的点进行子集化。
# First, make a rectangular grid over your `SpatialPolygonsDataFrame`
grd <- makegrid(spdf, n = 100)
colnames(grd) <- c("x", "y")
# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(spdf))
)
# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]
# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts_in))) +
geom_point(aes(x, y))
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
考虑这些坐标，这些坐标对应于测量降雨数据的气象站。
R 中 gstat 包的介绍使用了 meuse 数据集。在本教程的某个时刻：https://rpubs.com/nabilabd/118172，伙计们在这行代码中使用了 "meuse.grid"：
data("meuse.grid")
我没有这样的文件，我不知道如何创建它，我可以使用这些坐标创建一个吗？或者至少将我指向 material，它讨论了如何为自定义区域创建自定义网格（即不使用来自 GADM 的行政边界）。
可能措辞不对，甚至不知道这个问题对精通R的人是否有意义。仍然，很想听到一些指导，或者至少是提示。非常感谢！
R 的新手总数和统计数据。
编辑：查看我发布的教程中的示例网格，这就是我想要制作的东西。
编辑 2：这种方法可行吗？ https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html
如果您的研究区域是一个多边形，作为 SpatialPolygons 导入，您可以使用 package raster 对其进行栅格化，或者使用 sp::spsample 使用采样类型 [=12] 对其进行采样=].
如果您没有这样的多边形，您可以使用 expand.grid 创建规则分布在矩形 long/lat 区域上的点，使用 seq 生成一系列长和纬度值。
我将分享我为克里金法创建网格的方法。可能有更有效或更优雅的方法来完成相同的任务，但我希望这将成为促进一些讨论的开始。
发帖人想的是每 10 个像素 1 公里，但这可能太多了。我将创建一个像元大小等于 1 公里 * 1 公里的网格。此外，原发布者没有指定网格的原点，所以我会花一些时间来确定一个好的起点。我还假设 Spherical Mercator 投影坐标系是适合投影的选择。这是 Google 地图或开放街道地图的常见投影。
1。加载包
我将使用以下软件包。 sp、rgdal 和 raster 是为空间分析提供许多有用功能的包。 leaflet 和 mapview 是用于快速探索空间数据可视化的包。
# Load packages
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
2。站点位置的探索性可视化
我创建了一个交互式地图来检查四个站点的位置。因为楼主提供了这四个站的经纬度，所以我可以用Latitude/Longitude做一个SpatialPointsDataFrame投影。请注意，Latitude/Longitude 投影的 EPSG 代码是 4326。要了解有关 EPSG 代码的更多信息，请参阅本教程 (https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf)。
# Create a data frame showing the **Latitude/Longitude**
station <- data.frame(lat = c(7.16, 8.6, 8.43, 8.15),
long = c(124.21, 123.35, 124.28, 125.08),
station = 1:4)
# Convert to SpatialPointsDataFrame
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)
mapview 函数将创建交互式地图。我们可以使用这张地图来确定什么是适合网格原点的。
3。确定产地
查看地图后，我确定原点可能在经度 123 和纬度 7 附近。该原点将位于网格的左下角。现在我需要在球面墨卡托投影下找到代表同一点的坐标。
# Set the origin
ori <- SpatialPoints(cbind(123, 7), proj4string = CRS("+init=epsg:4326"))
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))
我首先根据原点的经纬度创建了一个SpatialPoints对象。之后我使用spTransform进行工程改造。对象 ori_t 现在是球面墨卡托投影的原点。请注意，球面墨卡托的 EPSG 代码是 3857.
要查看坐标值，我们可以使用coordinates函数，如下所示。
coordinates(ori_t)
coords.x1 coords.x2
[1,] 13692297 781182.2
4。确定网格的范围
现在我需要确定可以覆盖所有四个点的网格范围以及克里金法所需的区域，这取决于像元大小和像元数量。以下代码根据信息设置范围。我已决定像元大小为 1 km * 1 km，但我需要试验 x 和 y 方向的良好像元数。
# 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))
根据我创建的范围，我可以创建一个数字都等于0的栅格图层。然后我可以再次使用mapview函数来查看光栅和四个站点是否匹配。
# 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)
这个过程我重复了好几次。最后，我决定 x 和 y 方向的单元格数分别为 250 和 200。
5。创建空间网格
现在我已经创建了一个具有适当范围的栅格图层。我可以先将此光栅保存为 GeoTiff 以备将来使用。
# Save the raster layer
writeRaster(ras, filename = "ras.tif", format="GTiff")
最后，要使用包 gstat 中的克里金法函数，我需要将栅格转换为 SpatialPixels。
# Convert to spatial pixel
st_grid <- rasterToPoints(ras, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")
st_grid 是可用于克里金法的 SpatialPixels。
这是一个确定合适网格的迭代过程。在整个过程中，用户可以根据分析需要更改投影、原点、单元格大小或单元格数量。
@yzw 和@Edzer 提出了创建 规则矩形网格 的优点，但有时需要 在定义的多边形，通常用于克里金.
这是一个文档稀少的主题。 One good answer can be found here。我用下面的代码扩展它：
考虑内置的 meuse 数据集。 meuse.grid 是一个不规则形状的网格。 我们如何为我们独特的研究区域制作像 meuse.grid 这样的网格？
library(sp)
data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))
想象一个形状不规则的 SpatialPolygon 或 SpatialPolygonsDataFrame，称为 spdf。您首先在其上构建一个规则的矩形网格，然后通过不规则形状的多边形对规则网格中的点进行子集化。
# First, make a rectangular grid over your `SpatialPolygonsDataFrame`
grd <- makegrid(spdf, n = 100)
colnames(grd) <- c("x", "y")
# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(spdf))
)
# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]
# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts_in))) +
geom_point(aes(x, y))