Message-ID: <4c416e6a-6ea9-4b7b-97c6-907e5a9ddf11@unimib.it>
Date: 2025-08-06T17:56:07Z
From: Andrea Gilardi
Subject: How to retrieve the interior boundary (points) of a polygon sf object
In-Reply-To: <CAA=SVwHtRijb2ZFwAP3XxpASrNethh9z=0MsXxaFZPVHLTd_BQ@mail.gmail.com>
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