Hatena::ブログ(Diary)

自然環境保全のための周辺技術

2017-12-16

バイナリベクトルタイルをRで表示する

08:07 | バイナリベクトルタイルをRで表示するを含むブックマーク

これはベクトルタイル Advent Calendar 2017 16日目の記事です。


前回、QGISのVector Tiles Reader Pluginでバイナリベクトルタイルを表示してみましたので、今回は、そのコードを参考にRでバイナリベクトルタイルを表示できるようにしてみます。


できあがりイメージ


OpenStreetMapデータの場合

hfuさんのOSM日本ベクトルタイルデータを使わせてもらいました。)


f:id:tmizu23:20171214151819p:image:w600



植生データの場合


f:id:tmizu23:20171215091302p:image:w600




コード

ここにあります。

https://github.com/tmizu23/mvt.R



使い方

とりあえず使い方です。


コード全部を取ってきて、mvt.Rを上から順に実行します。前半部分は、バイナリベクトルタイルを読みこむための関数読み込みで、example以降が、メインの部分です。実行前に、以下のRライブラリのインストールも必要です。なお、今のところバイナリデータのデコードに、C++で書かれたdllを利用する関係で、動作するのはWindowsのみです。


install.packages(c("Rcpp", "sf", "curl","jsonlite", "lwgeom", "dplyr", "purrr"),dependencies=T)

コードの解説

以下、コードの解説です。まずは、メインの実行部分です。

  • "getMVT"がバイナリベクトルタイルを取得するための関数です。
  • 引数は、ベクトルタイルのURL、拡張子、取得したい場所の左下と右上の座標(緯度経度)、ズームレベル、になります。
  • 取ってきたデータは、ベクトルタイルに格納されているレイヤのリストになっています。
  • ただし、同じレイヤ名でも、ジオメトリのタイプが違う場合があるので、レイヤ名とジオメトリタイプの組み合わせになっています。(ラインの道路、ポリゴンの道路とか)
  • リストからレイヤを取り出すと、sfクラスになっているので、簡単にプロットやデータ処理ができます。
  • 複数タイルにまたがる場合は、単純にコンバインしてあります。
  • なので、タイル境界のデータは、オーバーラップしてて、データ処理に使いたい場合は、ディゾルブする必要があります。

#------------------------------------------------------
# example 
#------------------------------------------------------

#植生ベクトルデータ
baseurl<-"http://map.ecoris.info/mvt-tiles/veg/"
ext<-".pbf"
veg<-getMVT(baseurl,ext,140.74963,37.85554,140.75963,37.85554,14)
veg$vegPolygon["HANREI_N"] %>% plot
merged_veg<-veg$vegPolygon %>% group_by(HANREI_N) %>% summarise_all(first)

#OSMベクトルデータ tiled by hfu
baseurl<-"https://hfu.github.io/jp1710_"
ext<-".mvt"
osm<-getMVT(baseurl,ext,139.767052,35.65054,139.797052,35.68054,13)
osm$transportationLineString["class"] %>% plot(key.size = lcm(4),main="Rでバイナリベクトルタイル")
osm$buildingPolygon["id"] %>% plot(add=T,col="black",border=0)

ここからが、実装部分の説明です。とりあえず、必要なライブラリを読み込みます。


library(Rcpp)
library(sf)
library(curl)
library(jsonlite)
library(lwgeom)
library(dplyr)
library(purrr)

次に、cppファイルを読み込みます。初めて読み込む時は、コンパイルされるので少し時間がかかります。

mvtはバイナリデータなので、そのデコードには、C++で書かれたdllを利用しています。その部分は、Vector-Tiles-Reader-QGIS-Pluginのpbf2geojsonを利用させてもらいました。pbf2geojsonはdllになっているので、それをRから呼び出せるようにRcppを利用してラッパーを書いて、それを呼び出しています。今のところwindows用のdllを決め打ちで呼び出しているので、windowsでしか利用できないのですが、この部分を変更すれば、linuxとmacからも呼び出せるようになると思います。ただ、Rcppの仕組みを良く理解していないので、未対応です。

https://github.com/geometalab/Vector-Tiles-Reader-QGIS-Plugin/tree/master/ext-libs/pbf2geojson


# cppファイルのコンパイルと読み込み
sourceCpp("pbf2geojsonWrapper.cpp")

タイルのx,y,zoomからタイルのWEBメルカトル座標と解像度を返す関数です。


tileinfo<-function(x,y,zoom){
  #タイル座標からタイルの左下、右上のWEBメルカトル座標とm/pixcelを返す
  A<-20037508.342789244
  L<- 2*A/2^zoom
  x0<-L*x-A
  y1<-A-L*y
  x1<-x0+L
  y0<-y1-L
  (list(bbox=c(x0,y0,x1,y1),res=L/256))
}

緯度、経度、ズームレベルからタイル座標を返す関数です。この関数とRでタイル画像を取得するアイデアは、uriboさんのブログから頂きました。uriboさんには今年のFOSS4G Hokkaidoのハンズオンでお世話になり、そこでsfとdplyrの存在を教えてもらいました!


latlon2tile <- function(lon, lat, zoom) {
  #緯度、経度、ズームレベルからタイル座標を返す
  x = trunc((lon/180 + 1) * 2^zoom/2)
  y = trunc(((-log(tan((45 + lat/2) * pi/180)) + 
                pi) * 2^zoom/(2 * pi)))
  return(c(x,y,zoom))
}

タイル一つ分のベクトルタイルをデコードする部分です。

処理手順は、ベクトルタイルをダウンロード > バイナリデータをデコード > データ整形してgeojsonに変換 > sfクラスに変換 です。

バイナリデータのデコードは、pbf2geojsonWrapper中のdecodePBFを呼び出し、実際の処理はpbf2geojsonでしてもらってます。


loadMVT<-function(baseurl,ext,x,y,zoom){
  #レイヤ名&ジオメトリタイプごとのgeojsonをリストで返す
 
  myurl<-paste0(baseurl, paste(zoom, x, y, sep = "/"), ext)
  info<-tileinfo(x,y,zoom)
  tileX<-info$bbox[1] #タイル左上のWEBメルカトルX座標
  tileY<-info$bbox[4] #タイル左上のWEBメルカトルY座標
  tileSpanX<-info$res*256 #タイルのX方向長さ
  tileSpanY<-info$res*256*-1 #タイルのY方向長さ(負にする)
  #バイナリデータ読み込み
  print(myurl)
  con <- gzcon(curl(myurl))
  mvt <- readBin(con, raw(), n=1e6, endian = "little")
  close(con)
  mvt<-paste(mvt,collapse = "")
  #バイナリデータをjsonに変換(C++のライブラリを呼び出し)
  mvtjson<-decodePBF(zoom,x,y,tileX,tileY,tileSpanX,tileSpanY,mvt)
  Encoding(mvtjson)<-"UTF-8"
  mvtjson<-gsub("_row","__row",mvtjson)
  #jsonをリストに変換して、レイヤ名&ジオメトリタイプごとにgeojsonを作成
  mvtlist<-fromJSON(mvtjson)
  layer_names<-names(mvtlist)
  geo_types = c("Point","MultiPoint","Polygon","MultiPolygon","LineString","MultiLineString")
  layer_list<-list()
  for(layer_name in layer_names){
    layer<-mvtlist[[layer_name]]
    for(geo_type in geo_types){
      features = layer[[geo_type]]
      if(!is.null(features)&&length(features)!=0){
        empty_geojson<-paste('{"source": "", "tiles": [], "features": [], "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:EPSG::3857"}}, "type": "FeatureCollection", "scheme": "xyz", "zoom_level":', zoom, ',"layer":"',layer_name,'"}',sep="")
        geojsonlist<-fromJSON(empty_geojson)
        features$id<-NULL
        features$properties$`id`<-seq(1:nrow(features))
        features$properties$`_col`<-NULL
        features$properties$`__row`<-NULL
        features$properties$`_zoom`<-NULL
        geojsonlist$features<-features
        geojson<-toJSON(geojsonlist,auto_unbox = TRUE)
        geo<-read_sf(geojson)
        geo<-st_make_valid(geo)
        layer_type<-paste(layer_name,geo_type,sep="")
        layer_list[[layer_type]]<-geo
      }
    }
  }
  return(layer_list)
}

緯度、経度、ズームレベルから、ダウンロードするタイルを計算して、処理を呼び出す部分です。


getMVT<-function(baseurl,ext,lon0,lat0,lon1,lat1,zoom,combine=TRUE){
  #緯度経度の左下〜右上の範囲のベクトルタイルを取得

  xyz0 <- latlon2tile(lon0,lat0,zoom)
  xyz1 <- latlon2tile(lon1,lat1,zoom)
  features<-list()
  for(y in xyz0[2]:xyz1[2]){
    for(x in xyz0[1]:xyz1[1]){
      layer_list<-loadMVT(baseurl,ext,x,y,zoom)
      tileid<-paste(x,y,sep="-")
      features[[tileid]]<-layer_list
    }
  }
  #tileidとlayer_typeの階層を入れ替える  
  tmp_fet<-unlist(unname(features), recursive=FALSE)
  features <- split(setNames(tmp_fet, rep(names(features), lengths(features))), names(tmp_fet))

    #タイルを1つにする
  if(combine==TRUE){
    features<-features %>% map(bind_rows.sf)
  }
  return(features)
}

sfクラスに変換してリスト形式に格納している各タイルを一つのsfクラスにバインドする部分です。


bind_rows.sf<-function(x){
  #sfクラスのリストをバインドする。列にミスマッチがある場合はNAにする。
  a<-x %>% map(~dplyr::select(.,geometry)) %>% reduce(rbind)
  b<-x %>% map(~st_set_geometry(.,NULL)) %>% reduce(bind_rows)
  geom<-st_sf(bind_cols(a,b))
  return(geom)
}

だいたい、こんな感じです。


やり残し

  • Windows以外への対応(dll呼び出し部分をOSごとに変更すれば出来そうですが)
  • 自動ディゾルブ処理(indexキーがmvtの仕様にないので、自動でディゾルブ処理するには位置関係から判定する必要があってめんどいです)
  • mapbox glとかleafletとかopenlayersのバイナリベクトルデータのデコード&ディゾルブ処理を参考にする(pbf2geojsonよりも速くて上手いこと処理してそうな気がする)
  • style.jsonへの対応
  • パッケージ化(気が向いたら、誰か改良してパッケージ化してくれると嬉しいです)

ホントは、全部C++で書いて、Rcppで呼び出すだけの方が速いんだろーなという気もしますが、それは言わないお約束;-)


まとめ

RでGISをするときに、色々なデータをネット経由で簡単に取得できると良いなーと思います。ベクトルタイルの普及が進めば、それが可能になるかもしれません。あと、sfパッケージとdplyrで、GISデータ処理が楽にできるようになりましたが、どのコマンドを組み合わせれば望む処理ができるのか、発見するまでに時間がかかってしまいます。なので、ベクトルタイルのように簡単に取ってこれるデータを使った、R&GIS逆引き辞典を作りたい作って欲しいなーと思います。


こんな資料も作っておりましたので、参考まで

https://docs.google.com/presentation/d/1JMaeGPFVTEqNfmIZ1puMPOJM6zTPy9vcytCl2qf9U3k/edit

Connection: close