Rで解析:OpenStreetMapのデータをsfオブジェクトとして取得できる「osmdata」パッケージの紹介
地図データを使った空間解析では、対象とする地域の道路・建物・施設などの情報を集める作業が起点になります。しかし、OpenStreetMapから必要な範囲と種類のデータだけを取り出し、解析に使える形式へ整えるには手間がかかります。「osmdata」パッケージは、Overpass APIへの問い合わせを組み立て、OpenStreetMapのデータをsfオブジェクトやspオブジェクトとして取得します。地名からの境界ボックスの取得、複数のタグ条件を組み合わせた検索、地点の周囲や特定の領域内に限定した抽出が可能です。また、取得したデータの行名の整理や、送信前のクエリ文字列の確認も可能です。本パッケージの利用で、地図データの取得から解析までの流れを効率化できるのではないかと考えます。
パッケージバージョンは0.4.1。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。
パッケージのインストール
下記コマンドを実行してください。
# パッケージのインストール
install.packages("osmdata")
# パッケージの読み込み
library("osmdata")
# コマンド例で使用するパッケージの読み込み
# install.packages("sf")
library("sf")コマンド例
詳細はコメント、パッケージのヘルプを確認してください。
OpenStreetMap(オープンストリートマップ)は、世界中の地図情報を誰でも編集・利用できる形で蓄積している共同プロジェクトです。 Overpass API(オーバーパスAPI)は、そのデータベースから条件に合う地物だけを取り出すための問い合わせサービスです。
osmdataでは、まず取得したい範囲(境界ボックス)とタグの条件を指定してクエリを組み立て、Overpassサーバーへ送信してデータを受け取ります。受け取ったデータは、sfパッケージのsfオブジェクトなど、Rで扱いやすい形式へ変換されます。以下では、範囲や地点を指定してクエリを組み立てる関数と、取得後のデータを整える関数を取り上げます。
# 解析対象とする京都府内の地点(説明用の擬似データ)
uji_bashi <- c(lon = 135.8074, lat = 34.8917)
# 亀岡駅前の位置
kameoka_eki <- c(lon = 135.5776, lat = 35.0117)
# 嵐山周辺の検索範囲(左下と右上の経度・緯度)
arashiyama_bb <- c(135.665, 35.005, 135.685, 35.020)地名から境界ボックスを取得:getbbコマンド
地名から、その範囲を囲む境界ボックスの座標を取得します。取得先はNominatimという地名検索サービスで、出力形式は行列や文字列、データフレームなどから選べます。
| オプション | 意味 | 初期値 |
|---|---|---|
| place_name | 検索対象の地名またはWikidata ID | なし |
| display_name_contains | Nominatimが返すdisplay_nameと照合する文字列 | NULL |
| viewbox | 検索を絞り込む範囲。”x1,y1,x2,y2″ の文字列またはc(x1, y1, x2, y2)で指定する | NULL |
| format_out | 出力形式(matrix、data.frame、string、polygon、sf_polygon、osm_type_id) | c(“matrix”, “data.frame”, “string”, “polygon”, “sf_polygon”, “osm_type_id”) |
| base_url | データ取得先のベースURL | “https://nominatim.openstreetmap.org” |
| featuretype | OSMのfeature種別(settlementが既定) | “settlement” |
| limit | APIが返す結果の最大件数 | 10 |
| key | APIキーが必要なサービスで使用するキー | NULL |
| silent | 問い合わせURLを画面に表示するかどうか | TRUE |
# 「宇治市」の境界ボックスを行列形式で取得
getbb("宇治市 京都府")
min max
x 135.75965 135.87992
y 34.85792 34.95737
# Overpassクエリにそのまま渡せる文字列形式で取得
getbb("亀岡市 京都府", format_out = "string")
[1] "34.915167,135.37765,35.10129,135.646745"
# featuretypeを既定のsettlementのまま「南丹市」を取得
getbb("南丹市 京都府")
min max
x 135.36873 135.79516
y 35.02455 35.37787
複数の地物条件をまとめて指定:add_osm_featuresコマンド
1つのクエリに複数のタグ条件をまとめて追加します。キーと値の組をリストやベクトルで渡すことで、異なる種類の地物を一度に検索できます。
| オプション | 意味 | 初期値 |
|---|---|---|
| opq | overpass_queryオブジェクト | なし |
| features | list(“<key>” = “<value>”) または c(“<key>” = “<value>”) 形式の名前付きリスト・ベクトル、もしくはキーと値をエスケープした引用符で囲んだ文字列ベクトル | なし |
| bbox | 地物検索の対象範囲。opq側に範囲が設定されていない場合は指定が必須 | NULL |
| key_exact | FALSEにするとキーを完全一致で解釈しない | TRUE |
| value_exact | FALSEにすると値を完全一致で解釈しない | TRUE |
# 嵐山周辺の検索範囲から、飲食店と喫茶店の2条件を持つクエリを作成
q_food <- opq(bbox = arashiyama_bb) |>
add_osm_features(features = list(
"amenity" = "restaurant",
"amenity" = "cafe"))
# 組み立てたクエリ文字列を確認
cat(opq_string(q_food))
[out:xml][timeout:25];
(
node ["amenity"="restaurant"] (35.005,135.665,35.02,135.685);
way ["amenity"="restaurant"] (35.005,135.665,35.02,135.685);
relation ["amenity"="restaurant"] (35.005,135.665,35.02,135.685);
node ["amenity"="cafe"] (35.005,135.665,35.02,135.685);
way ["amenity"="cafe"] (35.005,135.665,35.02,135.685);
relation ["amenity"="cafe"] (35.005,135.665,35.02,135.685);
);
(._;>;);
out body;# キーと値を完全一致にせず、名称に「嵐山」を含む地物を対象にする
q_name <- opq(bbox = arashiyama_bb) |>
add_osm_features(features = list("name" = "嵐山"),
key_exact = FALSE, value_exact = FALSE)
# 完全一致でない指定がクエリにどう反映されるかを確認
cat(opq_string(q_name))
[out:xml][timeout:25];
(
node [~"name"~"嵐山"] (35.005,135.665,35.02,135.685);
way [~"name"~"嵐山"] (35.005,135.665,35.02,135.685);
relation [~"name"~"嵐山"] (35.005,135.665,35.02,135.685);
);
(._;>;);
out body;地点の周囲を半径で検索:opq_aroundコマンド
指定した地点を中心とする円の範囲に対象を絞ってクエリを組み立てます。半径はメートル単位で指定します。
| オプション | 意味 | 初期値 |
|---|---|---|
| lon | 対象地点の経度 | なし |
| lat | 対象地点の緯度 | なし |
| radius | データを抽出する対象地点からの半径(メートル)。大きな値を指定すると問い合わせが失敗することがある | 15 |
| key | (任意)囲む対象データのOSMキー | NULL |
| value | (任意)keyに対応する囲む対象データのOSM値 | NULL |
| timeout | 大きな問い合わせで必要に応じて増やす秒数。サーバーが全データ送信前に応答を打ち切ることがある | 25 |
# 宇治橋の位置の周囲、半径200メートルからベンチを探すクエリを作成
q_bench <- opq_around(uji_bashi[["lon"]], uji_bashi[["lat"]],
radius = 200, key = "amenity", value = "bench")
# 戻り値はOverpassクエリの文字列
cat(q_bench)
[out:xml][timeout:25];(node(around:200,34.8917, 135.8074)[amenity=bench];way(around:200,34.8917, 135.8074)[amenity=bench];relation(around:200,34.8917, 135.8074)[amenity=bench];);
(._;>;);
out;
地点を囲む領域を取得:opq_enclosingコマンド
指定した地点を内部に含む領域(行政区域や建物など)を取得するクエリを組み立てます。返す型は関係(relation)か経路(way)から選べます。
| オプション | 意味 | 初期値 |
|---|---|---|
| lon | 対象地点の経度 | NULL |
| lat | 対象地点の緯度 | NULL |
| key | (任意)囲む対象データのOSMキー | NULL |
| value | (任意)keyに対応する囲む対象データのOSM値 | NULL |
| enclosing | 囲む要素として返す型。’relation’(多くはマルチポリゴン)または’way’(多くはポリゴン) | “relation” |
| timeout | 大きな問い合わせで必要に応じて増やす秒数。サーバーが全データ送信前に応答を打ち切ることがある | 25 |
# 亀岡駅前の位置を内部に含む行政区域を取得するクエリを作成
q_area <- opq_enclosing(kameoka_eki[["lon"]], kameoka_eki[["lat"]],
key = "boundary", value = "administrative")
# クエリ文字列を確認
cat(opq_string(q_area))
[out:xml][timeout:25];
(
is_in(35.0117,135.5776)->.a;relation(pivot.a)["boundary"="administrative"];);
(._;>;);
out;sfオブジェクトから行名を除去:unname_osmdata_sfコマンド
osmdata_sf()が返すオブジェクトでは、ジオメトリの座標行列にOSMのノードIDが行名として付きます。この関数は、その行名と列名を取り除きます。
| オプション | 意味 | 初期値 |
|---|---|---|
| x | osmdata_sf()関数が返したosmdata_sfオブジェクト | なし |
# 嵐山周辺の歩道をsfオブジェクトとして取得
arashiyama_paths <- opq(bbox = arashiyama_bb, timeout = 60) |>
add_osm_feature(key = "highway", value = "footway") |>
osmdata_sf()
# ジオメトリの座標行列にはOSMのノードIDが行名として残る
head(as.matrix(arashiyama_paths$osm_lines$geometry[[1]]), 3)
lon lat
4277759095 135.6747 35.01589
4264086521 135.6745 35.01590
4264086522 135.6745 35.01590
# 行名(ノードID)と列名を取り除く
arashiyama_clean <- unname_osmdata_sf(arashiyama_paths)
# 行名が外れ、位置だけの行列になる
head(as.matrix(arashiyama_clean$osm_lines$geometry[[1]]), 3)
[,1] [,2]
[1,] 135.6747 35.01589
[2,] 135.6745 35.01590
[3,] 135.6745 35.01590
この記事が誰かの役に立ちますように。