polygonpoint

Why is st_join not returning a polygon for this point?


I'm trying to assign each snow station in my feature class from ArcGIS Online a district using sf::st_join.

All of the points in Dryden aren't getting a district and I cannot figure it out... they are clearly inside the district!

Here's a reproducible example - notice that SNOW-WLF-DY doesn't get a district. I tried setting sf_use_s2 to false, and also tried st_intersection, st_contains, etc.

enter image description here

library(tidyverse)
library(arcgisbinding)
library(sf)

# Check whether you are connected to your license
arc.check_product()

# Check whether you are connected to your login. 
arc.check_portal()

# Bring in and format data ----
sf_use_s2(TRUE)

# Get new snow stations layer and transform it to Lambert Conformal Conic 
snowlocs <- arc.select(arc.open("https://ws.lioservices.lrc.gov.on.ca/arcgis2/rest/services/LIO_OPEN_DATA/LIO_Open08/MapServer/31")) %>% 
  arc.data2sf() %>% 
  st_transform(3162) %>% 
  filter(STATION_NUMBER %in% c("SNOW-WLF-MK", "SNOW-WLF-DY", "SNOW-WLF-PK")) %>%
  group_by(STATION_NUMBER) %>% 
  slice(1)


# Get old districts layer and transform it to Lambert Conformal Conic 
districts_old <- arc.select(arc.open("https://services1.arcgis.com/TJH5KDher0W13Kgo/arcgis/rest/services/MNR_Administrative_Boundaries/FeatureServer/1")) %>% 
  arc.data2sf() %>% 
  st_transform(3162) %>% 
  filter(DISTRICT_NAME %in% c("Kenora", "Dryden", "Sioux Lookout"))

ggplot(districts_old) +
  geom_sf() +
  geom_sf(data = snowlocs) +
  geom_sf_text(data = snowlocs, aes(label = STATION_NUMBER), nudge_y = -10000)

# Why doesn't SNOW-WLF-DY get a district?
check <- st_join(snowlocs %>% select(STATION_NUMBER), districts_old %>% select(DISTRICT_NAME))
check

Solution

  • I posted this as a bug on github here: https://github.com/r-spatial/sf/issues/2408

    and found that my "districts" feature class was invalid.

    The operation runs correctly if I fix it with st_make_valid(districts_old).