Identifying lines that touch the edges of a polygon using sf package

I am doing spatial joins using the sf package in R. Specifically, I am attempting to find the roads (lines) that represent the borders of each census tract (polygons) in NJ. Note: for many census tracts in NJ, the borders of census tracts are all roads.

I have attempted to do a join between the lines (roads) and the polygons (census tracts) using st_touches, which I believe identifies geometries that share a common point but their interiors do not intersect. Using st_touches did successfully identify the roads the comprise the borders of the census tract, but it also identified roads inside the census tract that touch the edge of the census tract.

Using the image below as an example, when I performed the line-polygon join using st_touches, the roads that were identified as touching were the 4 roads in red, which represent the actual borders of the census tract. However, it also identified the 3 roads in blue, which are in the interior of the census tract polygon.

This could just be a misunderstanding of polygons, but I’m wondering whether a spatial join between a polygon and a line treats the polygon itself as a sort of “closed” line. In other words, I’m thinking that the join is treating the coordinates that comprise the edges of the polygon as lines, such that any of all of the actual lines (roads) that touch the edges of polygon are flagged as touching edges of the polygon, regardless of whether the roads are inside or on top of the polygon. When using st_touches, is it the case, that lines on the interior of polygon should not be identified as touching the polygon because they are on the inside of the polygon?

I’m sorry if this question is a little (or very) confusing, so if a reproducible example is needed, I can produce one. However, I just wanted to pose the question to see if someone could answer it without sample data. Thanks!

It looks suspicious that all three incorrect roads are inside the polygon, none are outside it. I would expect the opposite – you should get boundary roads, plus outside roads (orthogonal to the boundary) that share a single point with the boundary.

A code and data sample could help. But what I think could have happened:

you are using sf with spherical s2 mode (it is now default).
the track polygons are inverted here – i.e. they describe all the surface of the Earth except the census track.
having the inverted polygon, the result is technically correct, the “inner” roads are outside of the polygon, and touch it.

I’m assuming the question is still about

/../ find the roads (lines) that represent the borders of each census tract /../

As roads and tracts from tigris / cencus.gov seem to be aligned at vertex level, we can use st_covered_by predicate to join line sections of the tract polygon to streets:

split tract polygon to segments with nngeo::st_segments()
match each segment with the street it is covered by
optionally group segments by resulting tract and street name, combine

As cencus.gov practices some weird geo-blocking polices and I’m not able to use tigris from my local system, I’ve included relevant example datasets ( nj_roads is cropped by extended bbox of nj_tracts_172 )

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr)
library(mapview)

# check input data 
mapview(list(nj_tracts_172, nj_roads))

street_sections <- 
  nj_tracts_172[,"NAME"] |> 
  nngeo::st_segments() |> 
  st_set_geometry("geometry") |> 
  st_join(nj_roads[, "FULLNAME"], join = st_covered_by) |> 
  summarise(across(geometry, st_union), .by = c(NAME, FULLNAME))

street_sections
#> Simple feature collection with 4 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -74.03679 ymin: 40.76193 xmax: -74.03105 ymax: 40.76612
#> Geodetic CRS:  NAD83
#>   NAME     FULLNAME                       geometry
#> 1  172      West St LINESTRING (-74.03679 40.76...
#> 2  172      21st St LINESTRING (-74.03503 40.76...
#> 3  172 Palisade Ave LINESTRING (-74.03105 40.76...
#> 4  172      17th St LINESTRING (-74.03287 40.76...

mapview(street_sections, zcol = "FULLNAME")

Input data:

nj_roads <- 
  readr::read_csv(
'LINEARID,FULLNAME,RTTYP,MTFCC,WKT
1105646772963,Palisade Ave,M,S1400,"LINESTRING (-74.0330693190423 40.761629043074,-74.032868 40.761929,-74.032597 40.762332,-74.032576 40.762362,-74.03246 40.762536,-74.032427 40.762584,-74.03216 40.762981,-74.032035 40.763167,-74.03168 40.763695,-74.031614 40.763797,-74.031046 40.764678,-74.030752 40.765135,-74.030647 40.765308,-74.030646 40.7653096875007)"
110460669543,18th St,M,S1400,"LINESTRING (-74.036229 40.764303,-74.036932 40.764553,-74.037146 40.764629,-74.037188 40.7646439211607)"
110460669544,18th St,M,S1400,"LINESTRING (-74.03246 40.762536,-74.033806 40.763052,-74.035079 40.763539,-74.036399 40.764045)"
110460669554,19th St,M,S1400,"LINESTRING (-74.032035 40.763167,-74.033408 40.763684,-74.034677 40.764161,-74.035996 40.764657)"
110460669545,18th St,M,S1400,"LINESTRING (-74.0309355232216 40.7616290078124,-74.031666 40.761925,-74.032597 40.762332)"
110460669555,19th St,M,S1400,"LINESTRING (-74.030646 40.7623152151955,-74.031168 40.762563,-74.03216 40.762981)"
110460669308,Co Rd 687,C,S1400,"LINESTRING (-74.031891080755 40.7616290284621,-74.031666 40.761925,-74.031168 40.762563,-74.03216 40.762981)"
110460669569,21st St,M,S1400,"LINESTRING (-74.031046 40.764678,-74.032462 40.765188,-74.033732 40.765646,-74.035035 40.766116,-74.0356563323247 40.7664160331199)"
110460671762,Monastery Pl,M,S1400,"LINESTRING (-74.031614 40.763797,-74.033014 40.764311,-74.034285 40.764777,-74.035599 40.76526)"
110460670092,West St,M,S1400,"LINESTRING (-74.037188 40.762835019253,-74.036788 40.763452,-74.036399 40.764045,-74.036229 40.764303,-74.035996 40.764657,-74.035599 40.76526,-74.035035 40.766116,-74.0348412764384 40.7664160424894)"
1105646773006,Bergenline Ave,M,S1400,"LINESTRING (-74.036345182014 40.7616290207297,-74.036328 40.761655,-74.035916 40.762276,-74.035479 40.762935,-74.035079 40.763539,-74.034909 40.763796,-74.034677 40.764161,-74.034285 40.764777,-74.033732 40.765646,-74.03329 40.766341,-74.0332428880122 40.7664160442151)"
110460669538,17th St,M,S1400,"LINESTRING (-74.032868 40.761929,-74.034195 40.76244,-74.035479 40.762935,-74.036788 40.763452,-74.037188 40.7636105099737)"
1105646774597,New York Ave,M,S1400,"LINESTRING (-74.0347119936964 40.7616290434476,-74.034631 40.761756,-74.034195 40.76244,-74.033806 40.763052,-74.033408 40.763684,-74.033014 40.764311,-74.032462 40.765188,-74.032219 40.765575,-74.03204 40.76585,-74.031669368442 40.7664160243737)"', 
show_col_types = FALSE) |> 
  st_as_sf(wkt = "WKT", crs = "NAD83") |> 
  st_set_geometry("geometry")

nj_tracts_172 <- 
  readr::read_csv(
'STATEFP,COUNTYFP,TRACTCE,GEOID,NAME,NAMELSAD,MTFCC,FUNCSTAT,ALAND,AWATER,INTPTLAT,INTPTLON,WKT
34,017,017200,34017017200,172,Census Tract 172,G5020,S,125563,0,40.7640375,-74.033921,"MULTIPOLYGON (((-74.036399 40.764045,-74.036229 40.764303,-74.035996 40.764657,-74.035599 40.76526,-74.035035 40.766116,-74.033732 40.765646,-74.032462 40.765188,-74.031046 40.764678,-74.031614 40.763797,-74.03168 40.763695,-74.032035 40.763167,-74.03216 40.762981,-74.032427 40.762584,-74.03246 40.762536,-74.032576 40.762362,-74.032597 40.762332,-74.032868 40.761929,-74.034195 40.76244,-74.035479 40.762935,-74.036788 40.763452,-74.036399 40.764045)))"',
show_col_types = FALSE) |> 
  st_as_sf(wkt = "WKT", crs = "NAD83") |> 
  st_set_geometry("geometry")

nj_roads
#> Simple feature collection with 13 features and 4 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -74.03719 ymin: 40.76163 xmax: -74.03065 ymax: 40.76642
#> Geodetic CRS:  NAD83
#> # A tibble: 13 × 5
#>         LINEARID FULLNAME       RTTYP MTFCC                             geometry
#>  *         <dbl> <chr>          <chr> <chr>                     <LINESTRING [°]>
#>  1 1105646772963 Palisade Ave   M     S1400 (-74.03307 40.76163, -74.03287 40.7…
#>  2  110460669543 18th St        M     S1400 (-74.03623 40.7643, -74.03693 40.76…
#>  3  110460669544 18th St        M     S1400 (-74.03246 40.76254, -74.03381 40.7…
#>  4  110460669554 19th St        M     S1400 (-74.03203 40.76317, -74.03341 40.7…
#>  5  110460669545 18th St        M     S1400 (-74.03094 40.76163, -74.03167 40.7…
#>  6  110460669555 19th St        M     S1400 (-74.03065 40.76232, -74.03117 40.7…
#>  7  110460669308 Co Rd 687      C     S1400 (-74.03189 40.76163, -74.03167 40.7…
#>  8  110460669569 21st St        M     S1400 (-74.03105 40.76468, -74.03246 40.7…
#>  9  110460671762 Monastery Pl   M     S1400 (-74.03161 40.7638, -74.03301 40.76…
#> 10  110460670092 West St        M     S1400 (-74.03719 40.76284, -74.03679 40.7…
#> 11 1105646773006 Bergenline Ave M     S1400 (-74.03635 40.76163, -74.03633 40.7…
#> 12  110460669538 17th St        M     S1400 (-74.03287 40.76193, -74.03419 40.7…
#> 13 1105646774597 New York Ave   M     S1400 (-74.03471 40.76163, -74.03463 40.7…

nj_tracts_172
#> Simple feature collection with 1 feature and 12 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -74.03679 ymin: 40.76193 xmax: -74.03105 ymax: 40.76612
#> Geodetic CRS:  NAD83
#> # A tibble: 1 × 13
#>   STATEFP COUNTYFP TRACTCE     GEOID  NAME NAMELSAD MTFCC FUNCSTAT  ALAND AWATER
#> *   <dbl> <chr>    <chr>       <dbl> <dbl> <chr>    <chr> <chr>     <dbl>  <dbl>
#> 1      34 017      017200    3.40e10   172 Census … G5020 S        125563      0
#> # ℹ 3 more variables: INTPTLAT <dbl>, INTPTLON <dbl>,
#> #   geometry <MULTIPOLYGON [°]>

Filed under: Kiến thức lập trình - @ 04:59

Thẻ: rgeospatial

Thiết kế website giá rẻ

Danh mục

Identifying lines that touch the edges of a polygon using sf package