1

I'm trying to create a hexagonal grid from a biome boundaries (cerrado, but it could be able to apply to any other political boundaries). In sf documentation I saw that it is possible to do it with the st_make_grid function, but does not clear how the size cell are defined. I would like to make a hexagonal grid with 25km2 size cell.

I've tried to use this function:

grid = sf::st_make_grid(x= shpCerrado, cellsize = 0.4, square = FALSE)

I used the size 0.4 because a friend compared in QGis and told that this size had a 25m side. But we dont have sure about the total area of the cell (which we would like that was 25km2).

I tried to find some package docs or another tutorial that explained how to define the cell area of a hexagonal grid, but I couldn't. Does anyone know how this calculation could be done? Or is there another way to make this grid in R?

Karoline
  • 17
  • 2

1 Answers1

1

You can derive the target cell size from the generic formula for area of hexagon; with a little algebra you will arrive at sqrt(2*cell_area/sqrt(3))

Or consider this piece of code; it covers county Ashe from the well known nc.shp that ships with {sf} in a 25 km² grid.

library(sf)
library(dplyr)

nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  slice(1) %>% # coiunty Ashe only
  st_geometry() %>% # no need for data anymore
  st_transform(3857) # transform to a metric CRS 

cell_area <- units::as_units(25, "km^2") # target grid size

grid_spacing <- sqrt(2*cell_area/sqrt(3)) # size of hexagon calculated from area

grid <- nc %>% 
  st_make_grid(square = F, cellsize = grid_spacing)  # make the grid
  
plot(nc, col = "red")
plot(grid, add = T)

enter image description here

Jindra Lacko
  • 7,814
  • 3
  • 22
  • 44
  • When I run to make the grid (the st_make_grid) I get this error: `Error in Ops.units(e1, int * e2) : both operands of the expression should be "units" objects` Also, I would like to ask if that a way to rasterize this polygon when it makes? – Karoline May 06 '23 at 18:36
  • strange... I can't reproduce that error. You may be able to sidestep it by specifying the cell area as 25000000 (i.e. 25 * 1000 ^ 2) and make sure you are in a metric CRS (3857 is not required, but a lat long such as 4326 would be a problem). As for the rasterizing - that is really a separate issue, `terra::rasterize()` would be my first consideration – Jindra Lacko May 06 '23 at 20:27
  • It works!!! When run with the example (25 * 1000 ^ 2). Thanks! I will try to rasterize this polygon now... – Karoline May 07 '23 at 16:51
  • I am glad it does! – Jindra Lacko May 07 '23 at 18:30