# How to create a "random" shapefile

I have the following shapefile in R:

#https://stackoverflow.com/questions/75845500/r-find-out-list-of-boundary-coordinates
library(sf)
library(leaflet)
library(dplyr)

# built in shapefile
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>%
st_transform(st_crs(4326)) %>%
st_cast('POLYGON')

nc <- nc[nc\$NAME %in% c("Cleveland", "Gaston", "Rutherford", "Burke"), ]

nc_boundary_coords <- nc %>%
st_union() %>%
st_boundary() %>%
st_coordinates() %>%
data.frame()

plot(nc_boundary_coords\$X, nc_boundary_coords\$Y, type = "l")

I am trying to divide this area into "n" random (non-overlapping) subsections - for example, here n = 8:

In the end, I am trying to create a shapefile corresponding to these subsections (i.e. the geographical coordinates).

I found some posts in which I learned how to generate random points within a shapefile (e.g. https://stackoverflow.com/questions/57819700/creating-random-points-in-a-polygon-shapefile, https://stackoverflow.com/questions/75669453/r-check-to-see-if-coordinates-fall-inside-outside-of-a-shapefile-polygon) - but nothing similar to what I am trying to do.

How can I do this?

Note: I understand that this is an abstract problem and there are likely many different ways that this problem can be solved (i.e. using different algorithms) - I am interested in learning about different ways to solve this!

• There's an infinite number of ways to partition a polygon. This is an NP Hard problem, not suitable for resolution in a Q&A model. Your Note argues that we can't help you, because the Question is off-charter. Commented Mar 26, 2023 at 16:00
• @Vince I'm not sure there's enough definition in the question to say this is "NP Hard". The Q is intentionally vague though. You could tighten up the definition by defining what you mean by "random" and specifying any constraints on shape, area, perimeter, connectedness etc. There are methods for this eg using clustered voronoi polygons. Commented Mar 26, 2023 at 16:54
• Hey everyone, OP here - I know this is not fully random and the polygons are overlapping, but I created a simpler version of this problem here: stackoverflow.com/questions/75849216/r-checking-matrix-sums Commented Mar 26, 2023 at 16:57
• I don't see what that Q has to do with this one. Its more like your other recent question, about spatial joins... Commented Mar 26, 2023 at 18:17
• Can you generate the random points and then generate Voronoi polygons? If your points are random, then the output polygons would be more or less random. Commented Apr 14, 2023 at 18:17

Well, although the question is vague with several possible answers, let me share one potential approach. This consists on sampling 8 randoms points over the polygon and then compute Voronoi polygons based on that points. See:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#>     filter, lag
#> The following objects are masked from 'package:base':
#>
#>     intersect, setdiff, setequal, union

# built in shapefile
nc <- st_read(system.file("gpkg/nc.gpkg", package = "sf"), quiet = TRUE) %>%
st_transform(st_crs(4326)) %>%
st_cast("POLYGON")
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all
#> sub-geometries for which they may not be constant

nc <- nc[nc\$NAME %in% c("Cleveland", "Gaston", "Rutherford", "Burke"), ]

nc_boundary_pol <- nc %>%
st_union()

library(ggplot2)

sfplot <- ggplot(nc_boundary_pol) +
geom_sf()

sfplot

# Sample 8 points

set.seed(1234)
pp <- st_sample(nc_boundary_pol, size = 8) %>%
st_union()

# MULTIPOINT
pp
#> Geometry set for 1 feature
#> Geometry type: MULTIPOINT
#> Dimension:     XY
#> Bounding box:  xmin: -81.89617 ymin: 35.18283 xmax: -81.11418 ymax: 35.53374
#> Geodetic CRS:  WGS 84
#> MULTIPOINT ((-81.11418 35.18283), (-81.14613 35...

sfplot +
geom_sf(data = pp)

# Create Voronoi polygons
vors <- st_voronoi(pp, envelope = nc_boundary_pol) %>%
st_collection_extract("POLYGON") %>%
st_intersection(nc_boundary_pol)
#> Warning in st_voronoi.sfc(pp, envelope = nc_boundary_pol): st_voronoi does not
#> correctly triangulate longitude/latitude data

vors
#> Geometry set for 8 features
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -82.27856 ymin: 35.14917 xmax: -80.92607 ymax: 35.99478
#> Geodetic CRS:  WGS 84
#> First 5 geometries:
#> POLYGON ((-81.87042 35.18324, -81.76518 35.1826...
#> POLYGON ((-81.22789 35.15973, -81.22667 35.2162...
#> POLYGON ((-81.90607 35.86827, -81.67272 35.3516...
#> POLYGON ((-81.26856 35.40798, -81.28047 35.3080...
#> POLYGON ((-81.48643 35.47727, -81.67272 35.3516...
sfplot +
geom_sf(data = vors, fill = "blue", color = "white", alpha = 0.3) +
geom_sf(data = pp, color = "red")

Created on 2023-03-27 with reprex v2.0.2

• Just saw Spacedman comment, note this is in line with one of his suggestions. Commented Mar 27, 2023 at 7:08
• @diegherman: thank you so much for your answer! I am trying to use this question as a testing case for another question here - can you please take a look at it if you have time later? gis.stackexchange.com/questions/456233/… thank you so much! Commented Mar 27, 2023 at 7:15