最終更新:2026-03-30

1 はじめに

  • rayshader パッケージを使って東京都(島嶼部除く)の人口分布を3Dビジュアライゼーションとして表現する。
  • Spencer Schien のチュートリアルおよび RayshaderPortraits の手法に準拠している。
  • 人口データを「高さ(山の標高)」に見立てたheightmapを作成し、光・影・色を加えて3Dレンダリングする。
  • データソースとして、世界全体をカバーするKontur Population(400m解像度のH3ヘックスグリッド)を使用する。
  • ラスター化には stars::st_rasterize() を使用する。ベクター形式のヘックスデータを直接ラスターグリッドに変換できる。

参考コード:https://github.com/andrea-io/RayshaderPortraits/tree/main/Japan

1.1 使用データ

データ 説明
Kontur Population(日本) 400m解像度の人口データ(GeoPackage形式)
国土数値情報 行政区域(東京都) 東京都の行政区域ポリゴン(Shapefile形式)

1.2 Kontur Population の仕様

項目 仕様
ファイル形式 GeoPackage (.gpkg)
空間解像度 400m(H3解像度8)
座標参照システム EPSG:3857(Pseudo-Mercator)
集計基準年 2023年
カバレッジ 世界全域

2 データのダウンロード

2.1 Kontur Population(日本)

  1. HDX(Humanitarian Data Exchange)にアクセスする。
  2. 検索ボックスで「Japan」を検索する。
  3. kontur_population_JP_20231101.gpkg(45.8MB)をダウンロードして作業フォルダに配置する。

2.2 国土数値情報 行政区域(東京都)

  1. 国土数値情報ダウンロードサービスにアクセスする。
  2. 「行政区域」→ 東京都 → Shapefile形式でダウンロードする。
  3. 解凍して N03-20250101_13_GML/ フォルダを作業フォルダに配置する。

3 パッケージの読み込み

library(sf)          # ベクターデータの入出力・空間演算
library(stars)       # ラスターデータの操作(st_rasterize を提供)
library(ggplot2)     # データの可視化・地図描画
library(MetBrewer)   # 美術館にインスパイアされたカラーパレット
library(colorspace)  # 色彩ユーティリティ(swatchplot を提供)
library(rgl)         # rayshader が内部で使う3Dグラフィクスバックエンド
library(rayshader)   # heightmap の3Dレンダリング

4 境界データの読み込みとクリップ

# 国土数値情報の行政区域データ(東京都)を読み込む
tokyo <- sf::st_read("N03-20250101_13_GML/N03-20250101_13.shp", quiet = TRUE)

# 市区町村コードで島嶼部(13361以上)を除いた本土部分を抽出する
# 13101(千代田区)〜 13308(西多摩郡檜原村)が本土の範囲
tokyo_main <- tokyo[
  as.numeric(tokyo$N03_007) >= 13101 &
  as.numeric(tokyo$N03_007) <= 13308, ]

# 複数ポリゴンを1つに統合する(同一市区町村の複数ポリゴンを結合)
tokyo_main_union <- sf::st_union(tokyo_main)

# Kontur と同じ座標系(EPSG:3857)に変換する
tokyo_main_3857 <- sf::st_transform(tokyo_main_union, 3857)
# Kontur Population(日本)の GeoPackage を読み込む
JP <- sf::st_read("kontur_population_JP_20231101.gpkg", quiet = TRUE)

# 東京都本土の境界に含まれるヘックスのみ抽出する(空間クリップ)
pop_main <- JP[tokyo_main_3857, ]

cat("Clipped hexagons:", nrow(pop_main), "\n")
## Clipped hexagons: 2265
# データが正しくクリップされたかを概観プロットで確認する
pop_main |>
  ggplot2::ggplot() +
  geom_sf(color = NA, fill = "steelblue", alpha = 0.4) +
  labs(
    title   = "Kontur Population: Tokyo (Mainland)",
    caption = "Source: Kontur Population 2023"
  ) +
  theme_void()

5 バウンディングボックスとアスペクト比の計算

  • バウンディングボックス(bounding box):データ全体を囲む最小の長方形領域のこと。左下・右上の2点(xmin, ymin, xmax, ymax)で定義される。地図の範囲を把握するために使う。
  • アスペクト比(aspect ratio):画像の幅と高さの比率のこと。地図の実際の縦横比に合わせて出力画像を設定しないと、地形が歪んで表示されてしまう。
# クリップ後のデータからバウンディングボックスを取得する
bb <- sf::st_bbox(pop_main)

# 地図の幅と高さをメートル単位で計測するためのコーナーポイントを作成する
bottom_left <- sf::st_point(c(bb[["xmin"]], bb[["ymin"]])) |>
  sf::st_sfc(crs = sf::st_crs(pop_main))   # 左下コーナー

bottom_right <- sf::st_point(c(bb[["xmax"]], bb[["ymin"]])) |>
  sf::st_sfc(crs = sf::st_crs(pop_main))   # 右下コーナー

top_left <- sf::st_point(c(bb[["xmin"]], bb[["ymax"]])) |>
  sf::st_sfc(crs = sf::st_crs(pop_main))   # 左上コーナー

# 東西方向の距離(地図の幅)を計測する
width  <- sf::st_distance(bottom_left, bottom_right)

# 南北方向の距離(地図の高さ)を計測する
height <- sf::st_distance(bottom_left, top_left)

cat("Map width :", round(as.numeric(width)  / 1000, 1), "km\n")
## Map width : 110.1 km
cat("Map height:", round(as.numeric(height) / 1000, 1), "km\n")
## Map height: 51 km
# 出力画像が歪まないようにアスペクト比を計算する
# 長辺を 1 とし、短辺をそれに対する比率で表す
if (width > height) {
  w_ratio <- 1
  h_ratio <- height / width
} else {
  h_ratio <- 1
  w_ratio <- width / height
}

cat("w_ratio:", round(as.numeric(w_ratio), 4), "\n")
## w_ratio: 1
cat("h_ratio:", round(as.numeric(h_ratio), 4), "\n")
## h_ratio: 0.4633

6 人口データのラスター化

  • ヘックスデータ(hex data):地図を正六角形(ヘックス)のタイルで均等に分割し、各タイルに統計値を割り当てたベクターデータのこと。正方形グリッドより隣接セル間の距離が均一になるという利点がある。Kontur Population では H3(Uber が開発した六角形地理インデックス)の解像度8(1セル≒400m)を採用している。
  • ラスター化(rasterization):点・線・ポリゴンなどのベクターデータを、一定間隔のグリッド(画素の集合)に変換する処理のこと。ここでは六角形(ヘックス)ポリゴンの人口値を規則的なグリッドセルに焼き付け、rayshader が扱える行列(matrix)形式にする。
# 長辺のピクセル数を設定する(source.R では 5000;HTMLプレビュー用に 1200 を使用)
size <- 1200

# アスペクト比に基づいて各軸のピクセル数を計算する
nx <- floor(size * as.numeric(w_ratio))  # 列数(東西方向)
ny <- floor(size * as.numeric(h_ratio))  # 行数(南北方向)

cat("Raster size:", nx, "x", ny, "pixels\n")
## Raster size: 1200 x 555 pixels
# stars::st_rasterize でヘックスベクターデータを規則格子ラスターに変換する
# nx / ny が解像度を決める;セル内の人口値は合計される
pop_main_rast <- stars::st_rasterize(pop_main, nx = nx, ny = ny)

# 東京本土の境界ポリゴンを同じグリッドでラスター化してマスクを作成する
# template を指定することで人口ラスターと同一のグリッドになる
tokyo_mask_sf  <- sf::st_sf(mask = 1L, geometry = tokyo_main_3857)
tokyo_mask_rast <- stars::st_rasterize(tokyo_mask_sf, template = pop_main_rast)
mask_mat <- matrix(tokyo_mask_rast$mask, nrow = nx, ncol = ny)

# rayshader が処理できる通常の行列(matrix)に変換する
# nrow = nx:画像の列数=行列の行方向に対応する
mat <- matrix(pop_main_rast$population, nrow = nx, ncol = ny)

# 境界内でデータなし(NA)のセルのみ 0 に置換する(人口ゼロとして扱う)
# 境界外のセルは NA のまま残す(solid = FALSE で透明に描画される)
mat[!is.na(mask_mat) & is.na(mat)] <- 0

cat("Matrix size:", nrow(mat), "rows x", ncol(mat), "cols\n")
## Matrix size: 1200 rows x 555 cols

7 カラーパレットの設定

# MetBrewer の "Hiroshige" パレットを使用する(黄・橙から青系のグラデーション)
color <- MetBrewer::met.brewer(name = "Hiroshige")

# 256色に拡張し、低い値(人口が少ないエリア)の色域を広げる
# bias > 1 にすると低値側の色幅が広がり、高値のピークが際立つ
tx <- grDevices::colorRampPalette(color, bias = 2)(256)

# ベースパレット(8色)をスウォッチで確認する
colorspace::swatchplot(color)

8 3Dレンダリング

  • 3Dレンダリング(3D rendering):数値データ(ここでは人口を高さとした行列)をもとに、光・影・カメラ視点を計算して2次元画像として出力する処理のこと。rayshader はパストレーシング(光線の経路を物理的に追跡する手法)によって高品質な画像を生成する。
  • 3Dシーン(3D scene):カメラ位置・光源・オブジェクト(heightmap)の配置をまとめた仮想空間のこと。plot_3d() でシーンを構築し、render_highquality() で画像として書き出す。

8.1 3Dシーンの構築

# ヘッドレスの rgl デバイスに3Dシーンを構築する
mat |>
  rayshader::height_shade(texture = tx) |>
  rayshader::plot_3d(
    heightmap   = mat,
    zscale      = 10,        # 高さの強調係数(小さいほど山が高くなる)
    solid       = TRUE,     # FALSEだと底面(ソリッド)を非表示にしてピークのみ表示する
    shadowdepth = 0,         # 表面下の影の深さ
    windowsize  = c(nx, ny)  # ウィンドウを出力解像度に合わせる
  )

# カメラ位置を設定する
# theta:垂直軸まわりの回転角(20 = やや右に回転)
# phi  :仰角(30 = 低めの視点で横長の地形を見渡す)
# zoom :ズーム倍率(0.85 = 横長マップ全体を収める)
rayshader::render_camera(theta = 20, phi = 30, zoom = 0.85)

8.2 高品質レンダリング

# パストレーシングによる高品質レンダリングを実行する
# lightdirection = 120   : 南東方向から光を当てる
# lightcolor[1]          : パレットの color[4] を使った暖色系の光
# lightcolor[2] = "white": 上方からの寒色系の補助光
# samples = 100          : 出版品質にするには 450 に増やす(時間がかかる)
rayshader::render_highquality(
  filename       = "output/tokyo_population_3d.png",
  interactive    = FALSE,
  lightdirection = 120,
  lightaltitude  = c(20, 80),
  lightcolor     = c(color[4], "white"),
  lightintensity = c(900, 150),
  samples        = 100,
  width          = nx,
  height         = ny
)

8.3 3Dマップの表示

# レンダリング済みの PNG を HTML ドキュメント内に埋め込む
knitr::include_graphics("output/tokyo_population_3d.png")

Appendix: 2D人口マップ

  • ここでは rayshader を使わず、ggplot2 で東京都(島嶼部除く)の人口分布を2Dコロプレスマップとして描画する。
  • Kontur のヘックスポリゴンをそのまま geom_sf() で描画し、人口を Hiroshige パレットで色付けする。
# 市区町村ポリゴンを EPSG:3857 に変換する(レイヤーの CRS を統一するため)
# tokyo_main は st_union 前の個別ポリゴンなので区市町村境界として使える
tokyo_main_sf_3857 <- sf::st_transform(tokyo_main, 3857)

# Hiroshige パレット(10色)を取得する
palette_cols <- MetBrewer::met.brewer("Hiroshige", n = 10, type = "continuous")

# ggplot2 で人口コロプレスマップを描画する
ggplot2::ggplot() +
  # 下地として東京本土全体のポリゴンを薄いグレーで描く
  ggplot2::geom_sf(
    data      = tokyo_main_sf_3857,
    fill      = "grey92",
    color     = "grey60",
    linewidth = 0.4
  ) +
  # Kontur ヘックスを人口の大きさで色分けして重ねる
  ggplot2::geom_sf(
    data  = pop_main,
    ggplot2::aes(fill = population),
    color = NA                      # ヘックスの境界線は非表示
  ) +
  # 市区町村の境界線を上から重ねて区分を明確にする
  ggplot2::geom_sf(
    data      = tokyo_main_sf_3857,
    fill      = NA,                 # 塗りつぶしなし
    color     = "white",
    linewidth = 0.5
  ) +
  # カラースケールを設定する
  ggplot2::scale_fill_gradientn(
    colors = palette_cols,
    name   = "Population\n(400m grid)",
    labels = scales::comma           # 3桁区切りの数値表示
  ) +
  # タイトル・キャプションを設定する
  ggplot2::labs(
    title    = "Population Distribution: Tokyo (Mainland)",
    subtitle = "Kontur Population 2023 (400m H3 hex grid)",
    caption  = "Data: Kontur Population 2023, MLIT Japan | Author: Ayumu Tanaka"
  ) +
  # 背景・軸を削除してシンプルなレイアウトにする
  ggplot2::theme_void() +
  ggplot2::theme(
    plot.title    = ggplot2::element_text(size = 18, face = "bold",
                                          hjust = 0.5, margin = ggplot2::margin(b = 5)),
    plot.subtitle = ggplot2::element_text(size = 11, hjust = 0.5,
                                          color = "grey40", margin = ggplot2::margin(b = 10)),
    plot.caption  = ggplot2::element_text(size = 8, hjust = 1, color = "grey50"),
    legend.position = "right",
    plot.margin   = ggplot2::margin(20, 20, 20, 20)
  )

参考文献