R语言sf与osmdata工具：街道按所属行政区着色的实现问题

I'm trying to map the streets of Zaprešić and color them by their administrative district, using this code:

library(osmdata)
library(sf)
library(lwgeom)
library(tidyverse)
library(magrittr)
theme_set(theme_classic())
bb <- getbb("Zaprešić", featuretype = "city", format_out = "sf_polygon") %>% st_buffer(.05)
town <- getbb("Zaprešić", featuretype = "city") %>% opq() %>% add_osm_feature(key = "boundary", value = "administrative") %>% osmdata_sf() %>% trim_osmdata(bb)
streets <- getbb("Zaprešić", featuretype = "city") %>% opq() %>% add_osm_feature(key = "highway", value = c("residential", "primary", "secondary", "tertiary", "unclassified")) %>% osmdata_sf() %>% trim_osmdata(bb)
boundary <- getbb("Zaprešić", featuretype = "city") %>% opq() %>% add_osm_feature(key = "admin_level", value = "9") %>% osmdata_sf() %>% trim_osmdata(bb)
st_split(streets$osm_lines, boundary$osm_lines %>% filter(admin_level == 9) %>% st_buffer(.05)) %>% st_join(town$osm_multipolygons %>% filter(admin_level == 9)) -> street_split
ggplot() + geom_sf(data = street_split, aes(color = name.y)) + geom_sf(data = boundary$osm_lines %>% filter(admin_level == 9) %>% st_buffer(.05))

The issue is that some blue street lines are still visible under the black administrative boundary lines. I thought st_split would cut the streets at the boundary intersections to fix this cross-border coloring error, but it didn't work as expected. How can I adjust this to get the correct street coloring by district?

Let's break down what's going wrong and fix it step by step:

1. The Core Problem with Your Original Approach #

• You're using buffered boundary lines to split streets, not the actual administrative polygons. Buffering lines can introduce inaccuracies (like over-splitting or missing splits), and lines don't form closed shapes, so the split won't cleanly separate streets into district-specific segments.
• The default st_join uses a loose intersection check, which can still associate street segments with multiple districts or leave cross-border segments partially colored.

2. Revised Code with Explanations #

Here's the adjusted code that will correctly split streets by administrative districts and color them properly:

library(osmdata)
library(sf)
library(tidyverse)
library(magrittr)

theme_set(theme_classic())

# 1. Get bounding box and trim area
bb <- getbb("Zaprešić", featuretype = "city", format_out = "sf_polygon") %>% 
  st_buffer(.05)

# 2. Fetch streets (trimmed to bounding box)
streets <- getbb("Zaprešić", featuretype = "city") %>% 
  opq() %>% 
  add_osm_feature(key = "highway", value = c("residential", "primary", "secondary", "tertiary", "unclassified")) %>% 
  osmdata_sf() %>% 
  trim_osmdata(bb) %>% 
  .$osm_lines

# 3. Fetch administrative district POLYGONS (not lines)
district_polys <- getbb("Zaprešić", featuretype = "city") %>% 
  opq() %>% 
  add_osm_feature(key = "admin_level", value = "9") %>% 
  osmdata_sf() %>% 
  .$osm_multipolygons %>% 
  filter(admin_level == "9") %>% 
  st_transform(crs = st_crs(streets)) # Ensure matching CRS to avoid errors

# 4. Split streets using district polygons (not buffered lines)
street_segments <- st_split(streets, district_polys) %>% 
  st_collection_extract("LINESTRING") # Convert split collections to individual line features

# 5. Associate each street segment with its district (strictly within the polygon)
street_with_district <- street_segments %>% 
  st_join(district_polys, join = st_within, left = FALSE) %>% 
  # Remove any segments that didn't match a district (rare, but clean up)
  drop_na(name)

# 6. Plot the results
ggplot() +
  # Color streets by district name
  geom_sf(data = street_with_district, aes(color = name), linewidth = 0.8) +
  # Add district boundaries (no fill, thick black line)
  geom_sf(data = district_polys, fill = NA, color = "black", linewidth = 1.2) +
  # Optional: Improve legend and labels
  labs(title = "Zaprešić Streets by Administrative District", color = "District") +
  theme(plot.title = element_text(hjust = 0.5))

Key Improvements: #

• Using Administrative Polygons: Instead of buffering boundary lines, we use the actual closed district polygons to split streets. This ensures clean, accurate splits along district edges.
• Strict Spatial Join: st_within ensures a street segment is fully contained within a district (unlike the default intersection check), eliminating cross-border mismatches.
• CRS Alignment: We explicitly set the same coordinate reference system for streets and districts to avoid spatial operation errors.
• Cleaning Extracted Segments: st_collection_extract turns the split output into individual line features, which are easier to color and manipulate.

Troubleshooting Tips: #

• If you still see minor gaps or overlaps, try adjusting the st_buffer value on the bounding box slightly (e.g., .03 instead of .05) to refine the trimmed area.
• Ensure the admin_level value is correct for Zaprešić's districts—double-check OSM to confirm if level 9 is indeed the right administrative boundary.

内容的提问来源于stack exchange，提问作者J. Doe