Grile discrete globale | R-bloggeri

URMĂREȘTE-NE
16,065FaniÎmi place
1,142CititoriConectați-vă

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.

Harta mostrelor de grilă

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)
Harta grupării statistice pe o grilă hexagonală a populației francezeHarta grupării statistice pe o grilă hexagonală a populației franceze

Figura 2: Populația pe dimensiunea rețelei 20 km

build_map(100)
Harta grupării statistice pe o grilă hexagonală a populației francezeHarta grupării statistice pe o grilă hexagonală a populației franceze

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

Dominic Botezariu
Dominic Botezariuhttps://www.noobz.ro/
Creator de site și redactor-șef.

Cele mai noi știri

Pe același subiect

LĂSAȚI UN MESAJ

Vă rugăm să introduceți comentariul dvs.!
Introduceți aici numele dvs.