# Divide spatial hexagon in equal areas (pies) in r

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.

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.

``````library(sf)
library(dplyr)

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") %>%
st_cast("POLYGON")

# 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:
st_area(g_int[1,])/st_area(g[1,])
#0.3333333 [1]

``````

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