在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

Answer 1

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

最初的海报是每10个像素考虑1公里,但这可能太多了.我打算创建一个单元格大小等于1 km*1 km的网格.另外,原始海报没有指明网格的原点,所以我会花一些时间来确定一个好的起点.我还假设球形墨卡托投影坐标系是投影的合适选择.这是Google Map或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 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)

我重复了几次这个过程.最后,我决定细胞的数量是250和200分别用于x和y方向.

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可用于克里金法.

这是确定合适网格的迭代过程.在整个过程中,用户可以根据他们的分析需要更改投影,原点,单元格大小或单元格数量.

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