1

假设我们有一个数组modeldata（数据来自陆地模型），其维度为：

> dim(modeldata)
[1] 67420   518

第一个维度包括网格的单元格，第二个维度是时间序列1500:2017
。第一个维度的不寻常长度是由于单独存在陆地单元格以节省空间。

在raster包中，我通过以下方式处理它：

> coords
        [,1]   [,2]
[1,] -179.75 -16.25
[2,] -179.75  65.25
[3,] -179.75  65.75
[4,] -179.75  66.25
[...,]  ...     ...
[67420,] 179.75  71.25

> wgs84 <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
> modeldata_spdf <- sp::SpatialPixelsDataFrame(coords,
                                               data = data.frame(modeldata),
                                               proj4string = wgs84)

> modeldata_brick <- raster::brick(modeldata_spdf)

请不要以这种方式评判我，
我对使用 terra 包的可比较（高性能）方法更感兴趣。
另一种很好的方法是使用SpatRaster遮罩层而不是坐标。

谢谢 ：-）

4

1 回答 1

0

以下方法适用于在已知常规栅格上计算位置子集（即感兴趣像元）的值的情况。对于其他情况values<-，请参阅rast(, xyz=TRUE)和rasterize。

示例数据

library(terra)
r <- rast(res=10)
set.seed(0)
cells <- sample(ncell(r), 5)
xy <- xyFromCell(r, cells)
xy
#        x   y
#[1,] -165 -25
#[2,]   25  55
#[3,] -135 -55
#[4,] -155 -45
#[5,]  -75   5

在这种情况下，我们有 5 个位置，下面的模型数据md, 有 2 个时间步长（层），每个时间步长（层）有 5 个值

md <- cbind(A=1:5, B=5:1)

# compute cell numbers (not needed here as we have them already)
cells <- cellFromXY(r, xy)

# create empty raster
x <- rast(r, nlyr=2)

# assign values
x[cells] <- md

在您的情况下，保留单元格编号而不是坐标可能更有效。

上述内容也应该适用于raster稍作修改。

于 2021-09-01T16:25:50.133 回答