Search code examples
rgeometrygislatitude-longitudeshapefile

How can I map a lat/long to the nearest point of my polygon/shapefile in R

I am trying to map a set of lats/longs to my shape file to obtain CSDUID/census divisions information.

However, there are some outliers or lat/longs that are outside of the shapefile. I was wondering if it would be possible to map those points to the nearest point of the shapefile so that I can obtain division information for the later part of my analysis?

Reproducible example

Let's say I have a shapefile like gr

nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264)                
gr = st_sf(
    label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
    geom = st_make_grid(nc))

# gr
# Simple feature collection with 100 features and 1 field
# geometry type:  POLYGON
# dimension:      XY
# bbox:           xmin: 406262.2 ymin: 48374.87 xmax: 3052887 ymax: 1044158
# epsg (SRID):    2264
# proj4string:    +proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
# First 10 features:
#    label                           geom
# 1   J  1 POLYGON ((406262.2 48374.87...
# 2   J  2 POLYGON ((670924.7 48374.87...
# 3   J  3 POLYGON ((935587.2 48374.87...
# 4   J  4 POLYGON ((1200250 48374.87,...
# 5   J  5 POLYGON ((1464912 48374.87,...
# 6   J  6 POLYGON ((1729575 48374.87,...
# 7   J  7 POLYGON ((1994237 48374.87,...
# 8   J  8 POLYGON ((2258900 48374.87,...
# 9   J  9 POLYGON ((2523562 48374.87,...
# 10  J 10 POLYGON ((2788225 48374.87,...

And I have a point that is outside the bounding box:

point = st_point(c(106160,28370))

library(ggplot2)
ggplot() + geom_sf(data = gr) + geom_sf(data = point)

enter image description here

How can I map points to the closest polygon ?

Solution

  • You have great functions for spatial analysis, especially spatial joins, in sf package. Assuming your point data are named point and your polygons are named.... polygons (both sf objects), you can use sf::st_join for spatial joins.

    st_join accepts different types of spatial mapping. I think you are interested in st_nearest_feature merge (first remove points that fall within your polygons, e.g. by a first st_join using st_within)

    library(sf)

st_join(
  points,
  polygons,
  join = st_nearest_feature
)

    In the edited example, that would give:

    point = st_sf(st_sfc(point), crs = st_crs(gr))
st_join(point, gr, join = st_nearest_feature)

    see ?st_join doc.