Hexagonal

This is an R markdown detailing how to construct hexagonal bins in R with several packages.

Here we load libraries:

library(raster)

## Loading required package: sp

library(maps)
library(mapdata)
library(tidyverse)

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.5.0 
## ✔ readr   2.1.2      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ purrr::map()     masks maps::map()
## ✖ dplyr::select()  masks raster::select()

library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

First we get a map of the USA via the map package. I also reprojected the map with the EPSG:5070 projection (NAD83/Conus Albers). It’s the most accurate I think for the North America.

set.seed(1) 

#get US map and convert to hex grids
USA <- sf::st_as_sf(maps::map("usa", plot = FALSE, fill = TRUE)) %>%
  filter(ID == "main")

#reproject the map
hex <- st_transform(USA, crs=st_crs("EPSG:5070"))
ggplot() + geom_sf(data =hex)

Below is a function I made using some previous attempts by othersd with the sf package. This time I use only the sf package for ease and speed

#leo's version of making the gris

make_grid2 <- function(x, cell_diameter, cell_area, clip = FALSE) {
  if (missing(cell_diameter)) {
    if (missing(cell_area)) {
      stop("Must provide cell_diameter or cell_area")
    } else {
      cell_diameter <- sqrt(2 * cell_area / sqrt(3))
    }
  }
  ext <- as(extent(x) + cell_diameter, "SpatialPolygons")
  projection(ext) <- projection(x)
  # generate array of hexagon centers
  g <- st_make_grid(ext, cellsize = cell_diameter, crs = st_crs(ext),
                    square = FALSE)
  
  # clip to boundary of study area
  if (clip) {
    g <- st_intersection(g, st_make_valid(x))
  } else {
    g <- g[x, ]
  }
  # clean up feature IDs
  
  return(g)
}

NOTE! The unit for cell_area is in meters! So 1,000km2 is 1000000000 m2!

hex_grid<-make_grid2(hex, cell_area = 1000000000, clip = TRUE)

plot(hex_grid)

#You can save your shapefile here.
#saveRDS(hex_grid, "Data/cells/Grids/US_100KM2.RDS")

Hexagonal_Bins

Leo_Ohyama

2022-12-19

Hexagonal