我的csv文件A中有表示事件空间坐标的数据,还有另一个文件shapefile B,其中包含
head(A)
month longitude latitude lsoa_code crime_type
1 2014-09 -1.550626 53.59740 E01007359 Anti-social behaviour
2 2014-09 -1.550626 53.59740 E01007359 Public order
3 2014-09 -1.865236 53.93678 E01010646 Anti-social behaviour
head(B@data)
code name altname
0 E05004934 Longfield, New Barn and Southfleet <NA>
1 E05000448 Lewisham Central <NA>
2 E05003149 Hawcoat <NA>
我想将犯罪数据A加入到我的shapefile B中,以映射发生在我区域A中的犯罪事件。无法执行基于属性连接的
code
,因为A中的代码引用的单元不同于B中的代码。我阅读了许多教程和文章,但找不到答案。我尝试了:
joined = over(A, B)
和
overlay
,但是没有实现我想要的功能。是否有办法直接进行此连接或将从A到其他格式的中间转换是否需要?
有人遇到这个问题并解决了吗?
#1 楼
spaceEco包中的point.in.poly函数返回与sp多边形对象相交的点的SpatialPointsDataFrame对象,并可以选择添加多边形属性。首先让我们添加require包并创建一些示例数据。
require(spatialEco)
require(sp)
data(meuse)
coordinates(meuse) = ~x+y
sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349)))),'1')
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373)))),'2')
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791)))),'4')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4, y=runif(4)))
现在,让我们快速查看数据并将其绘制出来。
head(srdf@data) # polygons
head(meuse@data) # points
plot(srdf)
points(meuse, pch=20)
最后,我们可以将点与多边形相交。结果将是一个SpatialPointsDataFrame对象,在这种情况下,该对象具有srdf多边形数据中包含的两个额外属性(PIDS,y)。
pts.poly <- point.in.poly(meuse, srdf)
head(pts.poly@data)
不是多边形数据中的唯一标识列,您可以轻松地添加一个。srdf@data$poly.ids <- 1:nrow(srdf)
一旦我们将点和多边形相交,我们就可以使用唯一多边形对点进行汇总在多边形数据中是属性的ID。
# Number of points in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=length)
# Mean lead in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=mean)
评论
@ arvi1000,是的,但是sp :: point.in.polygon产生逻辑。 spatialEco:point.in.poly是over的包装,但是返回一个sp SpatialPointsDataFrame和快捷方式,用于关联多边形属性的某些步骤,就像raster:intersect对rgeos :: gIntersect所做的一样。
–杰弗里·埃文斯(Jeffrey Evans)
16年4月14日在16:35
sp :: point.in.polygon实际上返回一个数值(0 =点在外部,1 =内部,2 =在边缘,3 =在顶点)。在某些情况下可能是正确的事情。认为在此处进行注释很有帮助,因为这是相关条款的Google最佳搜索结果
– arvi1000
16年4月14日在18:02
#2 楼
软件包over()
中的sp
可能有点令人困惑,但效果很好。我假设您已经使用coordinates(A) <- ~longitude+latitude
将“ A”空间化了:同样没有作为A行,并从B中的每个相交多边形中获得一个变量“代码”。# Overlay points and extract just the code column:
a.data <- over(A, B[,"code"])
评论
我发现over()在多边形的顶点处存在点问题,尽管我认为这是到目前为止我发现的最简单的解决方案。
– JMT2080AD
16年7月20日在19:51
你有什么问题?
–Simbamangu
16年7月21日在10:48
排除。我需要进一步探索。今天晚些时候我会为您提供一些数据,如果您感兴趣,我们可以一起查看。我可能是错的,但是我很确定算法中有一些退化需要处理,至少对于我的数据而言。
– JMT2080AD
16年7月21日在17:33
没关系。我的数据一定有问题。这个实验装置工作正常。 r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
– JMT2080AD
16年7月21日在22:59
这是比接受的答案简单得多的方法,并且不需要安装除rgdal之外的其他软件包。
–布莱斯·弗兰克(Bryce Frank)
18年2月9日在20:35
#3 楼
这是类似dplyr的解决方案: library(spdplyr)
ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
what = "sp")
pop <- read_excel("data/SAPE20DT7-mid-2017-parlicon-syoa-estimates-unformatted.xls",sheet = "data")
pop <- janitor::clean_names(pop)
ukcounties_pop <- ukcounties %>% inner_join(pop, by = c("pcon18nm" = "pcon11nm"))
人口数据来自:https:/ /www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates
我必须将下载的形状文件转换为geoJson:
https://geoportal.statistics。 gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1
您可以通过以下方法做到:
<预备班=“ lang-r prettyprint-override”>
uk_constituencies <- readOGR("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC.shp")
uk_constituencies # this is in tmerc format. we need to convert it to WGS84 required by geoJson format.
# First Convert to Longitude / Latitude with WGS84 Coordinate System
wgs84 = '+proj=longlat +datum=WGS84'
uk_constituencies_trans <- spTransform(uk_constituencies, CRS(wgs84))
# Convert from Spatial Dataframe to GeoJSON
uk_constituencies_json <- geojson_json(uk_constituencies_trans)
# Save as GeoJSON file on the file system.
geojson_write(uk_constituencies_json, file = "data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson")
#read back in:
ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
what = "sp")
评论
您是否查看过sp包中的point.in.polygon()?@ arvi1000我有,将再试一次。我对point.in.polygon的想法是,这是否会保留变量month和Crime_type。你知道吗?
我对point.in.poly进行了更多尝试,最后选择了属于相关多边形的那些点。谢谢。
然后,也许您应该用解决方案回答自己的问题。请记住,这个网站的目的就是好的答案。