在每個多邊形中查找一組多邊形的最大點R

Question

我確定這個問題已經在其他地方得到了解答，但是我一直沒能通過搜索來想出它。

我有代表一個國家的城市的點數以及每個城市的人口。我也有縣的多邊形文件。我想找到每個縣內最大城市的位置。

這怎麼辦？

下面是一些數據

結構（名單（國家= C（ 「我們」， 「我們」， 「我們」， 「我們」， 「我們」， 「我們」， 「我們」，「我們我們「，」我們「，」我們「，」我們「，」我們「，」我們「，」我們「，」我們「 「，」我們「，」我們「，」我們「，」我們「，」我們「），城市= c（」cabarrus「，」cox商店「，」cal- vel「，」briarwood排屋「，」barker高度「，」davie
十字路口「，」螃蟹點村「，」杜鵑花「，」切斯特菲爾德「，」查爾斯蒙特「，」康納「，」三葉草花園「，」corriher高地「，」callisons「 ，「clegg」，「canaan park」，「chantilly」，「belgrade」，「brices crossroads」，「bluff」，「butner」，「bottom」，「bandy」，「bostian heights」），AccentCity = c（「 Cabarrus「，」C牛扒店「，」Cal-Vel「，」Briarwood Townhouses「，」Barker Heights「，」Davie Crossroads「，」Crab Point Village「，」Azalea「，」Chesterfield「，」Charlesmont「，」Connor「 「，」Corriher Heights「，」Callisons「，」Crestview Acres「，」Clegg「，」Canaan Park「，」Chantilly「，」Belgrade「，」Brices Crossroads「，」Bluff「，」Butner「，」Bottom「 「Bandy」，「Bostian Heights」），Region = c（「NC」，「NC」，「NC」，「NC」，「NC」，「NC」，「NC」，「NC」 NC，NC，NC，NC，NC，NC，NC，NC，NC，NC，NC 「，」NC「，」NC「，」NC「），人口= c（NA_整體，NA整體，NA整體，NA整體，NA整體NA_整體NA_整體NA_整體NA整體NA_整體NA_整體NA整體NA整體NA_整體NA整體，NA_整體__，NA_整體__，NA_整體__，NA_整體__，整數_整數_整數_整數_整數_整數_整數_整數_整數_整數_整數n_整數_整數_整數_ 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， - ），.Names = c（「Country」，「City」，「AccentCity」，「Region」，「Population」，「Latitude」，「Longitude」），row.names （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」）

and some code to read in nort h carolina

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

plot(xx) 

我想找到每個縣內人口最多的城市。對不起，我沒有一個可重複的例子。如果我這樣做，我會得到答案！

Answer 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.counties提取縣名

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

：我們處理矩陣行方向的下一條語句。

您在問題中提供的數據有一些問題，但這些問題仍具有指導意義。首先，maptools包中的NC shapefile分辨率相對較低。特別是這意味着一些沿海島嶼完全缺失，因此其中一個島嶼上的任何城市都不會映射到一個縣。你可能與你的真實數據有同樣的問題，所以要留意它。其次，將原始USGS數據集中的COUNTY列與我們添加的county列進行比較，其中有3個（不超過865個）不同意的縣。事實證明，在這些情況下，USGS數據庫是錯誤的（或過時的）。你可能有同樣的問題，所以要小心。

三，另外三個城市沒有映射到任何縣。這些都是沿海城市，可能反映了北卡羅萊納州造型文件中的不準確之處。你晚上也有這個問題。