本文へスキップ
からだにいいもの

Rのトピックスを中心に『まだ、まだ、知らない、役に立つ情報?』を発信します。

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_containsNominatimが返す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”
featuretypeOSMのfeature種別(settlementが既定)“settlement”
limitAPIが返す結果の最大件数10
keyAPIキーが必要なサービスで使用するキー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つのクエリに複数のタグ条件をまとめて追加します。キーと値の組をリストやベクトルで渡すことで、異なる種類の地物を一度に検索できます。

オプション意味初期値
opqoverpass_queryオブジェクトなし
featureslist(“<key>” = “<value>”) または c(“<key>” = “<value>”) 形式の名前付きリスト・ベクトル、もしくはキーと値をエスケープした引用符で囲んだ文字列ベクトルなし
bbox地物検索の対象範囲。opq側に範囲が設定されていない場合は指定が必須NULL
key_exactFALSEにするとキーを完全一致で解釈しないTRUE
value_exactFALSEにすると値を完全一致で解釈しない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が行名として付きます。この関数は、その行名と列名を取り除きます。

オプション意味初期値
xosmdata_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

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

スポンサーリンク
価格および配送状況は変更される場合があります。購入時は商品ページをご確認ください。
当サイトに表示されている商品情報はAmazonから提供されたものであり、更新または削除される場合があります。
karada-goodはAmazonアソシエイトとして、適格販売により収入を得ています。