Create hexagon grid for Finland

hexfin.R

library(magick)
library(raster)
library(rgeos)
library(sf)
library(sp)
library(tidyverse)

set.seed(1)

# Helper functions --------------------------------------------------------

# From http://strimas.com/spatial/hexagonal-grids/

make_grid <- 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(raster::extent(x) + cell_diameter, "SpatialPolygons")
  raster::projection(ext) <- raster::projection(x)
  # generate array of hexagon centers
  g <- sp::spsample(ext, type = "hexagonal", cellsize = cell_diameter,
                    offset = c(0.5, 0.5))
  # convert center points to hexagons
  g <- sp::HexPoints2SpatialPolygons(g, dx = cell_diameter)
  # clip to boundary of study area
  if (clip) {
    g <- rgeos::gIntersection(g, x, byid = TRUE)
  } else {
    g <- g[x, ]
  }
  # clean up feature IDs
  row.names(g) <- as.character(1:length(g))
  return(g)
}

dst_file <- "data/kuntajako_2017_maa_alueet_sote.gpkg"

# Get data (Finnish municipalities 2017) from Github if not already present
if (!file.exists(dst_file)) {
  download.file("https://github.com/avoindata/mml/raw/master/data/mml/kuntajako_2017_maa_alueet_sote.gpkg",
                dst_file)
}

# Read data. Note that the CRS (Coordinate Reference System) of the
# data is ETRS89 / ETRS-TM35FIN (epsg:3067)
fin_muni <- sf::read_sf("data/kuntajako_2017_maa_alueet_sote.gpkg")

# Dissolve municipalities to get just the Finnish borders
fin <- fin_muni %>%
  # Create a dummy variable which is TRUE everywhere for the dissolve
  mutate(is_fin = TRUE) %>%
  # Dissolve
  group_by(is_fin) %>%
  summarize()

# make_grid() requires a SpatilaPolygonsDataframe object
fin_sp <- sf::as_Spatial(fin)

# Create the hexagon grid. ETRS89 / ETRS-TM35FIN has meters as its unit, so
# requested cell area here is 10 000 x 10 000 = 1e+08 m2
hex_grid <- make_grid(fin_sp, cell_area = 1e+08, clip = FALSE)
# Convert the SpatialPolygonsDataframe back to a sf object
hex_grid <- sf::st_as_sf(hex_grid)

# Write out the hexagon grid
dsts_hex <- c("data/fin_hex_grid.gpkg",
             "data/fin_hex_grid.shp")
for (dst_hex in dsts_hex) {
  if (!file.exists(dst_hex)) {
    sf::write_sf(hex_grid, dsn = dst_hex)
  }
}

# Transform both data sets to WGS84 (lat/lon)
fin_wgs84 <- sf::st_transform(fin, "+init=epsg:4326")
hex_grid_wgs84 <- sf::st_transform(hex_grid, "+init=epsg:4326")

# Create an animated version including both plots
fig <- magick::image_graph(res = 96)

ggplot() +
  geom_sf(data = fin_wgs84, fill = "lightgrey", colour = "lightgrey") +
  geom_sf(data = hex_grid_wgs84, colour = "orange", alpha = 0) +
  ggtitle("WGS84") +
  theme_bw() + theme(legend.position = "none")

ggplot() +
  geom_sf(data = fin, fill = "lightgrey", colour = "lightgrey") +
  geom_sf(data = hex_grid, colour = "orange", alpha = 0) +
  ggtitle("ETRS89 / ETRS-TM35FIN") + 
  theme_bw() + theme(legend.position = "none")
dev.off()

animation <- magick::image_animate(fig, fps = 1)

# Write out the animation
magick::image_write(animation, "data/hexfin_animation.gif")