Thank you, Andrea! ________________________________ From: Andrea Gilardi - Unimib <[email protected]> Sent: Friday, August 8, 2025 22:20 To: Xiang Ye <[email protected]> Cc: [email protected] <[email protected]> Subject: Re: [R-sig-Geo] How to retrieve the interior boundary (points) of a polygon sf object
Dear Xiang, thanks for your kind words. 1. The underlying assumption behind st_duplicated is that the polygons used to construct the multipolygon share some vertices in the interior part of the multipolygon. Your explanation and understandings are correct. If you want to check the implementation, see here: https://github.com/luukvdmeer/sfnetworks/blob/c6eb10a05768f570e94981e826de8ee3f7fb7fce/R/utils.R#L21 2. Yes, you are correct, you can simply use st_cast(..., "POINT"). I'm not sure whether one of the two approaches is much faster than the other, but I guess it doesn't matter here. Kind regards Andrea On 8/8/2025 8:39 AM, Xiang Ye wrote: Dear Andrea, Thank you for your detailed response. I followed your instructions carefully, and was able to get the expected results. St_duplicated() is cool! I have two following questions. 1. Why st_duplicated() works? My guess is that for a typical (multi)polygon sf object, all vertices on the interior borders will be adopted twice to form polygons. Technically, this makes them appear twice in the geometry column of the data set. Meanwhile, vertices on the exterior borders only appear once. When checked by st_duplicated(), all the second appearance of the interior-border vertices will be marked true, and thus they are extracted. Is my understanding correct? 1. Can we simplify a sentence? In the script, there is a sentence: pts_sfc <- st_multipoint(pts_matrix[, c(1, 2)]) |> st_sfc(crs = st_crs(polygon)) |> st_cast("POINT") Is this sentence equal to the following one? pts_sfc_2 <- st_cast(polygon, 'point') I experimented by myself, and find identical(pts_sfc, pts_sfc_2)==FALSE. However, identical(st_duplicated(pts_sfc), st_duplicated(pts_sfc_2))==TRUE. Again, thank you so much for your awesome solution! Xiang ________________________________ From: Andrea Gilardi - Unimib <[email protected]><mailto:[email protected]> Sent: Thursday, August 7, 2025 01:56 To: Xiang Ye <[email protected]><mailto:[email protected]> Cc: [email protected]<mailto:[email protected]> <[email protected]><mailto:[email protected]> Subject: Re: [R-sig-Geo] How to retrieve the interior boundary (points) of a polygon sf object 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 > <[email protected]><mailto:[email protected]> 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 >> [email protected]<mailto:[email protected]> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > [email protected]<mailto:[email protected]> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
