Thank you Zivan! I guess I will take more time to digest the meaning of st_intersection(x). In particular, it contains points, polygons, multipolygons, and multilinestrings, it is such diverse an output that added complexity.
Xiang ________________________________ From: Zivan Karaman <zivan.kara...@gmail.com> Sent: Friday, August 8, 2025 22:42 To: Xiang Ye <yexiangm...@outlook.com> Cc: r-sig-geo@r-project.org <r-sig-geo@r-project.org> Subject: Re: [R-sig-Geo] How to retrieve the interior boundary (points) of a polygon sf object Hi Xiang, You are absolutely right; the first statement suffices here. As specified by the 'st_intersection' help: When called with missing y, the sfc method for st_intersection returns all non-empty intersections of the geometries of x; an attribute idx contains a list-column with the indexes of contributing geometries. So, you get all the intersections among the geometries within the sf(s) object (here, US states), and then you simply extract the lines to obtain the boundaries. Best, Zivan On Fri, Aug 8, 2025 at 4:27 PM Xiang Ye <yexiangm...@outlook.com<mailto:yexiangm...@outlook.com>> wrote: Dear Zivan, Thank you for your email. I followed your scripts, only to find it a genius solution. Actually, only the first sentence is enough: st_intersection(polygon) |> st_collection_extract('LINESTRING') And it has brought me the interior boundary already. However, I do not fully grasp the essence here. I know st_intersection(x, y) is to find the shared part of x and y, but what is st_intersection(x)? Conceptually it is "the shared part of x and x", but technically it brings me more features than x itself (and I feel confused about the output). But obviously you make the most of this feature. Thank you once again! Xiang ________________________________ From: R-sig-Geo <r-sig-geo-boun...@r-project.org<mailto:r-sig-geo-boun...@r-project.org>> on behalf of Zivan Karaman <zivan.kara...@gmail.com<mailto:zivan.kara...@gmail.com>> Sent: Thursday, August 7, 2025 03:07 To: Josiah Parry <josiah.pa...@gmail.com<mailto:josiah.pa...@gmail.com>> Cc: r-sig-geo@r-project.org<mailto:r-sig-geo@r-project.org> <r-sig-geo@r-project.org<mailto:r-sig-geo@r-project.org>> Subject: Re: [R-sig-Geo] How to retrieve the interior boundary (points) of a polygon sf object Hi, I'm not sure I understand what you want to achieve or why, but here's a try: # Load required libraries library(sf) library(spData) # Options to avoid messages from tigris polygon <- st_geometry(us_states) # Extract the shared borders only (where geometries touch more than once) internal_lines <- st_intersection(polygon) |> st_collection_extract("LINESTRING") # Now remove the external boundary # Step 1: Union all states into a single polygon (whole USA outline) national_outline <- st_union(polygon) |> st_boundary() # Step 2: Remove outer boundary from internal_lines internal_only <- st_difference(internal_lines, national_outline) # Plot result plot(st_geometry(national_outline), col = "blue", lwd = 0.5, main = "Interior Boundaries Only") plot(st_geometry(internal_only), col = "red", lwd = 0.5, add = TRUE) Best, Zivan On Wed, Aug 6, 2025 at 8:24 PM Josiah Parry <josiah.pa...@gmail.com<mailto:josiah.pa...@gmail.com>> wrote: > > I'm not sure I fully understand the objective. But this is how I would do it > with rsgeo—I'm not sure how to do it with sf > > library(sf) > library(spData) > library(rsgeo) > > union_geoms(as_rsgeo(us_states)) |> > expand_geoms() |> > .subset2(1) |> > head(1) |> > cast_geoms("multipoint") |> > st_as_sfc() |> > plot() > > > > > On Wed, Aug 6, 2025 at 10:56 AM Andrea Gilardi - Unimib > <andrea.gila...@unimib.it<mailto:andrea.gila...@unimib.it>> wrote: >> >> Hi Xiang! I would use the following approach: >> >> library(tidyverse) >> library(sf) >> library(spData) >> # You need the development version of sfnetworks: >> # remotes::install_github("luukvdmeer/sfnetworks") >> library(sfnetworks) # Exports st_duplicated >> >> polygon = st_geometry(us_states) >> pts_matrix <- st_coordinates(polygon) >> pts_sfc <- st_multipoint(pts_matrix[, c(1, 2)]) |> st_sfc(crs = >> st_crs(polygon)) |> st_cast("POINT") >> pts_dup <- st_duplicated(pts_sfc) >> pts_interior <- pts_sfc[pts_dup] >> >> plot(polygon, reset = FALSE) >> points(pts_interior, cex = 1, pch = 16, col = 2) >> >> You can also check the output here: >> https://gist.github.com/agila5/2b5d9b7e1453c4d3f292e7aae5fdb79e >> >> Hope that's helpful. >> >> Andrea >> >> On 8/6/2025 7:27 PM, Dexter Locke wrote: >> > library(tidyverse) >> > library(sf) >> > library(spData) >> > >> > polygon=st_geometry(us_states) >> > plot(polygon) >> > >> > polygon |> st_as_sf() |> summarise() |> plot() >> > >> > >> > On Wed, Aug 6, 2025 at 1:25 PM Xiang Ye >> > <yexiangm...@outlook.com<mailto:yexiangm...@outlook.com>> wrote: >> > >> >> Dear community, >> >> >> >> I am trapped by a seemingly easy task. I would like to write a short R >> >> script to retrieve the interior boundary (points) of a polygon sf object, >> >> like the state boundary of the Contiguous United States, but without the >> >> outline. >> >> >> >> I tried several methods, with none fulfilling my objective. Below is the >> >> code, with my observations in the comments: >> >> >> >> # Objective: Retrieve the internal boundary of a polygon data set >> >> >> >> library(tidyverse) >> >> library(sf) >> >> library(spData) >> >> >> >> polygon=st_geometry(us_states) >> >> plot(polygon) >> >> >> >> # Method 1: Use st_difference() >> >> st_union(polygon) -> border >> >> plot(border) >> >> st_difference(polygon, border) -> interior >> >> plot(interior) # You can tell that the result is weird. >> >> >> >> # Method 2: Retrieve vertices as points, then use st_difference() >> >> st_cast(polygon, 'POINT') -> points_polygon >> >> points_polygon # 3585 features >> >> st_cast(border, 'POINT') -> points_border >> >> points_border # 1275 features >> >> st_difference(points_polygon, points_border) -> points_interior # About 2 >> >> minutes to execute >> >> points_interior # 4.57 million features, about 2 GB. why? >> >> >> >> # Methods 3: Use st_disjoint instead of st_difference() >> >> polygon[border, op=st_disjoint] # Empty set >> >> points_polygon[border, op=st_disjoint] -> points_interior2 >> >> points_interior2 # 132 features; plot() it and doesn't look right >> >> points_polygon[points_border, op=st_disjoint] -> points_interior3 >> >> points_interior3 # 3585 features, meaning no points have been taken out >> >> >> >> # Methods 4: Retrieve vertices as coordinates, converting them into >> >> tibbles, then anti-join the border coordinates >> >> st_coordinates(polygon) %>% as_tibble -> coords_polygon >> >> coords_polygon# A tibble: 3,585 × 5 >> >> st_coordinates(border) %>% as_tibble -> coords_border >> >> coords_border # A tibble: 1,275 × 5 >> >> anti_join(coords_polygon, coords_border) -> coords_interior >> >> coords_interior # A tibble: 3,585 × 5, meaning no points have been taken >> >> out, and implying no coordinates in coords_border coincides with the >> >> coordinates in coords_polygon. This is just not reasonable. >> >> >> >> Any comments on why I fail and/or how to do it correctly will be helpful. >> >> Thank you in advance! >> >> >> >> 叶翔 YE, Xiang >> >> THINKING SPATIALLY<http://www.linkedin.com/in/spatialyexiang>. >> >> Ph.D. in Spatial Statistics >> >> >> >> >> >> [[alternative HTML version deleted]] >> >> >> >> _______________________________________________ >> >> R-sig-Geo mailing list >> >> R-sig-Geo@r-project.org<mailto:R-sig-Geo@r-project.org> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> >> > [[alternative HTML version deleted]] >> > >> > _______________________________________________ >> > R-sig-Geo mailing list >> > R-sig-Geo@r-project.org<mailto:R-sig-Geo@r-project.org> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org<mailto:R-sig-Geo@r-project.org> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org<mailto:R-sig-Geo@r-project.org> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org<mailto:R-sig-Geo@r-project.org> https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
