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