Grilă hexagonală?
Vom face binning statistic pe o grilă hexagonală, dar nu doar orice grila. Sisteme Geodezice Discrete Global Grid (Kimerling et al. 1999; Sahr, White și Kimerling 2003) permite utilizarea hexagonului ierarhic cu arii egale1 grile.
Pachetul {dggridR} (Barnes 2017) ne va ajuta să generăm aceste grile în Franța.
library(dggridR) library(dplyr) library(purrr) library(sf) library(units) library(glue) library(tibble) library(ggplot2) library(ggspatial) library(rnaturalearth) # devtools::install_github("ropensci/rnaturalearthhires")
Mai întâi vom construi aceste rețele globale discrete pentru Franța metropolitană cu dimensiunea celulei de 1, 5, 10, 20 și 100 km. Avem nevoie de o amprentă spațială.
# get it from # https://static.data.gouv.fr/resources/admin-express-cog-simplifiee-2024-metropole-drom-saint-martin-saint-barthelemy/20240930-094021/adminexpress-cog-simpl-000-2024.gpkg fx <- read_sf("~/data/adminexpress/adminexpress-cog-simpl-000-2024.gpkg", layer = "region") |> filter(insee_reg > "06") |> summarise()
Putem alege ce rezoluție a grilei vom face cu „dimensiunea” aproximativă a hexagonului dorit în km; eșantionarea ne permite să avem toate celulele (dacă sunt prea mari, unele celule vor fi absente). Durează câteva minute…
res <- tribble(~distance, ~sampling, 1, 0.005, 5, 0.010, 10, 0.010, 20, 0.010, 100, 0.100) #' Build a DGG #' #' Discrete Global Grid #' #' @param src (sf) : footprint #' @param dest (char) : output filename #' @param distance (num) : approximate cell size (km) #' @param sampling (num) : points sampling to create the cells (°) #' @param desc (char) : métadata (ex: layer description) #' #' @return (sf) : layer (+ file on disk with metadata) build_dgg <- function(src, dest, distance = 100, sampling = 0.1, desc = "") { msg <- capture.output( dg <- dgconstruct(spacing = distance, projection = "ISEA", aperture = 3, topology = "HEXAGON")) properties <- paste(glue_data(enframe(dg), "{name}: {value}"), collapse = "n") hex <- dg |> dgshptogrid(src, cellsize = sampling) hex |> st_join(src, left = FALSE) |> st_write(dest, layer = glue("dggrid_{distance}k"), layer_options = c( glue("IDENTIFIER=Discrete Global Grid - {distance} km"), glue("DESCRIPTION=ISEA3H {desc} {msg} {properties} {Sys.Date()}")), delete_layer = TRUE) } # execute for all resolutions res |> pwalk(build_dgg, src = fx, dest = "dggrid_fx.gpkg", desc = "France métropolitaine WGS84", .progress = TRUE)
Acum avem un geopachet cu toate grilele noastre:
st_layers("dggrid_fx.gpkg")
Driver: GPKG Available layers: layer_name geometry_type features fields crs_name 1 dggrid_1k Polygon 466080 1 WGS 84 2 dggrid_5k Polygon 17801 1 WGS 84 3 dggrid_10k Polygon 6079 1 WGS 84 4 dggrid_20k Polygon 2106 1 WGS 84 5 dggrid_100k Polygon 102 1 WGS 84
Și putem vedea frumosul cuib la scară spațială în Figura 1.

Figura 1: Eșantioane la 100, 20 și 10 km
Pentru a demonstra binning-ul vom folosi comuna populație disponibilă în GPKG descărcat mai devreme pe care o vom folosi ca strat punct.
Mai întâi avem nevoie de mai multe date…
# output projection # we chose an equal-area projection proj <- "EPSG:3035" fx_laea <- fx |> st_transform(proj) # population data pop <- read_sf("~/data/adminexpress/adminexpress-cog-simpl-000-2024.gpkg", layer = "commune") |> filter(insee_reg > "06") |> st_transform(proj) |> st_point_on_surface() # map zoom fx_bb <- st_bbox(st_buffer(fx_laea, 0)) # projection name lib_proj <- st_crs(fx_laea)$Name # regional boundaries reg_int <- read_sf("~/data/adminexpress/adminexpress-cog-simpl-000-2024.gpkg", layer = "region_int") |> st_transform(proj) # european base map eu <- ne_countries(continent = "europe", scale = 10, returnclass = "sf") |> st_transform(proj) |> st_intersection(st_buffer(fx_laea, 500000))
Apoi creăm o funcție pentru a genera hărțile cu rezoluții diferite ale grilei.
#' Create a map of population for a grid resolution #' #' @param resolution (int) : the grid resolution (among those available in #' dggrid_fx.gpkg) #' #' @return (ggplot) : population map build_map <- function(resolution) { dggrid <- read_sf("dggrid_fx.gpkg", layer = glue("dggrid_{resolution}k")) |> st_transform(proj) |> st_intersection(fx_laea) pop |> st_join(dggrid) |> st_drop_geometry() |> summarise(.by = seqnum, pop = sum(population, na.rm = TRUE)) |> left_join(x = dggrid, y = _, join_by(seqnum)) |> mutate(dens = drop_units(pop / set_units(st_area(geom), km^2))) |> ggplot() + geom_sf(data = eu, color = "#eeeeee", fill = "#eeeeee") + geom_sf(aes(fill = dens, color = dens)) + geom_sf(data = reg_int, color = "#ffffff", linewidth = 0.4, alpha = 0.5) + geom_sf(data = reg_int, color = "#555555", linewidth = .2) + scale_fill_fermenter(name = "inhabitants/km²", palette = "Greens", na.value = "white", breaks = c(20, 50, 100, 500), direction = 1) + scale_color_fermenter(name = "inhabitants/km²", palette = "Greens", na.value = "white", breaks = c(20, 50, 100, 500), direction = 1) + annotation_scale(location = "bl", height = unit(.15, "cm"), line_col = "#999999",text_col = "#999999", bar_cols = c("#999999", "white")) + annotation_north_arrow(pad_x = unit(.25, "cm"), pad_y = unit(.8, "cm"), style = north_arrow_fancy_orienteering( text_col = "#999999", text_size = 8, line_col = "#999999", fill = "#999999"), which_north = "true", height = unit(.5, "cm"), width = unit(.5, "cm")) + labs(title = "Population", subtitle = "Metropolitan France - 2022", caption = glue("data : Discrete Global Grid ISEA3H ≈{resolution} km, Natural Earth and IGN Admin Express 2024 proj. : {lib_proj}")) + coord_sf(xlim = fx_bb(c(1, 3)), ylim = fx_bb(c(2, 4))) + theme_void() + theme(text = element_text(family = "Marianne"), plot.caption = element_text(size = 6), plot.caption.position = "plot", panel.background = element_rect(color = "#E0FFFF", fill = "#E0FFFF55")) }
Acum putem demonstra diferitele noastre rezoluții:
build_map(20)


Figura 2: Populația pe dimensiunea rețelei 20 km
build_map(100)


Figura 3: Populația pe dimensiunea rețelei 100 km