使用R交集使用两个shapefile图层创建多边形内部多边形键 [英] Using R intersections to create a polygons-inside-a-polygon key using two shapefile layers

问题描述

数据

我有两个shapefile标记了 national 和<巴基斯坦的a href ="https://data.humdata.org/dataset/provincial-regional-constituency-boundaries-pakistan" rel ="nofollow noreferrer">省选区.

目标

我正在尝试使用R创建一个密钥，该密钥将基于此数据中的坐标生成哪些省级选区包含在其中"或与哪些国家级选区相交的列表.例如，NA-01对应于PA-01，PA-02，PA-03；NA-02对应于PA-04和PA-05等.(该密钥最终将用于链接包含国家和省级选举结果的单独数据框；我已经确定了这一部分.)

我只有基本/中级的R技能，这些技能主要是通过反复试验而获得的，没有使用R之外的GIS数据的经验.

尝试的解决方案

针对此问题我能找到的最接近的解决方案来自本指南来计算R中的交点面积.但是，我无法成功复制三种建议的方法中的任何一种(提问者对交点使用常规的TRUE/FALSE报告，或者更精确地计算重叠面积).

代码

 #导入地图文件NA_map<-readOGR(dsn ="./National_Constituency_Boundary"，layer ="National_Constituency_Boundary")PA_map<-readOGR(dsn ="./Provincial_Constituency_Boundary"，layer ="Provincial_Constituency_Boundary")#这两个现在都是分别具有273和577个元素的SpatialPolygonsDataFrame对象.#如果相关，我使用spdpylr调整了一些数据属性名称(稍后在加入选举数据帧时使用):NA_map<-NA_map％>％重命名(constituency_number = NA_Cons，district_name =地区，省=省)PA_map<-PA_map％>％重命名(省=省，district_name = DISTRICT，constituency_number = PA)#计算交叉点，取一个结果<-gIntersects(NA_map，PA_map，byid = TRUE)#这会创建一个包含157,521个元素的大型矩阵行名(结果)<-NA_map @ data $ constituency_number列名(结果)<-PA_map @ data $ constituency_number 

但是，尝试添加行名/名称标签会给我错误消息:

  dimnames(x)中的错误<-dn:'dimnames'[1]的长度不等于数组范围 

没有行名/公称标签，我无法读取覆盖矩阵，也不确定如何过滤它们，以便仅生成有助于创建NA-PA密钥的TRUE交集列表.

我还尝试复制其他两个建议的解决方案，以计算出精确的重叠区域:

 #计算交叉点，取两个pi<-相交(NA_map，PA_map)#这会生成一个包含273个元素的SpatialPolygons对象面积<-data.frame(area = sapply(pi @ polygons，FUN = function(x){slot(x，'area')})))#计算相交面积，但没有其他变量row.names(areas)<-sapply(pi @ polygons，FUN = function(x){slot(x，'ID')}) 

这会生成错误消息:

 `row.names<-.. data.frame`(`* tmp *`，value = c("2"，"1"，"4"，"5"，:不允许重复的"row.names"另外:警告消息:设置"row.names"时的非唯一值:"1" 

因此，当我尝试将区域附加到属性信息时，

  attArrea<-spCbind(pi，区域) 

我收到错误消息

  spCbind(pi，areas)中的错误:行名不相同 

尝试第三种建议的方法:

 #计算交叉点，取三个pi<-st_intersection(NA_map，PA_map) 

产生错误消息:

UseMethod("st_intersection")中的

 错误:没有适用于"st_intersection"的适用方法应用于类"c('SpatialPolygonsDataFrame'，'SpatialPolygons'，'Spatial'，'SpatialPolygonsNULL'，'SpatialVector')"的对象 

我知道我的SPDF映射不能用于第三种方法，但是从描述中不清楚，需要什么步骤来转换它并尝试这种方法.

寻求帮助

对于使用这些方法中的任何一种所必需的校正建议，或指向解决该问题的其他方法的指针，将不胜感激.谢谢！

解决方案

以下是一些示例数据

 库(光栅)p<-shapefile(system.file("external/lux.shp"，package ="raster"))p1<-汇总(p，by ="NAME_1")p2<-p [，'NAME_2'] 

所以我们有带区域的p1和带较低层划分的p2.

现在我们可以做

  x<-相交(p1，p2)#或x<-union(p1，p2)data.frame(x) 

应该(和)与原著相同

  data.frame(p)[，c('NAME_1'，'NAME_2')] 

要获取多边形的面积，可以执行

  x $ area<-area(x)/1000000#除以得到km2 

由于边界的细微变化，可能会出现许多条"，非常小的多边形.那可能对你来说没关系.

但是另一种方法可能是通过质心进行匹配:

  y<-p2e<-提取(p1，坐标(p2))y $ NAME_1<-e $ NAME_1data.frame(y) 

The data

I have two shapefiles marking the boundaries of national and provincial electoral constituencies in Pakistan.

The objective

I am attempting to use R to create a key that will generate a list of which provincial-level constituencies are "contained within" or otherwise intersecting with which national-level constituencies, based on their coordinates in this data. For example, NA-01 corresponds with PA-01, PA-02, PA-03; NA-02 corresponds with PA-04 and PA-05, etc. (The key will ultimately be used to link separate dataframes containing electoral results at the national and provincial level; that part I've figured out.)

I have only basic/intermediate R skills learned largely through trial and error and no experience working with GIS data outside of R.

The attempted solution

The closest solution I could find for this problem comes from this guide to calculating intersection areas in R. However, I have been unable to successfully replicate any of the three proposed approaches (either the questioner's use of a general TRUE/FALSE report on intersections, or the more precise calculations of area of overlap).

The code

# 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

Attempting to add the rowname/colname labels, however, gives me the error message:

Error in dimnames(x) <- dn : 
  length of 'dimnames' [1] not equal to array extent

Without the rowname/colname labels, I'm unable to read the overlay matrix, and unsure how to filter them so as to produce a list of only TRUE intersections that would help make a NA-PA key.

I also attempted to replicate the other two proposed solutions for calculating exact area of overlap:

# 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')})

This generates the error message:

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’ 

So that when I attempt to attach areas to attributes info with

attArrea <- spCbind(pi, areas)

I get the error message

Error in spCbind(pi, areas) : row names not identical

Attempting the third proposed method:

# calculate intersections, take three
pi <- st_intersection(NA_map, PA_map)

Produces the error message:

Error in UseMethod("st_intersection") : 
  no applicable method for 'st_intersection' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector')"

I understand that my SPDF maps can't be used for this third approach, but wasn't clear from the description what steps would be needed to transform it and attempt this method.

The plea for help

Any suggestions on corrections necessary to use any of these approaches, or pointers towards some other method of figuring this, would be greatly appreciated. Thanks!

解决方案

Here is some example data

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p1 <- aggregate(p, by="NAME_1")
p2 <- p[, 'NAME_2']

So we have p1 with regions, and p2 with lower level divisions.

Now we can do

x <- intersect(p1, p2)
# or  x <- union(p1, p2)
data.frame(x)

Which should be (and is) the same as the original

data.frame(p)[, c('NAME_1', 'NAME_2')]

To get the area of the polygons, you can do

 x$area <- area(x) / 1000000  # divide to get km2

There are likely to be many "slivers", very small polygons because of slight variations in borders. That might not matter to you.

But another approach could be matching by centroid:

y <- p2
e <- extract(p1, coordinates(p2))
y$NAME_1 <- e$NAME_1
data.frame(y)

这篇关于使用R交集使用两个shapefile图层创建多边形内部多边形键的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！