point in polygon with `sf`

July 26, 2022

i am trying to get that each point falls into the respective polygon. Given that I could not do it by aplying sp::over I am trying right now with the tidyverse::sf . I hope that someone can help me with this issue.

enter image description here


library(data.table)
library(sf)
library(sfheaders)
library(tidyverse)
library(mapview)
library(leaflet)
library(leafem)
library(tidyverse)

# Point data

coordinates = st_as_sf(data.frame(matrix(round(runif(n = 100, min = -10, max = 10),0), , 2), id = 1:(100)), coords = c("X1", "X2"))
mapview::mapview(coordinates)


# Polygon data


DT <- fread("ID       NW.X     NW.Y       NE.X     NE.Y       SE.X     SE.Y       SW.X       SW.Y value
  1 -9.5 9.5  -0.5 9.5 -0.5 0.5 -9.5 0.5 7
  2 -9.5 -0.5 -0.5 -0.5 -0.5 -9.5 -9.5 -9.5 14
  3 0.5 9.5 9.5 9.5 9.5 0.5 0.5 0.5 10
  4 0.5 -0.5 9.5 -0.5 9.5 -9.5 0.5 -9.5 8")

cols <- grep("^(NW|NE|SE|SW)\\.[XY]$", names(DT), value = TRUE)

DT[, (cols) := lapply(.SD, function(x) as.numeric(gsub(",", "\\.", x))), .SDcols = cols]
#set to workable format df
buffers <- setDF(DT) %>%
  # Melt to long, beep XY paired
  pivot_longer( cols = cols,
                names_to = c("point", ".value"),
                names_pattern = "(..)\\.(.)" ) %>%
  sfheaders::sf_polygon( x = "X", y = "Y", polygon_id = "ID" )

#visual incpection
mapview::mapview(buffers)


## Both spatial types

mapview::mapview(buffers) %>%
  leafem::addStaticLabels(
    label = buffers$ID,
    noHide = TRUE,
    direction = 'top',
    textOnly = TRUE,
    textsize = "20px")

mapview::mapview(coordinates) %>%
  leafem::addStaticLabels(
    label = coordinates$id,
    noHide = TRUE,
    direction = 'top',
    textOnly = TRUE,
    offset = c(0, 0),
    textsize = "12px")


mapview::mapview(buffers) +  
  mapview::mapview(coordinates)

I want that each points goes from their respective figure (point-in-polygon)



ggplot() +
  geom_sf(data=coordinates) +
  geom_sf(data=buffers) +
  theme_minimal()


points_sf_joined <- st_join(coordinates, buffers) %>%
  filter(!is.na(coordinates$id))

ggplot() +
  geom_sf(data=coordinates) +
  geom_sf(data=points_sf_joined) +
  theme_minimal()

With kind regards

>Solution :

The question as posted is somewhat unclear – what exactly is the problem?

The way I read it is that you have the point-in-polygon problem already solved (via the sf::st_join() call). The id of each polygon is stored in the ID column of your points_sf_joinded object (the point id is in id column, which may be kind of confusing btw.)

See the result of your last plot, slightly amended by me (I have removed the plotting of original coordinates not assigned to a polygon, and color coded the polygon IDs).

I have also added a geom_sf() call for the polygons to be plotted.

ggplot() +
  geom_sf(data = buffers) + # plot the polygons first
  geom_sf(data=points_sf_joined, aes(color = as.factor(ID))) +
  scale_color_brewer("polygon ID", palette = "Dark2") + 
  theme_minimal()