如何在R中对栅格对象进行分层

Question

我当前工作的目标是汇总一个包含 n 个变量的数据库： (1) 位置的纬度，(2) 位置的经度，(3...n-1) 协变量，(n) 地层。 在层变量中，我正在寻找每个点来自的隔间的数量。 但我正在努力获得地层。 我创建了两个栅格 R1 和 R2。

#'Creating the window of domaine D
D <- c(300, 300)  # Square Domaine D of side 300

Win <- owin(xrange =c(0, D[1]), yrange =c(0,D[2])) 

spatstat.options(npixel=c(D[1],D[2]))

# Creating Rasters
library(raster)
# First raster
ext <- extent(Win$xrange, Win$yrange)
v1 <- rnorm(n=seq(from = -10, to=50,length.out = 100), mean = 0, sd=1)
R1 <- raster(matrix(v1,byrow=T,nrow=10, ncol=10))
extent(R1) <- ext
crs(R1) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
names(R1) <- 'temp sim'
R1
# Second Raster
v2 <- rnorm(n=seq(from = 100, to=2000,length.out = 100), mean = 0, sd=1)
R2 <- raster(matrix(v1,byrow=T,nrow=10, ncol=10))
extent(R2) <- ext
crs(R2) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
names(R2) <- 'precip sim'
R2
# Raster List
Rasters.group <- stack(R1, R2)

完成之后，我将栅格转换为 .im 对象

#' Creating .im objects based on the previous rasters. 
#' Each .im object will represent an environmental condition
library(maptools)
IM1 <- as.im.RasterLayer(R1)
IM1 <- as.im.RasterLayer(R2)

#We can plot them as well
library(fields)
par(mfrow=c(2,2))
image.plot(list(x=IM1$xcol, y=IM1$yrow, z=IM1$v), main= "temp sim", asp=1)
image.plot(list(x=IM2$xcol, y=IM2$yrow, z=IM2$v), main= "precip sim", asp=1)

现在我已经组合了 .im 对象，以便从中生成空间点。

# Combination  of invironments

Com <- IM1 + IM1
image.plot(list(x=Com$xcol, y=Com$yrow, z=Com$v), main= "Global env", asp=1, col=c("green", "red", "blue", "yellow", "White", "Black"))

# Generating points Log - Gaussian Cox process

lgc.process <- rLGCP('matern', mu=Com, var=0.5, scale=0.05, nu=1)
lgc.process

xycoor <- cbind(lgc.process$x, lgc.process$y)

lgc.data <- cbind(xycoor, extract(Rasters.group, xycoor))

结果如下：

> head(lgc.data)
                           temp.sim precip.sim
[1,] 194.13086 169.76961 1.82615774 1.82615774
[2,]  31.07026 271.56473 0.09453324 0.09453324
[3,] 109.21386 232.53282 2.30704591 2.30704591
[4,]  27.16373 121.88652 0.83505001 0.83505001
[5,]  97.95982 105.71822 1.85318918 1.85318918
[6,] 167.63728  61.16835 1.14231913 1.14231913

现在它仍然是我希望看到每个点来自哪个兼容性（层）的最后一列。我至少需要10层。 我尝试过的所有东西都不会引导我找到我正在寻找的东西。我使用了split 函数，但没有得到结果。

我该怎么办？

Answer 1

您的代码中有几处看起来很奇怪。例如。当你 创建v1 和v2 你有一个参数n 的特殊向量，但只有 使用了长度，所以相当于：

v1 <- rnorm(100)
v2 <- rnorm(100)

一开始你说的是经度和纬度，但你 应该知道 spatstat 仅适用于平面投影坐标。 很高兴通过sf/maptools/raster/stars/terra如果你 需要使用不同的投影并在它们之间进行转换。在你的 问题中的代码，但是您从不执行任何操作，以下代码是 有点类似于你所做的：

library(spatstat)
Win <- square(300)
IM1 <- as.im(matrix(v1, 10, 10), W = Win)
IM2 <- as.im(matrix(v2, 10, 10), W = Win)
Com <- IM1 + IM1
lgc.process <- rLGCP('matern', mu=Com, var=0.5, scale=0.05, nu=1)
lgc.process # Note this has a lot of points
#> Planar point pattern: 792607 points
#> window: binary image mask
#> 10 x 10 pixel array (ny, nx)
#> enclosing rectangle: [0, 300] x [0, 300] units

将域划分为 nx*ny 个样方（镶嵌）并给它们标签（标记）

nx <- 4
ny <- 3
qtess <- quadrats(Win, nx, ny)
marks(qtess) <- 1:(nx*ny)
qtess
#> Tessellation
#> Tiles are equal rectangles, of dimension 75 x 100 units 
#> 4 by 3 grid of tiles
#> Tessellation is marked
#> window: rectangle = [0, 300] x [0, 300] units
plot(qtess) #' marks/labels not visible on this type of plot

将标记标签作为每个区域的值转换为因子值图像：

qim <- as.im(qtess, dimyx = c(ny, nx))
qim
#> factor-valued pixel image
#> factor levels:
#>  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12"
#> 3 x 4 pixel array (ny, nx)
#> enclosing rectangle: [0, 300] x [0, 300] units
plot(qim, valuesAreColours = FALSE)

最后提取图像在点阵位置的值：

val <- qim[lgc.process]
xycor <- coords(lgc.process)
lgc.data <- cbind(xycor, temp = IM1[lgc.process], precip = IM2[lgc.process], val = val)
head(lgc.data)
#>           x         y      temp     precip val
#> 1 260.93114 167.25626 1.6887990 -0.1590026   8
#> 2 253.33075 242.38444 2.1093809 -0.9992986   4
#> 3  23.57170 244.86029 2.3924167  0.3739248   1
#> 4  14.27467 263.48256 2.3924167  0.3739248   1
#> 5  91.47319  82.87186 1.6035174  1.8282477  10
#> 6 157.93998  26.90952 0.6159271  0.7130826  11