INPUTしたらOUTPUT!

忘れっぽいんでメモっとく

{osmdata}でエトセトラ(1/3) OpenStreetMapからデータ取得

R

こんな話がある。

qiita.com

dailyportalz.jp

これらをRでやりたい。{osmdata}パッケージを利用するとできたのでメモしとく。今回は道路ネットワークデータの取得編。



データ取得範囲の決定

前述のOSMnxで取得した道路ネットワークデータから極座標グラフを作ろう！ - QiitaではPythonのOSMnwモジュールでOpenStreetMapからデータを取得しているが、Rでは{osmdata}で取得できる。データを取得する範囲であるバウンディングボックスは以下の2通りの方法で指定することができる。

  • 都市名
  • 緯度・経度

前者はOSMデータの逆ジオコーディングツールであるNominatimでバウンディングボックスの座標を取得している。

江東区を例に取得してみる。

library(dplyr)

# Nominatimから江東区の行政境界を取得
city <-
  osmdata::getbb("江東区, 東京都, 日本", format_out = "sf_polygon")

# バウンディングボックスを取得
bbox <-
  city %>% 
  sf::st_bbox()
     xmin      ymin      xmax      ymax 
139.76805  35.55956 139.89250  35.70809


行政区域に海面が含まれておりバウンディングボックスも広い。{leaflet}で可視化してみると以下のようになる。

# 地図上に可視化
atr <- "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>"
m <-
  leaflet::leaflet(options = leaflet::leafletOptions(zoomControl = FALSE)) %>%
  leaflet::addTiles("http://cyberjapandata.gsi.go.jp/xyz/blank/{z}/{x}/{y}.png",
                    attribution = atr) %>%
  leaflet::addPolygons(data = city, weight = 1) %>%
  leaflet::addPolylines(
    data = sf::st_as_sfc(bbox),
    fillColor = "Transparent",
    color = "red"
  )
m

f:id:tak95:20200831225235p:plain
江東区バウンディボックス(Nominatim版)


国土数値情報 | 行政区域データでダウンロードした行政区域のデータからバウンディングボックスを取得するには次のようにする。

# 国土数値情報 行政区域データの読み込み
cities <-
  sf::read_sf("data/ksj/N03-200101_GML/N03-20_200101.shp",
              options = "ENCODING=CP932") %>% 
  # JGD2011からWGS84に変換
  sf::st_transform(4326)

# 江東区に絞り込み
city_code  <-  "13108"
city <-
  cities %>% 
  dplyr::filter(N03_007 == city_code)

# バウンディングボックスを取得
bbox <-
  city %>%  
  sf::st_bbox() 
bbox
     xmin      ymin      xmax      ymax 
139.77132  35.60184 139.84823  35.70801


同様に可視化してみると以下のようになる。（海の森公園が含まれていないが海面がない分こちらの方が良さそう）

f:id:tak95:20200830230853p:plain
江東区バウンディボックス(国土数値情報版)


道路ネットワークの取得

OpenStreetMapからデータを取得するにはまずosmdata::opq()で対象範囲のバウンディングボックスを指定する。次にosmdata::add_osm_feature()でタグを指定し、Overpass APIのリクエストを作成する。 以下を参考にkeyとvalueを指定し、道路のタグを設定する。

github.com

# 道路ネットワークを取得するリクエストクエリ
q <-
  osmdata::opq(bbox = sf::st_bbox(city)) %>%
  osmdata::add_osm_feature(
    key = 'highway',
    value = c(
      "trunk",           # 国道
      "trunk_link",
      "primary",         # 主要地方道
      "primary_link",
      "secondary",       # 一般都道府県道
      "secondary_link",
      "tertiary",        # 一般道(2車線以上)
      "unclassified",    # 一般道(2車線未満)
      "residential",     # 居住区域内道路
      "living_street",   # 歩行者優先道路
      "pedestrian",      # 歩行者専用道路
      "service",         # 敷地内道路
      "track"            # 砂利道
    )
  ) %>% 
  osmdata::add_osm_feature(key = "area", value = "!yes")
q

$bbox
[1] "35.60184045,139.77132484,35.70800764,139.84822678"

$prefix
[1] "[out:xml][timeout:25];\n(\n"

$suffix
[1] ");\n(._;>;);\nout body;"

$features
[1] " [\"highway\"~\"^(trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|unclassified|residential|living_street|pedestrian|service|track)$\"]"
[2] " [\"area\"!=\"yes\"]"                                                                                                                                          

attr(,"class")
[1] "list"           "overpass_query"


リクエストの結果は以下の3つの関数で取得できる。

  • osmdata:: osmdata_xml()
  • osmdata:: osmdata_sf()
  • osmdata:: osmdata_sp()

osmdata::osmdata_sf()によって返されるosmdataは他の空間オブジェクトで表されているかに関わらず全てのデータが含まれるためオブジェクトサイズが大きくなる（LINESTRING、POLYGONの頂点もPOINTに含まれる？）。osmdata:: osmdata_xml()で保存したxmlをsf::read_sf()でレイヤを指定して読み込んだ場合は固有のオブジェクトのみとなる。

lineオブジェクトでサイズを比較してみるとosmdataは3.4倍となった。メモリの節約には一度xmlに保存して読みこんだ方が良さそう。

# クエリの結果をxmlに保存
q %>% 
  osmdata::osmdata_xml(filename = glue::glue("data/osm/highway_{city_code}.osm"))

# クエリの結果をosmdataオブジェクトで取得
dat_sf <-
  q %>% 
  osmdata::osmdata_sf()

# オブジェクトサイズを比較
s1 <- object.size(dat_sf$osm_lines)
s2 <- object.size(sf::read_sf(glue::glue("data/osm/highway_{city_code}.osm"), 
                              layer = 'lines', quiet = TRUE))
as.numeric(s1 / s2)
[1] 3.566953


取得した道路ネットワークを可視化してみる。

# xmlからlinesレイヤーを読み込み
ways <-
  sf::read_sf(glue::glue("data/osm/highway_{city_code}.osm"), layer = 'lines')

# 地図で可視化
m %>% 
  leaflet::addPolylines(data = ways, color = "purple", weight = 2)

f:id:tak95:20200903001740p:plain
OSM道路ネットワーク

バウンディングボックス内の道路を取得しているため当然ながら江東区外の道路も含まれている。sf::st_intersection()で江東区の行政区域内の道路のみにする。

ways <-
  sf::read_sf(glue::glue("data/osm/highway_{city_code}.osm"), layer = 'lines') %>% 
  # 行政区域と交差する道路
  sf::st_intersection(city) %>% 
  # 国土数値情報 行政区域の列を除く
  dplyr::select(!dplyr::starts_with("N03_"))

f:id:tak95:20200903001818p:plain
OSM道路ネットワーク(江東区)


道路ネットワークデータの取得ができたので次は角度を求め、極座標グラフを作ってみる。