在R中创建栅格以在GSTATE中进行克里金法 [英] 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")

我没有这样的文件，我不知道如何创建它，我可以使用这些坐标创建一个吗？或者至少向我推荐讨论如何为自定义区域创建自定义网格(即不使用GADM的行政边界)的材料。

可能用错了这个词，甚至不知道这个问题对精明的人来说是否有意义。尽管如此，我还是希望听到一些方向，或者至少是一些建议。非常感谢！

R和统计的新手总数。

编辑：查看我发布的教程的样例网格，这就是我想要创建的内容。

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

推荐答案

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

最初的海报认为每10个像素需要1公里，但这可能太多了。我要创建一个网格，网格大小等于1千米*1千米。另外，最初的海报没有指定网格的原点，所以我会花一些时间来确定一个很好的起点。我还假设球墨卡托投影坐标系是投影的合适选择。这是Google地图或Open Street地图的常见投影。

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。

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

这篇关于在R中创建栅格以在GSTATE中进行克里金法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！