Rによる天気予報風のグラフ表現

3D棒グラフ＋地図の例

あきる（paithiov909）

2024年10月19日

誰？

今年の夏をふりかえる🎐

バテてたら夏は終わっていた

この夏にやったこと：

台風すごかった🌀

画像：台風10号まとめ　記録的な大雨・暴風　秋も台風の発生しやすい状況続く　動向に注意(気象予報士 吉田　友海 2024年09月02日) - 日本気象協会 tenki.jp

今日やること🎯

こういう表現をRでやりたい

画像：https://x.com/wni_jp/status/1425815449036476420 から

具体的に🥅

積算降水量の3D棒グラフ＋地図を描く

  • 直近48時間の積算降水量（48時間降水量）のグラフ
    • 3D表現（Z=0のXY平面上に地図を描き、Z軸方向に垂直に棒グラフを生やす）
    • 地図上で観測地点がある位置に棒グラフを描く
    • 棒は降水量で階級に分けて色をつける

どうやって❓

実はExcelでできるのだが…

  • CSVデータが気象庁｜最新の気象データあたりから取得できる
    • jmastatsパッケージで行けそうだが、ここでは手動でDLしておいた
  • 地図の表現にはMapbox GL JSというライブラリを使う
    • Rからはmapboxerというパッケージとして利用できる
    • ggplot2（rayshader、ggrglなど？）やplotlyでも地図は描けるが、地図のZ軸方向に棒グラフを立てる表現が難しい

降水量の色分け🚥

とりあえず適当に色分けする関数を用意しておく。 この関数の戻り値は#RRGGBBという形のカラーコード

val2colcode <- \(val) {
  dplyr::case_when(
    val >= 300 ~ "purple",
    dplyr::between(val, 200, 300) ~ "red",
    dplyr::between(val, 150, 200) ~ "yellow",
    dplyr::between(val, 100, 150) ~ "green",
    dplyr::between(val, 50, 100) ~ "blue",
    dplyr::between(val, 1, 50) ~ "cyan",
    .default = "gray"
  ) |>
    # 実際は7色だけ変換すればよいので、
    # この書き方は処理としては無駄が多い
    col2rgb() |>
    t() |>
    rgb(maxColorValue = 255)
}

データの読み込み📖

jmastats::stationsと結合して観測地点の座標を付けておく

pre48h <-
  readr::read_csv(
    here::here("pre48h00_rct.csv"), # 9/1 13時更新分のデータ
    locale = readr::locale(encoding = "Shift_JIS") # Why are you using Shift_JIS??
  ) |>
  dplyr::select(`観測所番号`, `地点`, `現在値(mm)`) |>
  dplyr::rename_with(~ c("station_no", "station_name_label", "val")) |>
  dplyr::filter(!is.na(val)) |> # 欠測は削除する
  dplyr::left_join(
    dplyr::distinct(
      jmastats::stations,
      station_no,
      .keep_all = TRUE
    ),
    by = dplyr::join_by(station_no == station_no)
  ) |>
  dplyr::distinct(geometry, .keep_all = TRUE) # 位置で再度`distinct`する

データの整形🔨

それぞれの観測地点に対応する棒グラフそのものを地物（POLYGON）として用意する。ここでは、観測地点をsf::st_bufferで広げる

pre48h <- pre48h |>
  dplyr::mutate(
    station_name_label = station_name_label,
    height = val * 600, # 地物の高さ。ここでは適当に`600`をかけておく
    colour = val2colcode(val), # さっきの色分け関数
    geometry = geometry,
    .keep = "used"
  ) |>
  sf::st_as_sf() |> # `sf`に変換しなおす
  sf::st_buffer(4000) # 棒グラフになるPOLYGONをつくる

データの確認🔎

pre48h
#> Simple feature collection with 1282 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 122.9385 ymin: 24.01834 xmax: 145.8088 ymax: 45.55648
#> Geodetic CRS:  WGS 84
#> # A tibble: 1,282 × 5
#>    station_name_label           val                       geometry height colour
#>  * <chr>                      <dbl>                  <POLYGON [°]>  <dbl> <chr> 
#>  1 宗谷岬（ソウヤミサキ）       7   ((141.9294 45.48422, 141.9299…   4200 #00FF…
#>  2 稚内（ワッカナイ）          11   ((141.6842 45.4508, 141.6839 …   6600 #00FF…
#>  3 礼文（レブン）              12.5 ((141.0685 45.33698, 141.0675…   7500 #00FF…
#>  4 声問（コエトイ）             5   ((141.7783 45.43538, 141.7776…   3000 #00FF…
#>  5 浜鬼志別（ハマオニシベツ）  21   ((142.1279 45.3135, 142.1292 …  12600 #00FF…
#>  6 本泊（モトドマリ）           4.5 ((141.2332 45.22688, 141.2347…   2700 #00FF…
#>  7 沼川（ヌマカワ）             1.5 ((141.861 45.21264, 141.8618 …    900 #00FF…
#>  8 沓形（クツガタ）             7.5 ((141.1679 45.14849, 141.1695…   4500 #00FF…
#>  9 豊富（トヨトミ）             6.5 ((141.8307 45.0978, 141.8313 …   3900 #00FF…
#> 10 浜頓別（ハマトンベツ）       6   ((142.3343 45.15935, 142.3335…   3600 #00FF…
#> # ℹ 1,272 more rows

ggplot2での地図表現🗾

これを素直にggplot2で描画すると、観測地点が全国各地に散らばっているようすを確認できる

library(ggplot2)

gp <- ggplot(pre48h) +
  geom_sf(aes(fill = factor(colour))) +
  theme_grey() +
  labs(title = "点の色は指定したものとは違うよ！")

gp

ggplot2での地図表現🗾

東海道新幹線のデータ🚅

東海道新幹線が運休して大変だったらしいので、新幹線の路線図を重ねて描いてみたい。jprailwayという自作パッケージからデータを読み込んでおく

# 鉄道区間のsf (LINESTRING)
tkd_line <- jprailway::polylines |>
  dplyr::filter(name == "東海道新幹線")

# 駅のsf (POLYGON). ただし、ここでは代表点の座標を使う
tkd_station <- jprailway::stations |>
  dplyr::semi_join(
    jprailway::lines |>
      dplyr::filter(name == "東海道新幹線") |>
      tidyr::unnest(station_list, names_sep = "_"),
    by = c("code" = "station_list_code")
  ) |>
  dplyr::select(code, name, lat, lng) |>
  sf::st_drop_geometry()

東海道新幹線の路線図🚅

これは、こんな感じのデータ

コード 
gp +
  geom_sf(data = tkd_line) +
  geom_point(aes(x = lng, y = lat), data = tkd_station) +
  ggrepel::geom_label_repel(aes(x = lng, y = lat, label = name), data = tkd_station) +
  coord_sf(xlim = c(135, 140), ylim = c(34.5, 35.8)) +
  labs(
    title = "東海道新幹線の停車駅",
    caption = paste(
      "出典：国土数値情報（鉄道データ）（国土交通省）",
      "https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html",
      sep = "\n"
    )
  )

Mapbox GL JSとは🎁

コード 
library(mapboxer)

jpmap <-
  mapboxer::mapboxer(
    style = basemaps$Carto$voyager,
    center = c(140, 40),
    zoom = 3.2,
    pitch = 60.00,
    bearing = -12.0
  )
コード 
jpmap <- jpmap |>
  mapboxer::add_line_layer(
    id = "tkd_line",
    source = tkd_line |>
      mapboxer::as_mapbox_source(),
    line_color = "#000000",
    line_opacity = 0.8
  ) |>
  mapboxer::add_circle_layer(
    id = "tkd_station",
    source = tkd_station |>
      mapboxer::as_mapbox_source(),
    circle_color = "#000000",
    circle_opacity = 0.8
  ) |>
  mapboxer::add_tooltips(
    layer_id = "tkd_station",
    tooltip = "{{name}}"
  )

本来はこんな雰囲気で地図を描画するためのライブラリ

この地図ウィジェットは、左クリックでのドラッグでパン、右クリックでのドラッグで回転、ホイールで倍率の変更ができる

Mapbox GL JSによる地図🎨

※ウィジェットの高さや幅が小さい場合、ページを再読込みしてみてください

高さがある地物を描く🗼

Mapbox GL JSには、地図にグラフを重ねることができる特別な機能があるわけではないが、fill-extrusionというレイヤーとして追加することで、高さがある地物を描画することができる

↑本来はこういった表現をするための機能らしい

3D棒グラフ＋地図の表現📊

コード 
jpmap <- jpmap |>
  mapboxer::add_source(
    id = "extrusion_src",
    source = pre48h |>
      mapboxer::as_mapbox_source()
  ) |>
  mapboxer::add_layer(
    list(
      id = "extrusion",
      type = "fill-extrusion",
      source = "extrusion_src",
      paint = list(
        "fill-extrusion-color" = list("get", "colour"),
        "fill-extrusion-height" = list("get", "height"),
        "fill-extrusion-base" = 0,
        "fill-extrusion-opacity" = 0.6
      )
    )
  )
コード 
jpmap |>
  mapboxer::add_tooltips(
    layer_id = "extrusion",
    tooltip = "観測所: {{station_name_label}}<br />降水量: {{val}} mm"
  )

こちらのブログ記事のやり方を参考に、棒グラフの棒そのものを地物（POLYGON）としてあらかじめ用意したので、それを描画することで、3D棒グラフのような表現ができる

3D棒グラフ＋地図の表現📊

まとめ🌂

Rで天気予報風のグラフを描けた！！

  • ただし、Mapbox GL JSを使う場合はfill-extrusionでの力技による
  • Rで3D棒グラフ＋地図の表現はなかなか大変
  • でも、変な図を自力で描くとR力の高まりを感じられる

Enjoy!!