如何通过sf找到一个点属于哪个多边形

Question

我有一个sf对象，包含一个城市地区的多边形信息（选民区），是通过一个.shp文件获得的。对于给定的纬度/经度对，我想确定它属于哪个选民区。我考虑可以利用sf::st_contains()，但是在获取正确格式的纬度/经度方面遇到了问题。

- kevinykuo

---

我发现使用 sp::point.in.polygon 很顺利（虽然只是用 sp，而不是 sf）。 - r2evans

如果您提供一些示例数据，将更容易帮助您。 - SymbolixAU

同时，在两个 sf 对象上使用 sf::st_join()。您可以指定 join 函数为 st_within 以获取多边形中的点，并且它也会返回一个 sf 对象。 - SymbolixAU

Answer 1

这可以被“向量化”。下面是一个例子：

library(sf)
library(tidyverse)

新加坡矢量文件：

singapore <- st_read("~/data/master-plan-2014-subzone-boundary-no-sea-shp/MP14_SUBZONE_NO_SEA_PL.shp", quiet=TRUE, stringsAsFactors=FALSE)
singapore <- st_transform(singapore, 4326)

回收中心的CSV文件：

centers <- read_csv("~/data/recycl.csv")
glimpse(centers)
## Observations: 407
## Variables: 10
## $ lng             <dbl> 104.0055, 103.7677, 103.7456, 103.7361, 103.8106, 103.962...
## $ lat             <dbl> 1.316764, 1.296245, 1.319204, 1.380412, 1.286512, 1.33355...
## $ inc_crc         <chr> "F8907D68D7EB64A1", "ED1F74DC805CEC8B", "F48D575631DCFECB...
## $ name            <chr> "RENEW (Recycling Nation's Electronic Waste)", "RENEW (Re...
## $ block_house_num <chr> "10", "84", "698", "3", "2", "1", "1", "1", "357", "50", ...
## $ bldg_name       <chr> "Changi Water Reclamation Plant", "Clementi Woods", "Comm...
## $ floor           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ post_code       <int> 498785, 126811, 608784, 689814, 159047, 486036, 39393, 55...
## $ street          <chr> "Changi East Close", "West Coast Road , Clementi Woods Co...
## $ unit            <chr> "(Lobby)", "#B1-01 (Management Office)", "(School foyer)"...

将^^转换为一个简单的特性对象：

map2(centers$lng, centers$lat, ~st_point(c(.x, .y))) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf(centers[,-(1:2)], .) -> centers_sf

这可能比逐行操作更快，但我会让其他人来进行有趣的基准测试:

bind_cols(
  centers,
  singapore[as.numeric(st_within(centers_sf, singapore)),]
) %>% 
  select(lng, lat, inc_crc, subzone_name=SUBZONE_N) %>% 
  mutate(subzone_name = str_to_title(subzone_name))
## # A tibble: 407 x 4
##         lng      lat          inc_crc               subzone_name
##       <dbl>    <dbl>            <chr>                      <chr>
##  1 104.0055 1.316764 F8907D68D7EB64A1             Changi Airport
##  2 103.7677 1.296245 ED1F74DC805CEC8B             Clementi Woods
##  3 103.7456 1.319204 F48D575631DCFECB              Teban Gardens
##  4 103.7361 1.380412 1F910E0086FD4798                 Peng Siang
##  5 103.8106 1.286512 55A0B9E7CBD34AFE             Alexandra Hill
##  6 103.9624 1.333555 C664D09D9CD5325F                      Xilin
##  7 103.8542 1.292778 411F79EAAECFE609                  City Hall
##  8 103.8712 1.375876 F4516742CFD4228E Serangoon North Ind Estate
##  9 103.8175 1.293319 B05B32DF52D922E7            Alexandra North
## 10 103.9199 1.335878 58E9EAF06206C772            Bedok Reservoir
## # ... with 397 more rows

- hrbrmstr

Answer 2

回复晚了，我自己在寻找答案。

最终得到了这个结果：

library(sf)
library(tidyverse)

nc = st_read(system.file("shape/nc.shp", package="sf"),
             stringsAsFactors = FALSE)
d <-
  data_frame(lon = runif(1e3, -84.5, -75.5),
             lat = runif(1e3,  34,  36.6),
             somevariable = rnorm(1e3, 1000, 100))

geo_inside <- function(lon, lat, map, variable) {

  variable <- enquo(variable)
  # slow if lots of lons and lats or big sf - needs improvement
  pt <-
    tibble::data_frame(x = lon,
                       y = lat) %>%
    st_as_sf(coords = c("x", "y"), crs = st_crs(map))
  pt %>% st_join(map) %>% pull(!!variable)

}

d <-
  d %>%
  mutate(county = geo_inside(lon, lat, nc, NAME))

glimpse(d)
Observations: 1,000
Variables: 4
$ lon          <dbl> -79.68728, -79.06104, -83.92082, -76.36866, -75.8635...
$ lat          <dbl> 36.11349, 35.67239, 35.08802, 35.78083, 36.55786, 34...
$ somevariable <dbl> 910.9803, 1010.6816, 919.3937, 924.0845, 1154.0975, ...
$ county       <chr> "Guilford", "Chatham", "Cherokee", "Tyrrell", NA, NA...

d %>%
  ggplot() +
  geom_sf(data = nc) +
  geom_point(aes(lon, lat, colour = county)) +
  theme(legend.position = "none")

虽然速度不是很快，但似乎能够胜任工作。

Einar

- einar

Answer 3

使用 st_point() 函数处理经纬度数据，就可以与其他sf函数配合使用。

示例：

find_precinct <- function(precincts, point) {
  precincts %>%
    filter(st_contains(geometry, point) == 1) %>%
    `[[`("WARDS_PREC")
}


ggmap::geocode("nicollet mall, st paul") %>%
  rowwise() %>%
  mutate(point = c(lon, lat) %>%
           st_point() %>%
           list(),
         precinct = find_precinct(msvc_precincts, point)
         )

- kevinykuo

---

有没有一种向量化的方法来实现这个，即不使用“逐行”方式？ - RoyalTS

Answer 4

如果您有一个坐标的数据框（mydf），请将其转换为一个sf对象，然后与一个多边形sf map相交：

mydf_sf <- sf::st_as_sf(mydf, coords=c("lon","lat"), crs=4326)
int <- sf::st_intersects(mydf_sf , map)
mydf$country <- map$country_name[unlist(int)]

这里有一个完整的工作示例，位于https://gis.stackexchange.com/a/318629/36710

- Berry Boessenkool