在R中创建网格以进行gstat中的克里金插值

Question

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")

我没有这样的文件，也不知道如何创建一个，我能用这些坐标创建一个吗？或者至少指点一下我如何为自定义区域创建自定义网格的资料（不使用 GADM 的行政边界）。

可能用词不当，不知道这个问题是否对 R 专家有意义。不过，仍然希望听到一些方向，或者至少是一些提示。非常感谢！

我在 R 和统计学中完全是个新手。

编辑：参见我发布的教程所示的网格示例，那就是我想要制作的东西。

编辑 2：这种方法可行吗？https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html

- ace_01S

---

你想创建的网格的投影和分辨率是多少？ - www

我不确定你的意思，但我猜测分辨率决定了比例尺，所以我希望每个...我不知道，可能是10个像素对应1公里。请原谅我的无知，我还在学习中。 - ace_01S

我在谈论你想要的网格单元大小。投影是将3D坐标（如纬度和经度）转换为2D坐标系统的过程。这份文档可能有助于解释这些概念：https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf。一般来说，我们可以为纬度和经度创建网格，但这可能是没有意义的，因为纬度和经度1度变化的长度不是恒定的。我建议您确定一个适合您需求的投影。 - www

下面，我将展示如何创建一个不规则形状的网格，例如用于克里金插值的meuse.grid。这是一个令人惊讶地鲜有人谈论的话题。 - Rich Pauloo

Answer 1

我将分享我的方法来创建 kriging 网格。可能有更有效或更优雅的方式来完成相同的任务，但我希望这将是一个开始，促进一些讨论。

原始贴主考虑了每 10 个像素为 1 公里，但这可能太多了。我将创建一个单元格大小为 1 公里 * 1 公里 的网格。此外，原始贴主没有指定网格的起点，所以我将花一些时间确定一个好的起始点。我还假设采用 球面墨卡托 投影坐标系是投影的适当选择。这是 Google 地图或 Open Street Maps 的常见投影。

1. 载入包

我将使用以下包。 sp、rgdal 和 raster 是提供许多有用的空间分析功能的包。 leaflet 和 mapview 是用于快速探索可视化空间数据的包。

# Load packages
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)

2. 探索性可视化：站点位置

我创建了一个交互式地图来检查这四个站点的位置。由于原始海报提供了这四个站点的纬度和经度，我可以使用纬度/经度投影创建一个SpatialPointsDataFrame。请注意，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公里*1公里，但需要尝试找出适合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，可用于克里金插值。

这是一个迭代的过程来确定适当的格网。在整个过程中，用户可以根据其分析需求更改投影、起点、单元格大小或单元格数量。

- www

---

你好！感谢您的回复并容忍我的无能。然而，当我尝试执行设置投影线的代码时，出现了以下错误：Geographical CRS given to non-conformant data: 125.08。这是怎么回事？ - ace_01S

我不知道是什么导致了这个错误。根据这里的文档：https://rdrr.io/cran/sp/src/R/Class-Spatial.R. 经度小于-180或大于360时，会出现此错误。但是125.08在此范围之间。 - www

然而，它应该是一个巨大的正方形吗？ - ace_01S

哦，它需要看起来像那样。那就算了吧。该死，这只是显示了我的愚蠢。你正在使用哪个版本的 R？ - ace_01S

@RemyM 感谢您检查我的方法。看起来您的问题值得一个新帖子。如果您不介意，请发布一个带有可重现数据和示例的新问题，以便人们可以查看。 - www

Answer 2

@yzw和@Edzer提出了创建正规矩形网格的好方法，但有时需要在定义的多边形上创建不规则网格，通常用于克里金插值。

这是一个文献记录稀少的话题。可以在这里找到一个很好的答案。我会在下面的代码中进行扩展:

考虑内置的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))

- Rich Pauloo

Answer 3

如果您有一个作为多边形导入的研究区域（使用SpatialPolygons），则可以使用raster软件包将其转换为栅格数据，或使用sp::spsample采用“regular”采样类型对其进行抽样。

如果没有这样的多边形，则可以使用expand.grid在长/纬矩形区域上定期生成点，使用seq生成一系列经度和纬度值。

- Edzer Pebesma