Population Map Plot

Authors

[Editor] Xuyang Yi;

Example

Shows the geographical distribution of disease incidence, with light and dark colors indicating higher and lower incidence rates, and map boundaries representing administrative divisions.

Setup

System Requirements: Cross-platform (Linux/MacOS/Windows)
Programming language: R
Dependent packages: sf; rnaturalearth; rnaturalearthdata; gstat; dplyr; ggplot2; ggspatial; ggnewscale; ggrepel; ggfx; doParallel; viridis

# Install packages
if (!requireNamespace("sf", quietly = TRUE)) {
  install.packages("sf")
}
if (!requireNamespace("rnaturalearth", quietly = TRUE)) {
  install.packages("rnaturalearth")
}
if (!requireNamespace("rnaturalearthdata", quietly = TRUE)) {
  install.packages("rnaturalearthdata")
}
if (!requireNamespace("gstat", quietly = TRUE)) {
  install.packages("gstat")
}
if (!requireNamespace("dplyr", quietly = TRUE)) {
  install.packages("dplyr")
}
if (!requireNamespace("ggplot2", quietly = TRUE)) {
  install.packages("ggplot2")
}
if (!requireNamespace("ggspatial", quietly = TRUE)) {
  install.packages("ggspatial")
}
if (!requireNamespace("ggnewscale", quietly = TRUE)) {
  install.packages("ggnewscale")
}
if (!requireNamespace("ggrepel", quietly = TRUE)) {
  install.packages("ggrepel")
}
if (!requireNamespace("ggfx", quietly = TRUE)) {
  install.packages("ggfx")
}
if (!requireNamespace("doParallel", quietly = TRUE)) {
  install.packages("doParallel")
}
if (!requireNamespace("viridis", quietly = TRUE)) {
  install.packages("viridis")
}

# Load packages
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(gstat) 
library(dplyr)
library(ggplot2)
library(ggspatial) 
library(ggnewscale)
library(ggrepel)
library(ggfx)
library(doParallel)
library(viridis)

Data Preparation

# Global geographic data
world <- ne_countries(scale = "medium", returnclass = "sf")
# Simulating epidemiological data
set.seed(123)
world$incidence <- runif(nrow(world), 0, 100)  # Randomly generate incidence data

Visualization

1. Basic Plot

# Basic map of global disease incidence distribution
p1 <- ggplot(data = world) +
  geom_sf(aes(fill = incidence)) +
  scale_fill_viridis(option = "C") +
  labs(title = "Global Disease Incidence",
       fill = "Incidence Rate\n(per 100k)")

p1

Figure 1: Basic map of global disease incidence distribution

Tip

Key parameter analysis: binwidth / bins aes(fill): Define color maps for epidemiological indicators

scale_fill_viridis(): Use accessible color gradients

ne_countries(): The scale parameter controls the map’s level of detail.（small/medium/large）

# Customized epidemiological maps
p2 <- ggplot(data = world) +
  geom_sf(aes(fill = incidence), color = "white", size = 0.2) +
  scale_fill_gradientn(colors = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
                       breaks = c(20, 40, 60, 80)) +
  # scale
  annotation_scale(
    location = "bl",
    plot_unit = "km",
    style = "ticks",
    width_hint = 0.1
  ) +
  # Compass (after restoration)
  annotation_north_arrow(
    location = "tr",
    which_north = "grid",  # Use Grid North
    style = north_arrow_minimal(
      line_width = 1,
      text_size = 10
    ),
    pad_x = unit(1.2, "cm"),
    pad_y = unit(1.2, "cm")
  ) +
  coord_sf(crs =  "+proj=robin",
           xlim = c(-1.6e7, 1.6e7),
           ylim = c(-7.5e6, 8.5e6),
           expand = FALSE) +  # Use Robinson projection
  theme_void() +
  guides(fill = guide_colorbar(barwidth = 1)) +
  labs(fill = "Incidence Rate")+ 
  theme(
    legend.position = c(0.1, 0.3)  # Relative position coordinates (the lower left corner is 0,0)
  )

p2

Figure 2: Customized epidemiological maps

2. Advanced plot

# Advanced plot
# Further process the data and use spatial points to infer the situation of nearby surfaces
# Reduce calculation time, the fewer points below, the closer the estimated distance, the higher the resolution, and the more accurate it is
# Parallel initialization
registerDoParallel(cores = 4)

# Data Preparation
world_proj <- st_transform(world, "+proj=eqc +units=m")
centroids <- st_centroid(world_proj) %>% 
  dplyr::select(incidence)

# Create a low-resolution grid (100x100)
grid <- st_make_grid(world_proj, n = c(100,100)) %>% 
  st_as_sf() %>%
  st_join(world_proj, join = st_intersects)

# Variogram model optimization
variogram_model <- vgm(
  psill = 30,
  model = "Exp",  # Switch to an exponential model
  range = 2e6,    # 2000 km related range
  nugget = 5
)

# Block parallel computing
grid_chunks <- split(grid, cut(st_coordinates(st_centroid(grid))[,1], 4))

krige_result <- foreach(i=1:4, .combine=rbind) %dopar% {
  krige(incidence ~ 1,
        locations = centroids,
        newdata = grid_chunks[[i]],
        model = variogram_model,
        nmax = 30)
} %>% st_as_sf()

# Convert back to WGS84 coordinate system
krige_result <- st_transform(krige_result, 4326)


advanced_map <- ggplot() +
  # Spatially interpolated surfaces
  geom_sf(data = krige_result, 
          aes(fill = var1.pred, color = var1.pred),
          alpha = 0.6) +
  scale_fill_gradientn(
    colors = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
    name = "Interpolated\nIncidence"
  ) +
  scale_color_gradientn(
    colors = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
    guide = "none"
  ) +
  ggnewscale::new_scale_fill() +
  
  # Original national borders
  geom_sf(data = world, 
          aes(fill = incidence),
          color = "white", 
          size = 0.1,
          alpha = 0.5) +
  scale_fill_gradientn(
    colors = c("#2c7bb6", "#abd9e9", "#ffffbf", "#fdae61", "#d7191c"),
    name = "Country\nIncidence",
    breaks = seq(0, 100, 20)
  ) +
  
  # 3D relief effect
  ggfx::with_shadow(
    geom_sf(data = world,
            aes(geometry = geometry),
            color = "grey30",
            fill = NA,
            size = 0.2),
    sigma = 5,
    x_offset = 3,
    y_offset = 3
  ) +
  
  # Hotspot annotation
  geom_sf(data = world %>% filter(incidence > quantile(incidence, 0.9)),
          color = "red",
          fill = NA,
          size = 0.5) +
  ggrepel::geom_label_repel(
    data = world %>% filter(incidence > quantile(incidence, 0.9)),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    size = 2.5,
    box.padding = 0.2,
    min.segment.length = 0,
    fill = alpha("white", 0.8)
  ) +
  
  # Map elements
  annotation_north_arrow(
    location = "tr",
    width = unit(1.2, "cm"),   # New width parameter
    height = unit(1.2, "cm"),  # Added height parameter
    style = north_arrow_minimal()  # Remove size parameters
  ) +
  
  # Coordinate projection
  coord_sf(crs = "+proj=robin") +
  
  # Theme Settings
  theme_void() +
  theme(
    legend.position = c(0.12, 0.3),
    legend.background = element_rect(fill = alpha("white", 0.8), color = NA),
    plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
    panel.background = element_rect(fill = "#F0F8FF")
  ) +
  labs(title = "Advanced Spatial Distribution of Disease Incidence")

advanced_map

Optionally, you can use callout-tip to add a detailed description of the parameter if needed.

Application

Population Map Application

This chart shows how common cancers vary across countries around the world. [1]

Reference

[1] Global Burden of Disease 2019 Cancer Collaboration; Kocarnik JM, Compton K, Dean FE, et al. Cancer Incidence, Mortality, Years of Life Lost, Years Lived With Disability, and Disability-Adjusted Life Years for 29 Cancer Groups From 2010 to 2019: A Systematic Analysis for the Global Burden of Disease Study 2019. JAMA Oncol. 2022;8(3):420-444. https://doi.org/10.1001/jamaoncol.2021.6987

[2] Fidler MM, Gupta S, Soerjomataram I, Ferlay J, Steliarova-Foucher E, Bray F. Cancer incidence and mortality among young adults aged 20-39 years worldwide in 2012: a population-based study. Lancet Oncol. 2017;18(12):1579-1589. https://doi.org/10.1016/S1470-2045(17)30677-0

[3] Wickham H, Chang W, Henry L, et al. ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics [Computer software]. (Version 3.4.0). 2022. https://ggplot2.tidyverse.org/reference/ggsf.html