How to retrieve the interior boundary (points) of a polygon sf object
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 <andrea.gilardi at unimib.it>
Sent: Thursday, August 7, 2025 01:56
To: Xiang Ye <yexiangmaup at outlook.com>
Cc: r-sig-geo at r-project.org <r-sig-geo at r-project.org>
Subject: Re: [R-sig-Geo] How to retrieve the interior boundary (points) of a polygon sf object
Sent: Thursday, August 7, 2025 01:56
To: Xiang Ye <yexiangmaup at outlook.com>
Cc: r-sig-geo at r-project.org <r-sig-geo at r-project.org>
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 <yexiangmaup at 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 at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo