I have a grid of hexagons that is being used to create a model per hexagon. I would like for each hexagon to be divided into 3 or 6 equal areas so that I can use 2/3 of the hexagon total area to train a model, and a 1/3 to test the model.

Here is an example of hexagon:

hexa <- st_sf(ID = 1, crs = 4326,
              geometry = st_sfc(st_polygon(list(cbind(c(-96.2379631928876, -94.0319170874439, -92.0048940427835, -92.224603469399,-94.3568579552715, -96.3420584229642, -96.2379631928876), c(41.3710446366938, 42.1956644084384, 41.1433453337767, 39.3108662668349, 38.521263223804, 39.5306362872028, 41.3710446366938))))))
hexa_m <- st_transform(hexa, crs = 2163)

I have not been able to use the st_subdivide function for this purpose.

One thought is to cast the hexagon to points, add a center point, then create lines between adjacent points, to ultimately transform the multilines object into multiple polygons, but I am not sure how to create lines only with adjacent points. Here is the center point if that is an idea you feel could work.

ctr_pt <- st_sf(ID = 1, 
                crs = 4326,
                geometry = st_sfc(st_point(c(-94.20018, 40.35845))))
ctr_pt_m <- st_transform(ctr_pt, crs = 2163)
hexa_9pts <- st_union(hexa_pts, ctr_pt_m)

I'm leaning toward using pies of a hexagon, but in reality, I think I could use any 1/3 area block in each hexagon.

Here is a image of what I would like at the end, i.e., one hexagon, divided into 2/3 - 1/3 polygons. To add a little bit of complexity, some hexagons, at the edge of a species' distribution are not complete. Not sure what to do about those but I would still like to use them. enter image description here

1 Answer 1


A hexagon may be divided in three equal parts by intersecting it with a same sized hexagon that has been displaced the length of the side (consecutive vertices); conveniently this can be queried with st_distance(hex_points[1,], hex_point[2,]), where hex_points is the vertices of the polygon casted with st_cast to points.

After that we add this distance to every Y coordinate and then intersect the two polygon layers. This distance is the same as (1/sin(60))*apothem-length, where the apothem is a line segment from the center to the midpoint of one of its sides.


nc = read_sf(system.file("shape/nc.shp", package="sf")) %>% st_transform(32616)
g = st_make_grid(nc, square = F) 

# let's check the distance from the bottom corner of the first polygon to the first center
# in reality is (1/sin(60))*apotheme (obviously transformed to radians)
# in hexagons, this distance is the same as the side length
poly_p = st_cast(g[1,], "POINT")
apothem_distance = st_distance(poly_p[1,], poly_p[2,])

g_apothem = st_coordinates(g) %>% as.data.frame() %>%
  mutate(Y = Y + units::drop_units(apothem_distance[[1]])) %>%
  st_as_sf(coords = c("X", "Y"), crs = 32616) %>%
  group_by(L2) %>% summarise(do_union = F) %>% st_cast("LINESTRING") %>% 

# give each hexagon a number
g_row = g %>% st_as_sf() %>% mutate(pol_num = as.factor(row_number()))
# intersect
g_int = st_intersection(g_row, g_apothem) 
# for some reason a point is included, let's filter that out
plot(g_int[,1] %>% filter(st_geometry_type(.) == "POLYGON"))

# a 3rd of the area:
#0.3333333 [1]

enter image description here

  • Beautiful solution: this works perfectly and can also be easily adapted if I model one hexagon at the time through a function. Thank you!
    – JoeBird
    Commented Dec 24, 2021 at 16:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.