Rで解析：背景地図を手軽に取得！44種類の地図タイルを自在に取得「maptiles」パッケージの紹介

解析結果を地図上にプロットすると、データの特徴が伝わりやすくなります。しかし、背景となる地図の用意には手間がかかります。「maptiles」パッケージは、その背景地図の準備を簡単におこなえるパッケージです。sfやterraの空間データを渡すだけで、指定範囲に必要なタイルの選択とダウンロード、一枚のラスタへの合成までを自動でおこなうコマンドが収録されています。OpenStreetMapやStadia、Esri、CARTO、Thunderforestなど、44種類の背景地図から選択が可能です。また、得られたデータはterraパッケージのSpatRasterなので、観測地点や調査範囲の重ね書きも手軽におこなえます。出典表記の取得や、国土地理院の地理院タイルのような独自のタイルサービスの登録も可能です。本パッケージの利用で、地図表現の手間を減らせるのではないかと考えます。

パッケージバージョンは0.12.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

誰も教えてくれなかった　医療統計の使い分け〜迷いやすい解析手法の選び方が，Rで実感しながらわかる！ | 吉田 寛輝, 山本 紘司

Amazonで吉田 寛輝, 山本 紘司の誰も教えてくれなかった　医療統計の使い分け〜迷いやすい解析手法の選び方が，Rで実感しながらわかる！。アマゾンならポイント還元本が多数。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("maptiles")

# パッケージの読み込み
library("maptiles")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

・地図タイルとは

世界地図をズームレベルごとに正方形の画像へ分割したもので、必要な範囲の画像だけを取り寄せて並べることで地図を表示する仕組みです。

・CARTOとは

位置情報データの分析と可視化を提供するクラウドサービスで、背景地図用のタイルも公開しています。

基本の流れは3ステップです。

1. get_providersコマンドで利用できるプロバイダを確認
2. sfやterraの空間データをget_tilesコマンドへ渡してタイルを取得
3. plot_tilesコマンドで描画

以下の例は、北海道の一部地域に6地点の観測点を置き、その観測値を付与したものです。

# 空間データを扱うためsfパッケージを読み込み
# install.packages("sf")
library("sf")

# 観測地点データを作成（地点名、経度、緯度、観測値）
obs_data <- data.frame(
  station = c("St_A", "St_B", "St_C", "St_D", "St_E", "St_F"),
  lon = c(141.36, 141.58, 141.25, 141.47, 141.19, 141.62),
  lat = c(42.82, 42.90, 42.75, 42.68, 42.95, 42.73),
  value = c(18.4, 22.1, 15.7, 27.3, 12.9, 20.6)
)

# 経緯度をWGS84（EPSG:4326）の点データとしてsfオブジェクトへ変換
obs_sf <- st_as_sf(obs_data, coords = c("lon", "lat"), crs = 4326)

# 距離をメートルで扱えるようUTM座標系（EPSG:32654）へ投影変換
obs_sf <- st_transform(obs_sf, crs = 32654)

# 全観測地点を囲む矩形に8kmの余白を加えて調査範囲を作成
area_sfc <- st_buffer(st_as_sfc(st_bbox(obs_sf)), dist = 8000)

利用できるタイルプロバイダの一覧取得：get_providersコマンド

オプションはなし。組み込みプロバイダの名前、URL、サブドメイン、出典表記をまとめたリストを取得できます。

リストの要素名がproviderに指定できる名前です。qがタイルのURL、citが出典表記です。出典表記は地図を利用する際に明記することが大切です。

# 組み込みプロバイダの一覧をリストで取得
provider_list <- get_providers()

# 登録されているプロバイダの総数を確認
length(provider_list)
[1] 44

# プロバイダ名からCARTO系のものだけを抽出
grep("CartoDB", names(provider_list), value = TRUE)
[1] "CartoDB.Positron"             "CartoDB.PositronNoLabels"     "CartoDB.PositronOnlyLabels"   "CartoDB.DarkMatter"          
[5] "CartoDB.DarkMatterNoLabels"   "CartoDB.DarkMatterOnlyLabels" "CartoDB.Voyager"              "CartoDB.VoyagerNoLabels"     
[9] "CartoDB.VoyagerOnlyLabels" 

# CartoDB.Positronの設定内容を確認
provider_list$CartoDB.Positron
$src
[1] "CartoDB.Positron"

$q
[1] "https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png"

$sub
[1] "a" "b" "c" "d"

$cit
[1] "© OpenStreetMap contributors © CARTO"

タイル提供元の出典表記の取得：get_creditコマンド

該当するプロバイダが見つからない場合はNULLが返ります。タイルの多くは出典表記の掲載が利用条件になっていますので、地図を自分の範囲外で公開する際は必ず明記してください。

# 組み込みプロバイダ名を指定して出典表記を取得
get_credit("OpenStreetMap")
[1] "© OpenStreetMap contributors"

# 背景地図に使用するCARTOのタイルの出典表記を取得
get_credit("CartoDB.Positron")
[1] "© OpenStreetMap contributors © CARTO"

# 存在しないプロバイダ名を指定した場合の挙動を確認
get_credit("NoSuchProvider")
NULL

空間データの範囲に合わせたタイル取得：get_tilesコマンド

verboseをTRUEにすると、決定されたズームレベルと出典表記、キャッシュ先が表示されます。キャッシュ先は既定でtempdirフォルダなのでRを終了すると消えますが、cachedirにフォルダを指定しておくと、次回以降はダウンロードなしで同じタイルを再利用できます。ズームレベルを1上げるとタイル枚数は約4倍になりますので、サーバへの負荷を考えて必要最小限の値を選んでください。

# 調査範囲に合わせてCARTOの淡色タイルを取得
bg_tile <- get_tiles(
  x = area_sfc,
  provider = "CartoDB.Positron",
  zoom = 10,
  crop = TRUE,
  verbose = TRUE
)

# 返り値のクラスを確認
class(bg_tile)

# タイルをプロット
plot_tiles(bg_tile)

取得したタイルの描画：plot_tilesコマンド

plot_tilesコマンドはterraパッケージのplotRGBコマンドのラッパーです。描画後の座標系は投影後の描写に揃うため、同じ座標参照系のsfオブジェクトであればadd = TRUEでそのまま重ねられます。adjustをTRUEにすると余白の追加や拡大が行われ、タイルが引き伸ばされません。

# 取得したタイルを背景地図として描画
plot_tiles(bg_tile, adjust = TRUE)

# 観測値の大きさを点の直径に反映して観測地点を重ね描き
plot(
  st_geometry(obs_sf),
  pch = 21,
  bg = "#D7263D",
  col = "white",
  cex = sqrt(obs_sf$value) / 2,
  add = TRUE
)

# 地点名を各点の上側に配置
text(
  st_coordinates(obs_sf),
  labels = obs_sf$station,
  pos = 3,
  cex = 0.8,
  font = 2
)

# 出典表記を図の右下に追加
mtext(
  get_credit("CartoDB.Positron"),
  side = 1,
  line = 0,
  adj = 0.9,
  cex = 2
)

独自タイルプロバイダの作成：create_providerコマンド

組み込みプロバイダには日本の地図サービスが含まれていません。ここでは国土地理院が公開している地理院タイル（淡色地図）を独自プロバイダとして登録し、背景地図として使ってみます。実際に利用する際は「地理院タイル利用規約」を確認のうえ、出典の表示を必ず行ってください。

get_creditコマンドは文字列だけでなく、create_providerコマンドが返したオブジェクトも適応可能です。独自プロバイダでも、出典表記を新たに記述する必要はありません。

# 地理院タイル（淡色地図）を独自プロバイダとして定義
gsi_pale <- create_provider( name = "GSI_pale",
                             url = "https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png",
                             citation = "地理院タイル（国土地理院）" )

# 地理院タイルを取得
gsi_tile <- get_tiles( x = area_sfc, provider = gsi_pale, zoom = 11, crop = TRUE )

# 現在の作図設定を退避
old_par <- par(no.readonly = TRUE)

# 取得したタイルを描画
plot_tiles(gsi_tile, adjust = FALSE)

# 調査範囲の外周を重ねて表示
plot(area_sfc, border = "#0B7285", lwd = 2, col = NA, add = TRUE)

# 作図デバイスに依存しない基準として背景地図の範囲を取得
map_ext <- as.vector(terra::ext(gsi_tile))
x_rng <- map_ext[2] - map_ext[1]
y_rng <- map_ext[4] - map_ext[3]

# 作図領域の実寸（インチ）から文字サイズの基準を決定
cex_base <- max(0.6, min(1.1, min(par("pin")) / 4.5))

# 観測地点の座標を行列として取得
obs_xy <- st_coordinates(obs_sf)

# 最大の棒が地図の高さの18%になるよう換算係数と棒の幅を地図範囲から計算
bar_scale <- y_rng * 0.18 / max(obs_sf$value)
bar_width <- x_rng * 0.015

# 各観測地点に観測値の大きさを高さとした棒を描画
rect( xleft = obs_xy[, 1] - bar_width, ybottom = obs_xy[, 2],
      xright = obs_xy[, 1] + bar_width, ytop = obs_xy[, 2] + obs_sf$value * bar_scale,
      col = "#D7263DCC", border = "white" )

# 棒の先端に観測値を表示
text( x = obs_xy[, 1], y = obs_xy[, 2] + obs_sf$value * bar_scale,
      labels = obs_sf$value, pos = 3, cex = cex_base * 0.8, font = 2 )

# インセットの地点名とタイトルの文字サイズを設定
cex_lab <- cex_base * 0.7
cex_ttl <- cex_base * 0.85

# 文字の実寸を地図座標の単位で測定
lab_w <- max(strwidth(obs_sf$station, cex = cex_lab))
lab_h <- strheight("A", cex = cex_lab)
ttl_w <- strwidth("地点別の観測値", cex = cex_ttl)
ttl_h <- strheight("A", cex = cex_ttl)

# 文字寸法から棒1本分の幅と外枠の余白を決定
step <- lab_w * 1.5
pad <- lab_h * 0.9

# 中身が収まる外枠の幅を算出
box_w <- max(step * nrow(obs_sf) + pad * 2, ttl_w + pad * 2)

# 外枠が地図幅の45%を超える場合の縮小率を求める
shrink <- min(1, x_rng * 0.45 / box_w)

# 縮小率を文字サイズと各寸法へ反映
cex_lab <- cex_lab * shrink
cex_ttl <- cex_ttl * shrink
lab_h <- lab_h * shrink
ttl_h <- ttl_h * shrink
step <- step * shrink
pad <- pad * shrink
box_w <- box_w * shrink

# 棒グラフ部分の高さと外枠全体の高さを幅から算出
bar_area_h <- box_w * 0.42
box_h <- pad + lab_h * 1.4 + bar_area_h + ttl_h * 1.8

# 外枠を地図の右下へ配置
box_x2 <- map_ext[2] - x_rng * 0.02
box_x1 <- box_x2 - box_w
box_y1 <- map_ext[3] + y_rng * 0.02
box_y2 <- box_y1 + box_h

# インセットの背景を半透明の白で描画
rect(box_x1, box_y1, box_x2, box_y2,
     col = "#FFFFFFE6", border = "grey50")

# 棒グラフの基線の高さと各棒の中心位置を計算
base_y <- box_y1 + pad + lab_h * 1.4
bar_cx <- box_x1 + (box_w - step * nrow(obs_sf)) / 2 + step * (seq_len(nrow(obs_sf)) - 0.5)

# 観測値をインセットの高さへ換算
in_scale <- bar_area_h / max(obs_sf$value)

# インセットに全地点の観測値を棒グラフで描画
rect( xleft = bar_cx - step * 0.3, ybottom = base_y,
      xright = bar_cx + step * 0.3, ytop = base_y + obs_sf$value * in_scale,
      col = "#D7263D", border = NA )

# 棒グラフの基線を描画
segments(box_x1 + pad, base_y, box_x2 - pad, base_y, col = "grey40")

# 各棒の下に地点名を表示
text(bar_cx, base_y - lab_h * 0.35, labels = obs_sf$station,
     adj = c(0.5, 1), cex = cex_lab)

# インセットのタイトルを上端に表示
text( x = (box_x1 + box_x2) / 2, y = box_y2 - ttl_h * 0.6,
      labels = "地点別の観測値", adj = c(0.5, 1), cex = cex_ttl, font = 2 )

# プロバイダオブジェクトから出典表記を取り出して図に追加
mtext(get_credit(gsi_pale), side = 1, line = 0, adj = 0.9, cex = cex_base)

# 作図設定を元に戻す
par(old_par)

＜おすすめのRに関する書籍です＞

Rをはじめよう生命科学のためのRStudio入門

AmazonでAndrew P. Beckerman, Dylan Z. Childs, Owen L. Petchey, 富永 大介のRをはじめよう生命科学のためのRStudio入門。アマゾンならポイント還元本が多数。

www.amazon.co.jp

この記事が誰かの役に立ちますように。