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

Question

我具有以下经过重新分类和批准的栅格，我正尝试对其进行拦截/叠加/重叠以获得新的栅格。这个想法是从栅格R1和R2的截取区域获得一个新的叠加栅格。完成此操作后，我将进行分区操作。HereR1，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_A1的情况下，简单的操作可能会产生错误：

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。

我知道这个答案并不是真正的答案，但我希望它能帮助您更紧密地联系。让我知道您的想法，并将所有反馈都放在这里。

一切顺利