からだにいいもの Rで解析：複数のパッケージでAPIキーなしで地図をプロットする！
「maptiles」パッケージは、OpenStreetMapなどの背景地図を取り寄せ、手元の地点を重ねてプロットすることができます。渡すものは、経度と緯度の列を持つデータフレームです。APIキーの登録や課金の設定は要りません。この記事では、地名から緯度経度を調べる「osmdata」パッケージと、2地点の道のりを求める「osrm」パッケージも合わせて、地図を1枚仕上げるまでをまとめました。手元の住所や座標を地図の上で確かめたい人向けです。
パッケージバージョンは「maptiles」が0.12.0、「osmdata」が0.4.1、「osrm」が5.0.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。
＜おすすめのRに関する書籍です＞
パッケージのインストール
下記コマンドを実行してください。
# パッケージのインストール
install.packages(c("maptiles", "sf", "osmdata", "osrm"))コマンド例
詳細はコメント、パッケージのヘルプを確認してください。
略語＝日本語の意味：
・タイル：地図を正方形の画像に切り分けたもの。必要な範囲の画像だけを取り寄せて並べると1枚の地図になる
・API：Application Programming Interfaceの略で、外部のサービスへプログラムから問い合わせる窓口
・OSM：OpenStreetMapの略で、誰でも編集できる世界の地図データ
・EPSG：座標系に付けられた番号。4326は世界測地系の緯度経度、6674は日本の平面直角座標系第VI系で近畿を担当する
以前は「ggmap」パッケージから、Google MapとStamen Mapのタイルを取得できました。現在はGoogle MapがAPIキーの登録と課金が必要になり、Stamen MapのタイルもStadia Mapsへ移ってキーが要るようになっています。ここではキーの要らないOpenStreetMapのタイルを使い、地点と経路を地図に載せるまでを順に紹介していきます。
まず、例で共通して使うデータを作ります。京都府宇治市にある茶の直売所4店の位置と、月の販売量をまとめた想定の架空データです。
# パッケージの読み込み
library("maptiles")
library("sf")
Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library("osmdata")
Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
connected to: https://overpass-api.de/api/interpreter
library("osrm")
Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
Routing: OSRM - http://project-osrm.org/
# 宇治市内にある茶の直売所4店の位置と月の販売量（kg）
Chokubaijo <- data.frame(
Mise = c("宇治橋通り店", "平等院表参道店", "JR宇治駅前店", "京阪宇治駅前店"),
Keido = c(135.8065, 135.8077, 135.8020, 135.8113),
Ido = c(34.8916, 34.8896, 34.8879, 34.8930),
Hanbai = c(320, 480, 260, 180)
)
# 経度と緯度の列を位置情報に変換：st_as_sfコマンド
# 経度、緯度の順に列名を指定：coordsオプション
# 世界測地系の緯度経度を指定：crsオプション
ChokubaijoSf <- st_as_sf(Chokubaijo, coords = c("Keido", "Ido"), crs = 4326)
# 中身を確認
ChokubaijoSf
Simple feature collection with 4 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 135.802 ymin: 34.8879 xmax: 135.8113 ymax: 34.893
Geodetic CRS: WGS 84
Mise Hanbai geometry
1 宇治橋通り店 320 POINT (135.8065 34.8916)
2 平等院表参道店 480 POINT (135.8077 34.8896)
3 JR宇治駅前店 260 POINT (135.802 34.8879)
4 京阪宇治駅前店 180 POINT (135.8113 34.893)
＜おすすめのRに関する書籍です＞
地名から緯度経度を調べる：getbbコマンド
手元に住所や地名しかないときは、osmdataパッケージのgetbbコマンドで緯度経度を調べることが可能です。取得はOpenStreetMapのNominatimで、APIキーは要りません。既定では市町村などの行政界を探すので、寺社や駅のような場所を探すときはfeaturetypeオプションにNULLを渡します。なお、サーバー負荷軽減のため、連続作業時の情報取得は1秒に1回までがおすすめです。
# 市町村の名前から範囲の緯度経度を取得：getbbコマンド
UjiHanni <- getbb("京都府宇治市")
# 西端、東端、南端、北端の順に確認
UjiHanni
min max
x 135.75965 135.87992
y 34.85792 34.95737
# 施設の名前から代表点の緯度経度を取得：format_outオプション
# 行政界以外も探す：featuretypeオプション
Byodoin <- getbb("平等院 宇治市", format_out = "data.frame", featuretype = NULL)
# 見つかった場所の名前と緯度経度を確認
Byodoin[1, c("display_name", "lat", "lon")]
display_name lat
1 平等院, あじろぎの道, 宇治, 宇治市, 京都府, 611-0021, 日本 34.8895001
lon
1 135.8074525| オプション | 意味 | 初期値 |
|---|---|---|
| place_name | 探す地名、またはWikidataのID。WikidataのIDのときはformat_out、base_url、silentだけが使われる | なし |
| display_name_contains | Nominatimが返すdisplay_nameに含まれる文字列で絞り込む | NULL |
| viewbox | 探す範囲。”x1,y1,x2,y2″の文字列かc(x1, y1, x2, y2)で指定し、xが経度、yが緯度 | NULL |
| format_out | 返す形式。matrixが既定で、stringは文字列、data.frameは該当したすべて、sf_polygonとpolygonは範囲の多角形、osm_type_idはOSMの種別と番号 | c(“matrix”, “data.frame”, “string”, “polygon”, “sf_polygon”, “osm_type_id”) |
| base_url | 問い合わせ先のURL | “https://nominatim.openstreetmap.org” |
| featuretype | 探すOSMの地物の種類 | “settlement” |
| limit | 返す件数の上限 | 10 |
| key | キーを求めるサービスへ渡すAPIキー | NULL |
| silent | TRUEで問い合わせ先を画面に出さない | TRUE |
背景地図のタイルを取得する：get_tilesコマンド
get_tilesコマンドは、渡した範囲を覆うタイルを取り寄せ、1枚につなげた画像を返します。範囲は緯度経度のままでも渡せますが、日本付近では東西が縮んで見えるため、ここでは平面直角座標系第VI系（EPSG:6674）へ直してから渡しています。この座標系の単位はメートルなので、店の周り400メートルの書き方が使用できます。
# 緯度経度を平面直角座標系第VI系へ変換：st_transformコマンド
ChokubaijoHeimen <- st_transform(ChokubaijoSf, 6674)
# 4店それぞれの周り400mを囲む：st_bufferコマンド
# 全体を囲む四角形を作る：st_bboxコマンド
HyojiHanni <- st_as_sfc(st_bbox(st_buffer(ChokubaijoHeimen, 400)))
# 背景地図のタイルを取得：get_tilesコマンド
# タイルの種類を指定：providerオプション
# 拡大の程度を指定：zoomオプション；数字が大きいほど細かい
# 指定した範囲で切り落とす：cropオプション
UjiHaikei <- get_tiles(HyojiHanni, provider = "OpenStreetMap", zoom = 15, crop = TRUE)
# 取り寄せた画像の大きさと座標系を確認
UjiHaikei
class : SpatRaster
size : 350, 423, 3 (nrow, ncol, nlyr)
resolution : 3.91, 3.91 (x, y)
extent : -18498.38, -16844.45, -123757.7, -122389.2 (xmin, xmax, ymin, ymax)
coord. ref. : JGD2011 / Japan Plane Rectangular CS VI (EPSG:6674)
source(s) : memory
colors rgb : 1, 2, 3
names : lyr.1, lyr.2, lyr.3
min values : 10, 0, 0
max values : 255, 255, 255| オプション | 意味 | 初期値 |
|---|---|---|
| x | 範囲を示すsf、sfc、bbox、SpatRaster、SpatVector、SpatExtentのオブジェクト。SpatExtentのときは世界測地系の緯度経度（EPSG:4326）で指定する | なし |
| provider | タイルを取り寄せるサーバ。組み込みの種類の名前か、create_providerコマンドで作った名前付きリストを指定 | “OpenStreetMap” |
| zoom | 拡大の程度 | なし |
| crop | TRUEでxの範囲に合わせて切り落とす。xが1点のsfのときはFALSEになる | FALSE |
| project | TRUEでxの座標系へ投影する。FALSEのときはEPSG:3857（Webメルカトル）のまま | TRUE |
| verbose | TRUEでタイルの保存先、拡大の程度、出典を表示 | FALSE |
| apikey | APIキー。ThunderforestとStadiaは、環境変数THUNDERFOREST_MAPSまたはSTADIA_MAPSを設定していれば要らない | なし |
| cachedir | タイルを貯めておくフォルダの名前。指定しないと一時フォルダに貯める | なし |
| forceDownload | TRUEで貯めてあるタイルを取り寄せ直す | FALSE |
| retina | TRUEで高解像度のタイルがあればそれを取り寄せる。StadiaとCARTOが対応 | TRUE |
＜おすすめのRに関する書籍です＞
www.amazon.co.jp
地点を地図に重ねる：plot_tilesコマンド
取り寄せたタイルはplot_tilesコマンドでプロットします。そのあとはaddオプションを付けた通常のプロットで、点や文字を重ねられます。タイルと重ねるデータは座標系をそろえてください。丸の大きさは面積が販売量に比例するよう、平方根を使っています。OpenStreetMapのタイルを載せるときは出典の表記が要るので、get_creditコマンドで文字列を受け取り、図の右下に記述しています。
# 背景地図をプロット：plot_tilesコマンド
plot_tiles(UjiHaikei)
# 直売所の位置を重ねる：addオプション
# 丸の面積が販売量に比例するよう平方根を使う：cexオプション
plot(st_geometry(ChokubaijoHeimen),
pch = 21, bg = "#2F5D2A", col = "white",
cex = sqrt(ChokubaijoHeimen$Hanbai) / 3, lwd = 3, add = TRUE
)
# 店の名前を丸の上に置く：posオプション
text(st_coordinates(ChokubaijoHeimen),
labels = ChokubaijoHeimen$Mise, pos = 3, cex = 1.1
)
# タイルの出典を右下に表記：get_creditコマンド
mtext(get_credit("OpenStreetMap"), side = 1, line = -1.5, adj = 1, cex = 0.8)背景地図をプロット：plot_tilesコマンド
| オプション | 意味 | 初期値 |
|---|---|---|
| x | SpatRasterのオブジェクト | なし |
| adjust | TRUEで画像を拡大縮小せずにプロットする。画像が画面より小さければ余白を足し、大きければ画面に合わせる。緯度経度のままの画像には効かない | FALSE |
| add | TRUEで既にあるプロットに重ねる | FALSE |
| … | bgalpha、smoothなど、plotRGBコマンドへ渡すオプション | なし |
地図の種類を選ぶ：get_providersコマンド
取り寄せられるタイルは44種類あります。ただしStadia Maps、Thunderforest、CartoDBのタイルはAPIキーの登録が要り、キーなしで取り寄せると「API KEY REQUIRED」と書かれた画像が返ってきます。キーなしで使えるのはOpenStreetMapの系統、OpenTopoMap、Esriの系統です。データを重ねる背景には、色の薄いEsri.WorldGrayCanvasが向いています。
# タイルの種類の数を確認：get_providersコマンド
length(get_providers())
[1] 44
# 種類の名前を先頭から10件確認
head(names(get_providers()), 10)
[1] "OpenStreetMap" "OpenStreetMap.DE"
[3] "OpenStreetMap.France" "OpenStreetMap.HOT"
[5] "OpenTopoMap" "Stadia.AlidadeSmooth"
[7] "Stadia.AlidadeSmoothDark" "Stadia.OSMBright"
[9] "Stadia.Outdoors" "Stadia.StamenToner"
# APIキーの要らない4種類を選ぶ
Shurui <- c("OpenStreetMap", "OpenTopoMap", "Esri.WorldGrayCanvas", "Esri.WorldImagery")
# 名前を置く場所に使う四角形の角
Waku <- st_bbox(HyojiHanni)
# 2行2列に並べる
par(mfrow = c(2, 2))
# 種類ごとに取り寄せてプロットし、左上に種類の名前を置く
for (Tairu in Shurui) {
Haikei <- get_tiles(HyojiHanni, provider = Tairu, zoom = 15, crop = TRUE)
plot_tiles(Haikei)
legend(
x = Waku["xmin"], y = Waku["ymax"], legend = Tairu,
bg = "white", xjust = 0, yjust = 1, cex = 1.1
)
}
# 並べる設定を元に戻す
par(mfrow = c(1, 1))種類ごとに取り寄せてプロットし、左上に種類の名前を置く
2地点の道のりを求める：osrmRouteコマンド
osrmパッケージのosrmRouteコマンドは、2地点をつなぐ道のりと、その所要時間と距離を返します。直線距離ではなく、実際の道に沿った形です。問い合わせ先は既定で公開のデモサーバなので、続けて何百回も呼ぶような使い方は避け、必要ならOSRMのサーバを自分で立ててosrm.serverオプションに指定してください。
# 徒歩の道のりを求める設定にする："car"、"bike"、"foot"を指定
options(osrm.profile = "foot")
# JR宇治駅前店から平等院表参道店までの道のりを取得：osrmRouteコマンド
# 細かい形で受け取る：overviewオプション
Michinori <- osrmRoute(src = ChokubaijoSf[3, ], dst = ChokubaijoSf[2, ], overview = "full")
# 所要時間（分）と距離（km）を確認
st_drop_geometry(Michinori[, c("duration", "distance")])
duration distance
3_2 11.49167 0.8613
# 背景地図と同じ座標系へ変換
MichinoriHeimen <- st_transform(Michinori, 6674)
# 背景地図に道のりを重ねてプロット
plot_tiles(UjiHaikei)
plot(st_geometry(MichinoriHeimen), col = "#C8102E", lwd = 5, add = TRUE)
# 出発地と目的地に丸を置く
plot(st_geometry(ChokubaijoHeimen[c(3, 2), ]),
pch = 21, bg = "white", col = "#C8102E", cex = 3, lwd = 4, add = TRUE
)
# タイルの出典を右下に表記
mtext(get_credit("OpenStreetMap"), side = 1, line = -1.5, adj = 1, cex = 0.8)背景地図に道のりを重ねてプロット
| オプション | 意味 | 初期値 |
|---|---|---|
| src | 出発地。経度と緯度のベクトル、データフレーム、行列、POINTのsfc、POINTのsfを指定する。複数あるときは先頭だけを使い、行の名前を識別に使う | なし |
| dst | 目的地。指定できる形はsrcと同じ | なし |
| loc | 出発地、経由地、目的地をまとめて指定する。先頭が出発地、続きが経由地、最後が目的地になる | なし |
| overview | 返す道のりの形。”full”で細かい形、”simplified”で簡略化した形、FALSEで時間と距離だけ | “simplified” |
| exclude | OSRMへ渡す除外の指定 | なし |
| osrm.server | 道のりを計算するサーバのURL | getOption(“osrm.server”) |
| osrm.profile | 道のりの種類。”car”、”bike”、”foot”など | getOption(“osrm.profile”) |
地図については、次の記事もあります。タイルの種類ごとの見え方はRで解析：背景地図を手軽に取得！44種類の地図タイルを自在に取得「maptiles」パッケージの紹介、市区町村ごとに色を塗り分けた地図はRで解析：日本地図もOK！インタラクティブなテーママップ作成「mapsf.gui」パッケージの紹介、道路や建物のデータそのものを取り寄せたいときはRで解析：OpenStreetMapのデータをsfオブジェクトとして取得できる「osmdata」パッケージの紹介、マウスで動かせる地図はRで解析：「leaflet」パッケージなどで使える日本の地図ファイルの作成例を確認してください。
この記事が誰かの役に立ちますように。
