R 语言｜另类地图可视化｜好玩的 rayshader 包制作

R 语言里，一种比较另类的地图可视化，使用柱图代表数据的数量。使用好玩的 rayshader 包制作。

library(tidyverse)
library(rnaturalearth)
library(sf)
library(rayrender)

1. Map

1.1 Get maps

map.sf <- rnaturalearth::ne_countries(
  returnclass = "sf",
  scale = "small"
)

na.map.sf <- rnaturalearth::ne_countries(
  continent = "north america",
  returnclass = "sf",
  scale = "large"
)

ca.map.sf <- rnaturalearth::ne_states(
  country =  "canada",
  returnclass = "sf"
)

usa.map.sf <- rnaturalearth::ne_states(
  country = "united states of america",
  returnclass = "sf"
) |> 
  filter(name != "Alaska")

sts.map.sf <- rnaturalearth::ne_states(
  country = "united states of america",
  returnclass = "sf"
) |> 
  filter(name %in% c("New Mexico", "Colorado", "Wyoming"))

1.2 The bounding

the.box <- sts.map.sf |> st_bbox()
r.long <- seq(the.box["xmin"], the.box["xmax"], len = 30)
r.lat <- seq(the.box["ymin"], the.box["ymax"], len = 10)
coords.sf <- expand.grid(x = r.long, y = r.lat) |>
  st_as_sf(coords = c("x", "y"), crs = "WGS84") |>
  st_intersection(sts.map.sf)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

1.3 Generating data

proj0 <- "+proj=geos +h=358000000 +sweep=y" # https://proj.org/en/9.4/operations/projections
proj.add <- "+lon_0=-100 +lat_0=0"
proj1 <- paste(proj0, proj.add)

set.seed(1)
coords.dat <- coords.sf |>
  st_coordinates() |>
  as.data.frame() |>
  sample_frac(0.5) |> 
  mutate(
    X = X,
    Y = Y
  ) |> 
  mutate(bar.height = rpois(n(), 5) / 5)

1.4 2D map (.png for rayrender)

当然可以直接在 rayrender 包里，制作网格线等。但是我认为那会更复杂（并且可能增加 rayrender 的渲染时间）

故而，直接使用 ggplot2 绘制网格地图，然后贴到 rayrender 渲染的地球的表面

sts.map.bx <- st_bbox(sts.map.sf)

wd <- (sts.map.bx[c(1, 3)] * 100) |> diff(); wd

##     xmax 
## 901.1816

ht <- (sts.map.bx[c(2, 4)] * 100) |> diff(); ht

##     ymax 
## 1367.271

img.bx <- st_bbox(sts.map.bx) |> as.numeric()
long.lat.gap <- 5

xy.lab <- expand.grid(
  x = seq(-125, -95, long.lat.gap),
  y = seq(25, 50, long.lat.gap)
)

p.sf0 <- ggplot() +
  geom_sf(
    data = map.sf,
    fill = NA,
    color = "grey0",
    lwd = 0.05
  ) +
  geom_sf(
    data = sts.map.sf,
    color = "grey0",
    fill = NA,
    lwd = 0.05
  ) +
  geom_sf(
    data = coords.dat |> st_as_sf(coords = c("X", "Y"), crs = "WGS84"),
    shape = 16,
    size = 0.01,
    color = "dodgerblue"
  ) +
  geom_sf_text(
    data = sts.map.sf,
    aes(label = name),
    size = 0.5,
    family = "Smiley Sans"
  ) +
  annotate(
    "segment",
    x = seq(-180, 180, long.lat.gap),
    xend = seq(-180, 180, long.lat.gap),
    y = -90,
    yend = 90,
    color = "grey50",
    lwd = 0.01,
    lty = 1
  ) +
  annotate(
    "segment",
    x = -180,
    xend = 180,
    y = seq(-90, 90, long.lat.gap),
    yend = seq(-90, 90, long.lat.gap),
    color = "grey50",
    lwd = 0.01,
    lty = 1
  ) +
  annotate(
    "text",
    x = xy.lab$x,
    y = xy.lab$y,
    label = paste0("(", xy.lab$x, ", ", xy.lab$y, ")"),
    size = 0.5,
    color = "grey50",
    family = "Smiley Sans"
  ) +
  coord_sf(
    crs = "+proj=longlat",
    xlim = c(-180, 180),
    ylim = c(-90, 90),
    clip = "off",
    expand = FALSE
  ) +
  theme_void()
p.sf0

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

# ggview::ggview(p.sf0, width = 10, height = 5, bg = "white")
ggsave("map.sf.png", p.sf0, width = 16, height = 8, dpi = 900, bg = "white")

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

1.5 2D map (for comparison)

# 各种投影方式
proj0 <- "+proj=geos +h=358000000 +sweep=y" # https://proj.org/en/9.4/operations/projections
proj.add <- "+lon_0=-120 +lat_0=0"
proj1 <- paste(proj0, proj.add)

coords.trans1 <- coords.dat |>
  st_as_sf(coords = c("X", "Y"), crs = "WGS84") |>
  st_transform(crs = proj1) |>
  st_coordinates() |>
  as.data.frame() |> 
  mutate(bar.height = coords.dat$bar.height * 5e4)

p.sf1 <- ggplot() +
  geom_sf(
    data = sts.map.sf,
    color = "grey50",
    fill = NA,
    lwd = 0.5
  ) +
  geom_sf(
    data = coords.dat |> st_as_sf(coords = c("X", "Y"), crs = "WGS84"),
    shape = 16,
    size = 1,
    color = "dodgerblue"
  ) +
  geom_sf_text(
    data = sts.map.sf,
    aes(label = name),
    color = "grey50",
    family = "Smiley Sans"
  ) +
  geom_segment(
    data = coords.trans1,
    aes(
      x = X,
      xend = X,
      y = Y,
      yend = Y + bar.height,
    ),
    color = "dodgerblue",
    lwd = 0.5
  ) +
  coord_sf(
    crs = proj1,
    clip = "off",
    expand = FALSE
  ) +
  theme_void() +
  theme(
    plot.margin = margin(10, 0, 10, 0)
  )
p.sf1