使用 R 交集使用两个 shapefile 层创建多边形内的多边形键
Using R intersections to create a polygons-inside-a-polygon key using two shapefile layers
数据
我有两个 shapefile 标记了巴基斯坦 national and provincial 个选区的边界。
objective
我正在尝试使用 R 创建一个密钥，该密钥将根据这些数据中的坐标生成哪些省级选区 "contained within" 或以其他方式与哪些国家级选区相交的列表。例如，NA-01对应PA-01、PA-02、PA-03； NA-02 对应于 PA-04 和 PA-05 等。（关键将最终用于 link 分离包含国家和省级选举结果的数据框；那部分我已经弄清楚了。）
我只有 basic/intermediate R 技能，主要是通过反复试验学到的，没有使用 R 之外的 GIS 数据的经验。
尝试的解决方案
对于这个问题，我能找到的最接近的解决方案来自 this guide 计算 R 中的交叉区域。但是，我无法成功复制所提出的三种方法中的任何一种（提问者使用一般 TRUE/FALSE 报告交点，或更精确的重叠面积计算）。
代码
# import map files
NA_map <- readOGR(dsn = "./National_Constituency_Boundary", layer = "National_Constituency_Boundary")
PA_map <- readOGR(dsn = "./Provincial_Constituency_Boundary", layer = "Provincial_Constituency_Boundary")
# Both are now SpatialPolygonsDataFrame objects of 273 and 577 elements, respectively.
# If relevant, I used spdpylr to tweak some of data attribute names (for use later when joining to electoral dataframes):
NA_map <- NA_map %>%
rename(constituency_number = NA_Cons,
district_name = District,
province = Province)
PA_map <- PA_map %>%
rename(province = PROVINCE,
district_name = DISTRICT,
constituency_number = PA)
# calculate intersections, take one
Results <- gIntersects(NA_map, PA_map, byid = TRUE)
# this creates a large matrix of 157,521 elements
rownames(Results) <- NA_map@data$constituency_number
colnames(Results) <- PA_map@data$constituency_number
尝试添加 rowname/colname 标签时，出现错误消息：
Error in dimnames(x) <- dn :
length of 'dimnames' [1] not equal to array extent
没有 rowname/colname 标签，我无法读取覆盖矩阵，并且不确定如何过滤它们以生成有助于制作 NA-PA 密钥的仅 TRUE 交叉点列表。
我还尝试复制其他两个建议的解决方案来计算准确的重叠区域：
# calculate intersections, take two
pi <- intersect(NA_map, PA_map)
# this generates a SpatialPolygons object with 273 elements
areas <- data.frame(area=sapply(pi@polygons, FUN = function(x) {slot(x, 'area')}))
# this calculates the area of intersection but has no other variables
row.names(areas) <- sapply(pi@polygons, FUN=function(x) {slot(x, 'ID')})
这会生成错误消息：
Error in `row.names<-.data.frame`(`*tmp*`, value = c("2", "1", "4", "5", :
duplicate 'row.names' are not allowed
In addition: Warning message:
non-unique value when setting 'row.names': ‘1’
因此，当我尝试使用
将区域附加到属性信息时
attArrea <- spCbind(pi, areas)
我收到错误消息
Error in spCbind(pi, areas) : row names not identical
正在尝试第三种方法：
# calculate intersections, take three
pi <- st_intersection(NA_map, PA_map)
产生错误信息：
Error in UseMethod("st_intersection") :
no applicable method for 'st_intersection' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector')"
我知道我的 SPDF 地图不能用于这第三种方法，但从描述中不清楚需要哪些步骤来转换它并尝试这种方法。
请求帮助
对于使用这些方法中的任何一种所必需的更正的任何建议，或指向其他计算方法的指示，将不胜感激。谢谢！
您的代码不是独立的，因此我没有尝试复制您报告的错误。
但是，使用 sf 包（旨在取代 rgeos、rgdal 和 sp 包，获得您想要的 'key' 非常简单在不远的将来）。看这里：
library(sf)
# Download shapefiles
national.url <- 'https://data.humdata.org/dataset/5d48a142-1f92-4a65-8ee5-5d22eb85f60f/resource/d85318cb-dcc0-4a59-a0c7-cf0b7123a5fd/download/national-constituency-boundary.zip'
provincial.url <- 'https://data.humdata.org/dataset/137532ad-f4a9-471e-8b5f-d1323df42991/resource/c84c93d7-7730-4b97-8382-4a783932d126/download/provincial-constituency-boundary.zip'
download.file(national.url, destfile = file.path(tempdir(), 'national.zip'))
download.file(provincial.url, destfile = file.path(tempdir(), 'provincial.zip'))
# Unzip shapefiles
unzip(file.path(tempdir(), 'national.zip'), exdir = file.path(tempdir(), 'national'))
unzip(file.path(tempdir(), 'provincial.zip'), exdir = file.path(tempdir(), 'provincial'))
# Read map files
NA_map <- st_read(dsn = file.path(tempdir(), 'national'), layer = "National_Constituency_Boundary")
PA_map <- st_read(dsn = file.path(tempdir(), 'provincial'), layer = "Provincial_Constituency_Boundary")
# Get sparse list representation of intersections
intrs.sgpb <- st_intersects(NA_map, PA_map)
length(intrs.sgpb) # One list element per national constituency
# [1] 273
print(intrs.sgpb[[1]]) # Indices of provnicial constituencies intersecting with first national constituency
# [1] 506 522 554 555 556
print(PA_map$PROVINCE[intrs.sgpb[[1]]])[1] # Name of first province intersecting with first national constituency
# [1] KHYBER PAKHTUNKHWA
这是一些示例数据
library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p1 <- aggregate(p, by="NAME_1")
p2 <- p[, 'NAME_2']
所以我们有 p1 区域和 p2 较低级别的分区。
现在我们可以做
x <- intersect(p1, p2)
# or x <- union(p1, p2)
data.frame(x)
应该（并且）与原来的一样
data.frame(p)[, c('NAME_1', 'NAME_2')]
要获得多边形的面积，您可以
x$area <- area(x) / 1000000 # divide to get km2
可能有很多 "slivers"，非常小的多边形，因为边界略有不同。这对你来说可能无关紧要。
但另一种方法可能是通过质心匹配：
y <- p2
e <- extract(p1, coordinates(p2))
y$NAME_1 <- e$NAME_1
data.frame(y)
数据
我有两个 shapefile 标记了巴基斯坦 national and provincial 个选区的边界。
objective
我正在尝试使用 R 创建一个密钥，该密钥将根据这些数据中的坐标生成哪些省级选区 "contained within" 或以其他方式与哪些国家级选区相交的列表。例如，NA-01对应PA-01、PA-02、PA-03； NA-02 对应于 PA-04 和 PA-05 等。（关键将最终用于 link 分离包含国家和省级选举结果的数据框；那部分我已经弄清楚了。）
我只有 basic/intermediate R 技能，主要是通过反复试验学到的，没有使用 R 之外的 GIS 数据的经验。
尝试的解决方案
对于这个问题，我能找到的最接近的解决方案来自 this guide 计算 R 中的交叉区域。但是，我无法成功复制所提出的三种方法中的任何一种（提问者使用一般 TRUE/FALSE 报告交点，或更精确的重叠面积计算）。
代码
# import map files
NA_map <- readOGR(dsn = "./National_Constituency_Boundary", layer = "National_Constituency_Boundary")
PA_map <- readOGR(dsn = "./Provincial_Constituency_Boundary", layer = "Provincial_Constituency_Boundary")
# Both are now SpatialPolygonsDataFrame objects of 273 and 577 elements, respectively.
# If relevant, I used spdpylr to tweak some of data attribute names (for use later when joining to electoral dataframes):
NA_map <- NA_map %>%
rename(constituency_number = NA_Cons,
district_name = District,
province = Province)
PA_map <- PA_map %>%
rename(province = PROVINCE,
district_name = DISTRICT,
constituency_number = PA)
# calculate intersections, take one
Results <- gIntersects(NA_map, PA_map, byid = TRUE)
# this creates a large matrix of 157,521 elements
rownames(Results) <- NA_map@data$constituency_number
colnames(Results) <- PA_map@data$constituency_number
尝试添加 rowname/colname 标签时，出现错误消息：
Error in dimnames(x) <- dn :
length of 'dimnames' [1] not equal to array extent
没有 rowname/colname 标签，我无法读取覆盖矩阵，并且不确定如何过滤它们以生成有助于制作 NA-PA 密钥的仅 TRUE 交叉点列表。
我还尝试复制其他两个建议的解决方案来计算准确的重叠区域：
# calculate intersections, take two
pi <- intersect(NA_map, PA_map)
# this generates a SpatialPolygons object with 273 elements
areas <- data.frame(area=sapply(pi@polygons, FUN = function(x) {slot(x, 'area')}))
# this calculates the area of intersection but has no other variables
row.names(areas) <- sapply(pi@polygons, FUN=function(x) {slot(x, 'ID')})
这会生成错误消息：
Error in `row.names<-.data.frame`(`*tmp*`, value = c("2", "1", "4", "5", :
duplicate 'row.names' are not allowed
In addition: Warning message:
non-unique value when setting 'row.names': ‘1’
因此，当我尝试使用
将区域附加到属性信息时attArrea <- spCbind(pi, areas)
我收到错误消息
Error in spCbind(pi, areas) : row names not identical
正在尝试第三种方法：
# calculate intersections, take three
pi <- st_intersection(NA_map, PA_map)
产生错误信息：
Error in UseMethod("st_intersection") :
no applicable method for 'st_intersection' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector')"
我知道我的 SPDF 地图不能用于这第三种方法，但从描述中不清楚需要哪些步骤来转换它并尝试这种方法。
请求帮助
对于使用这些方法中的任何一种所必需的更正的任何建议，或指向其他计算方法的指示，将不胜感激。谢谢！
您的代码不是独立的，因此我没有尝试复制您报告的错误。
但是，使用 sf 包（旨在取代 rgeos、rgdal 和 sp 包，获得您想要的 'key' 非常简单在不远的将来）。看这里：
library(sf)
# Download shapefiles
national.url <- 'https://data.humdata.org/dataset/5d48a142-1f92-4a65-8ee5-5d22eb85f60f/resource/d85318cb-dcc0-4a59-a0c7-cf0b7123a5fd/download/national-constituency-boundary.zip'
provincial.url <- 'https://data.humdata.org/dataset/137532ad-f4a9-471e-8b5f-d1323df42991/resource/c84c93d7-7730-4b97-8382-4a783932d126/download/provincial-constituency-boundary.zip'
download.file(national.url, destfile = file.path(tempdir(), 'national.zip'))
download.file(provincial.url, destfile = file.path(tempdir(), 'provincial.zip'))
# Unzip shapefiles
unzip(file.path(tempdir(), 'national.zip'), exdir = file.path(tempdir(), 'national'))
unzip(file.path(tempdir(), 'provincial.zip'), exdir = file.path(tempdir(), 'provincial'))
# Read map files
NA_map <- st_read(dsn = file.path(tempdir(), 'national'), layer = "National_Constituency_Boundary")
PA_map <- st_read(dsn = file.path(tempdir(), 'provincial'), layer = "Provincial_Constituency_Boundary")
# Get sparse list representation of intersections
intrs.sgpb <- st_intersects(NA_map, PA_map)
length(intrs.sgpb) # One list element per national constituency
# [1] 273
print(intrs.sgpb[[1]]) # Indices of provnicial constituencies intersecting with first national constituency
# [1] 506 522 554 555 556
print(PA_map$PROVINCE[intrs.sgpb[[1]]])[1] # Name of first province intersecting with first national constituency
# [1] KHYBER PAKHTUNKHWA
这是一些示例数据
library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p1 <- aggregate(p, by="NAME_1")
p2 <- p[, 'NAME_2']
所以我们有 p1 区域和 p2 较低级别的分区。
现在我们可以做
x <- intersect(p1, p2)
# or x <- union(p1, p2)
data.frame(x)
应该（并且）与原来的一样
data.frame(p)[, c('NAME_1', 'NAME_2')]
要获得多边形的面积，您可以
x$area <- area(x) / 1000000 # divide to get km2
可能有很多 "slivers"，非常小的多边形，因为边界略有不同。这对你来说可能无关紧要。
但另一种方法可能是通过质心匹配：
y <- p2
e <- extract(p1, coordinates(p2))
y$NAME_1 <- e$NAME_1
data.frame(y)