如何通過sf找到點屬於哪個多邊形

Question

我有一個sf對象，其中包含通過.shp文件獲取的城市區域的多邊形信息（選區）。對於給定的經緯度對，我想確定它屬於哪個分區。我在想我可以利用sf::st_contains()，但我無法以正確的格式獲取緯度/經度。

Answer 1

在lon/lat上使用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) 
     ) 

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

Answer 3

這可以是「向量化的」。這裏有一個例子：

library(sf) 
library(tidyverse) 

新加坡shape文件：

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 

這可能比行方式更快OP但我會讓別人有樂趣的標杆：

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