为一组多边形 R 在每个多边形中找到最大点答案

【问题标题】：Find maximum point with in each polygon for a set of polygons R为一组多边形 R 在每个多边形中找到最大点
【发布时间】：2014-11-22 19:32:48
【问题描述】：

我确信这个问题已经在其他地方得到了回答，但我无法通过搜索找到它。

我有代表一个国家内城市的点以及每个城市的人口。我还有一个县的多边形文件。我想找到每个县内最大城市的位置。

如何做到这一点？

这是一些数据

结构（列表（国家= c（“我们”，“我们”，“我们”，“我们”，“我们”，“我们”，“我们”，“我们”，“我们”，“我们” , “我们”,
“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们”、“我们” ", "我们"), 城市 = c("cabarrus", "cox store", "cal-vel", "briarwood townhouses", "barker heights", "davie
十字路口”、“蟹点村”、“杜鹃花”、“切斯特菲尔德”、“查尔斯蒙特”、“康纳”、“三叶草花园”、“科里赫高地”、“卡利森”、“克雷斯特维尤英亩”、“克莱格”、“迦南”公园”，“尚蒂伊”，“贝尔格莱德”，“布里斯十字路口”，“虚张声势”，“巴特纳”，“底部”，“班迪”，“博斯蒂安高地”），AccentCity = c（“卡巴鲁斯”，“考克斯商店” , “Cal-Vel”, “Briarwood Townhouses”, “Barker Heights”, “Davie Crossroads”, “Crab Point Village”, “Azalea”, “Chesterfield”, “Charlesmont”, “Connor”, “Clover Garden”, “ Corriher Heights”、“Callisons”、“Crestview Acres”、“Clegg”、“Canaan Park”、“Chantilly”、“Belgrade”、“Brices Crossroads”、“Bluff”、“Butner”、“Bottom”、“Bandy” , "波斯蒂安高地"), 区域 = c("NC", "NC", "NC", "NC", "NC", "NC", "NC", "NC", "NC", "NC" 、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“NC”、“ NC", "NC", "NC"), 人口 = c(NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, A _integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_，NA_integer_），纬度=（35.2369444，35.275，36.4291667，35.295，35.3111111，35.8319444，34.7602778，35.58，35.81，5.9341667， 35.7419444，36.1883333，35.5605556，35.0841667，35.0213889，35.8655556，36.2761111，36.3016667，34.88，34.8186111，35.8377778，36.1319444，36.4747222，35.6419444，35.7544444），经度= C（-80.5419444，-82.0352778，-78.9694444，-81.5238889，-82.4441667， -80.535, -76.7305556, -82.4713889, -81.6611111, -81.5127778, -78.1486111, -79.4630556, -80.635, -76.7255556, -80.5427778, -78.8497222, -79.7852778, -76.1711111, -77.2352778, -78.1016667, -82.8580556, -78.7569444, -80.7741667, -81.09, -80.9294444)), .Names = c("Country", "City", "Accentity" "人口", "纬度", "经度"), row.names = c(544L, 889L, 551L, 434L, 190L, 975L, 894L, 147L, 717L, 700L, 831L, 773L, 862L, 559L, 915L, 753L, 584L, 695L, 262L, 437L, 372L, 537L, 406L, 178L, 02L), class= "data.frame")

还有一些在北卡罗来纳州需要阅读的代码

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
                IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))

plot(xx)

我想找到每个县内人口最多的城市。对不起，我没有可重复的例子。如果我这样做了，我会得到答案！

【问题讨论】：

如果您提供 reproducible example 来显示您正在尝试完成的工作，这可能是最好的回答
发布dput(data) 的输出，其中data 是您的“一些数据”。
抱歉，感谢您提供的信息。我不知道如何在问题中输入数据。
是的，但你还没有这样做。我所做的只是将您发布的数据放入代码块中。这还不够。输入 dput(data) 并将结果复制/粘贴到问题中。
我不明白：在你提供的数据中，Population 都是NA??顺便说一句：你不应该向我发表评论（使用@jlhoward），否则我不会收到通知。

标签： r gis spatial maptools sp

【解决方案1】：

简短的回答是您应该在包rgeos 中使用gContains(...)。

这里是长答案。

在下面的代码中，我们从 GADM 数据库中获取北卡罗来纳州县的高分辨率 shapefile，并从美国地质调查局数据库中获取北卡罗来纳州城市的地理编码数据集。后者已经有县信息，但我们忽略了这一点。然后我们使用gContains(...) 将城市映射到相应的县，将该信息添加到城市数据框中，并使用 data.table 包确定每个县中最大的城市。大部分工作在接近尾声的 4 行代码中完成。

library(raster)   # for getData(...);   you may not need this
library(foreign)  # for read.dbf(...);  you may not need this
library(rgeos)    # for gContains(...); loads package sp as well

setwd("< directory for downloaded data >")
# get North Carolina Counties shapefile from GADM database
USA         <- getData("GADM",country="USA",level=2)   # level 2 is counties
NC.counties <- USA[USA$NAME_1=="North Carolina",]      # North Carolina Counties
# get North Carolina Cities data from USGS database
url <- "http://dds.cr.usgs.gov/pub/data/nationalatlas/citiesx010g_shp_nt00962.tar.gz"
download.file(url,"cities.tar.gz")
untar("cities.tar.gz")
data      <- read.dbf("citiesx010g.dbf",as.is=TRUE)
NC.data   <- data[data$STATE=="NC",c("NAME","COUNTY","LATITUDE","LONGITUDE","POP_2010")]
## --- evverything up to here is just to set up the example

# convert cities data.frame to SpatialPointsDataFrame
NC.cities <- SpatialPointsDataFrame(NC.data[,c("LONGITUDE","LATITUDE")],
                                    data=NC.data,
                                    proj4string=CRS(proj4string(NC.counties)))
# map cities to counties
city.cnty   <- gContains(NC.counties,NC.cities,byid=TRUE)
# add county information to cities data
NC.data$county <- apply(city.cnty,1,function(cnty)ifelse(any(cnty),NC.counties@data[cnty,]$NAME_2,NA))
# identify largest city in each county
library(data.table)
result <- setDT(NC.data)[,.SD[which.max(POP_2010)],by="county"]
head(result)
#      county             NAME   COUNTY LATITUDE LONGITUDE POP_2010
# 1:  Jackson        Cullowhee  Jackson 35.31371 -83.17653     6228
# 2:   Graham     Robbinsville   Graham 35.32287 -83.80740      620
# 3:   Wilkes North Wilkesboro   Wilkes 36.15847 -81.14758     4245
# 4:    Rowan        Salisbury    Rowan 35.67097 -80.47423    33662
# 5: Buncombe        Asheville Buncombe 35.60095 -82.55402    83393
# 6:    Wayne        Goldsboro    Wayne 35.38488 -77.99277    36437

这里的主力是这条线：

city.cnty   <- gContains(NC.counties,NC.cities,byid=TRUE)

这会将 SpatialPointsDataFrame NC.Cities 中的每个点与 SpatialPolygonsDataFrame NC.counties 中的每个多边形进行比较，并返回一个逻辑矩阵，其中行代表城市，列代表县，如果是城市，[i,j] 元素是 TRUE i 在县j，否则FALSE。我们在下一条语句中逐行处理矩阵：

NC.data$county <- apply(city.cnty,1,function(cnty)ifelse(any(cnty),NC.counties@data[cnty,]$NAME_2,NA))

连续使用每一行索引NC.counties的属性表以提取县名。

您在问题中提供的数据存在一些问题，但具有指导意义。首先，maptools 包中的 NC shapefile 分辨率相对较低。特别是这意味着一些沿海岛屿完全消失了，因此其中一个岛屿上的任何城市都不会映射到一个县。您的真实数据可能存在同样的问题，因此请小心。

其次，将原始 USGS 数据集中的 COUNTY 列与我们添加的 county 列进行比较，有 3 个（865 个）县不同意。事实证明，在这些情况下，USGS 数据库是错误的（或过时的）。你可能也有同样的问题，所以也要小心。

第三，另外三个城市没有映射到任何县。这些都是沿海城市，可能反映了北卡罗来纳州形状文件中的小错误。你晚上也有这个问题。

【讨论】：

是的，这给了我一个结果，可以在数据框中找到记录该县的最大城市。然而，我的特定应用程序假设县是未知的，只能从多边形数据集中确定。
所以，我将拥有一个点和种群的坐标文件以及一个多边形的 shapefile。我想要每个多边形中人口最多的点的位置。我尝试使用 sp 包中的 over(pts, pol) ，但它只给了我位于整个文件中的点。它没有告诉我哪些点在哪些县。有没有办法做到这一点？感谢您抽出宝贵时间。我很感激。（在我的例子中，我在非洲国家的地区有城市，但我用 NC 作为一个易于理解的例子）。
在某处上传您的真实数据（坐标文件和完整的 shapefile，其中包含多个文件），然后发布链接。
好的。因此，我完全修改了答案，以演示如何使用城市坐标和县 shapefile 将城市映射到县。