使用来自其他两个已批准栅格的重叠批准栅格进行分区操作

Question

我有以下重新分类和批准的栅格，我试图截取/叠加/重叠以获得新的栅格。想法是从栅格 R1 和 R2 的截取区域中获取新的叠加栅格。完成此操作后，我将进行区域操作。 Here R1、R2、ED 栅格。

R1:
class      : RasterLayer 
dimensions : 1399, 1855, 2595145  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : UK_GDP_2010_PPP_percapita_km2 
values     : 1, 6  (min, max)
attributes :
       ID AG
 from:  1  a
  to :  6  f

R2:
class      : RasterLayer 
dimensions : 1399, 1855, 2595145  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 1, 5  (min, max)
attributes :
 ID AG
  1  A
  2  B
  3  C
  4  D
  5  E

这里是拦截/叠加/重叠的代码

1st approach
library(sf)
library(tidyverse)
R1_SPDF <- as(R1,'SpatialPolygonsDataFrame')
R1_SPDF <- st_as_sf(R1_SPDF)

R2_SPDF <- as(R2,'SpatialPolygonsDataFrame')
R2_SPDF <- st_as_sf(R2_SPDF)

R3 <- st_intersection(R1_SPDF, R2_SPDF)

R3:
Simple feature collection with 1174501 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -8.166667 ymin: 50.01667 xmax: 1.541667 ymax: 59.55
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
    UK_GDP_2010_PPP_percapita_km2 layer
1                               1     1
2                               1     1
2.1                             1     1
3                               1     1
3.1                             1     1
4                               1     1
4.1                             1     1
1.1                             1     1
5                               1     1
6                               1     1
                      geometry
1   POLYGON ((-1.641667 59.55, ...
2   POLYGON ((-1.633333 59.55, ...
2.1 POLYGON ((-1.633333 59.55, ...
3   POLYGON ((-1.625 59.55, -1....
3.1 POLYGON ((-1.625 59.55, -1....
4   POLYGON ((-1.616667 59.55, ...
4.1 POLYGON ((-1.616667 59.55, ...
1.1 POLYGON ((-1.641667 59.5416...
5   POLYGON ((-1.65 59.54167, -...
6   POLYGON ((-1.641667 59.5416...

但是，我不确定这个结果是否是我正在寻找的结果，因为我希望新的批准栅格 R3 具有由 R1 和 R2 区域的组合/截取形成的覆盖区域（R3 区域批准：Aa，. .Af,...,Ea,...Ef) 或类似的东西。

Expected R3:
class      : RasterLayer 
dimensions : #from intersection 
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -13.69167, 1.766667, 49.86667, 61.525  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : layer 
values     : 1, 30  (min, max) #aproximately 30 because R1: ID=6, and R2: ID=5.
attributes :
 ID AGnew
  1  Aa
  2  Ab
  .  .
  .  .
  .  .
  30 Ef

这里再次尝试使用光栅包：

2nd approach
library(raster)
R3.1 <- intersect(R1_SPDF, R2_SPDF)
R3.1: 
Simple feature collection with 485611 features and 0 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -8.65 ymin: 49.875 xmax: 1.766667 ymax: 60.84167
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
                     geometry
1  POLYGON ((-0.9 60.84167, -0...
2  POLYGON ((-0.8916667 60.841...
3  POLYGON ((-0.8833333 60.841...
4  POLYGON ((-0.875 60.84167, ...
5  POLYGON ((-0.85 60.84167, -...
6  POLYGON ((-0.8416667 60.841...
7  POLYGON ((-0.9 60.83333, -0...
8  POLYGON ((-0.8916667 60.833...
9  POLYGON ((-0.8833333 60.833...
10 POLYGON ((-0.875 60.83333, ...

获得 R3 后，我希望执行以下区域操作。对 R3 覆盖区域内的 ED 栅格值求和。

sum_R <- zonal(ED, R3, "sum")

非常欢迎任何建议。

Answer 1

我认为答案根本不清楚，尽管我也认为您已经有了答案。 也许您在这里的主要问题是您正在使用相对较大的数据集，并且您无法查看每个阶段发生的情况。

因此，为了简单起见，我提出了一个较小的问题实例。请注意，在 stackoverflow 中询问时，这种方法可能更有助于获得帮助。事实上，您的文件可能太大而无法在某些计算机上下载和使用，从而避免人们参与其中。

但是让我们看一下代码。

library(tidyverse)
library(raster)
library(sf)
library(tmap)

您基本上有两个具有相同分辨率、位置和范围的栅格，如下所示：

R1 <- raster(ncol=10, nrow=10, xmn = 0, xmx=10, ymn = 0, ymx = 10)
R2 <- raster(ncol = 10, nrow = 10, xmn = 0, xmx=10, ymn = 0, ymx = 10)
values(R1) <- c(rep(1,ncell(R1)/2), rep(2,ncell(R1)/2))
values(R2) <- rep(10,ncell(R2))

使用这些文件，您可以直接执行许多操作。相同的范围和分辨率是一个加号：

Rsum <- R1 + R2
plot(R1) #See the legend of the plot
plot(R2) #See the legend of the plot
plot(Rsum) #See the legend of the plot

我不明白为什么您需要执行许多其他操作才能进行区域操作。

如果您的栅格不同，那将是有意义的。例如：

R2_Alt <- raster(ncol = 5, nrow = 2, xmn = 0, xmx=10, ymn = 0, ymx = 10) #Alt for alternative
values(R2_Alt) <- rep(10,ncell(R2))

在使用 R2_Al 的情况下，一个简单的操作可能会产生错误：

Test <- R2_Alt + R1

所以转移到矢量文件是有意义的：

R1_sp <- as(R1, "SpatialPolygonsDataFrame")
R2_Alt_sp <- as(R2_Alt, "SpatialPolygonsDataFrame")

R1_sf <- st_as_sf(R1_sp)
R2_Alt_sf <- st_as_sf(R2_Alt_sp)

names(R1_sf)[1] <- "V1" #Just to get differente variable names
names(R2_Alt_sf)[1] <- "V2_A"


tm_shape(R1_sf) + tm_polygons(border.col = "black") +
 tm_shape(R2_Alt_sf) + tm_polygons(border.col = "red", lwd = 4, alpha = 0)

这里有一个较小的问题实例可以阐明您的问题：

Int <- st_intersection(R2_Alt_sf, R1_sf)
View(Int)
plot(st_geometry(Int))

您还可以使用以前的 R2 栅格查看结果：

R2_sp <- as(R2, "SpatialPolygonsDataFrame")
R2_sf <- st_as_sf(R2_sp)
names(R2_Alt_sf)[1] <- "V2_A"

Int2 <- st_intersection(R2_sf,R1_sf)

最后，确保您需要的操作是交集。看 This.

我知道这个答案并不是真正的答案，但我希望它可以帮助您更接近。 让我知道您的想法并将所有反馈都放在这里。

一切顺利