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

Question

我有以下重新分类和批准的栅格，我正在尝试 intercept/overlay/overlap 获得一个新栅格。这个想法是从栅格 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

这里是 intercept/overlay/overlap

的代码

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 区域的 combination/interception 形成的叠加区域（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

我觉得答案一点都不清楚，虽然我也觉得你已经有了答案。也许您的主要问题是您正在处理相对较大的数据集，并且您无法看到每个阶段的情况。

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

但让我们转到代码。

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.

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

祝一切顺利

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

Zonal operation using an overlapped-ratified raster which comes from two other ratified rasters

raster

r-raster

sp

sf