14

I have a sf object that contains polygon information (precincts) for a metro area, obtained through a .shp file. For a given lat/lon pair, I want to determine which precinct it belongs to. I'm thinking I can utilize sf::st_contains() but am having trouble getting the lat/lon in the right format.

3
  • I've found good luck using sp::point.in.polygon (though just with sp, not with sf).
    – r2evans
    Apr 17, 2017 at 19:14
  • If you provide some example data it will be easier to help you
    – SymbolixAU
    May 3, 2017 at 23:05
  • 4
    also, use sf::st_join() on two sf objects. You can specify the join function to be st_within to get the points in polygons, and it will return you an sf object too.
    – SymbolixAU
    Feb 6, 2018 at 0:24

4 Answers 4

8

This can be "vectorized". Here's an example:

library(sf)
library(tidyverse)

Singapore shapefile:

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 of recycling centers:

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)"...

Turn ^^ into a simple features object:

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

This is likely faster than the row-wise op but I'll let someone else have fun benchmarking:

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
2

Use st_point() on the lon/lat then it can work with other sf functions.

Example:

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)
         )
1
  • 1
    Any way to do this in a vectorized way, i.e. not rowwise?
    – RoyalTS
    Apr 28, 2017 at 20:11
2

Late response, was searching for an answer myself.

Ended up with this:

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")

Not happy with the speed though, but seems to do the job.

Einar

1

If you have a data.frame of coordinates (mydf), convert them to an sf object and then intersect with an sf map of polygons:

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)]

There's a complete working example at https://gis.stackexchange.com/a/318629/36710

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.