使用来自其他两个已批准栅格的重叠批准栅格进行分区操作
Zonal operation using an overlapped-ratified raster which comes from two other ratified rasters
我有以下重新分类和批准的栅格,我正在尝试 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")
非常欢迎任何建议。
我觉得答案一点都不清楚,虽然我也觉得你已经有了答案。
也许您的主要问题是您正在处理相对较大的数据集,并且您无法看到每个阶段的情况。
因此,为了简单起见,我提出了一个较小的问题实例。请注意,在 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.
我知道这个答案不是真正的答案,但我希望它能帮助你更接近。
让我知道您的想法,并将所有反馈放在这里。
祝一切顺利
我有以下重新分类和批准的栅格,我正在尝试 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")
非常欢迎任何建议。
我觉得答案一点都不清楚,虽然我也觉得你已经有了答案。 也许您的主要问题是您正在处理相对较大的数据集,并且您无法看到每个阶段的情况。
因此,为了简单起见,我提出了一个较小的问题实例。请注意,在 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.
我知道这个答案不是真正的答案,但我希望它能帮助你更接近。 让我知道您的想法,并将所有反馈放在这里。
祝一切顺利