コードの解説
以下、コードの解説です。まずは、メインの実行部分です。
- "getMVT"がバイナリベクトルタイルを取得するための関数です。
- 引数は、ベクトルタイルのURL、拡張子、取得したい場所の左下と右上の座標（緯度経度）、ズームレベル、になります。
- 取ってきたデータは、ベクトルタイルに格納されているレイヤのリストになっています。
- ただし、同じレイヤ名でも、ジオメトリのタイプが違う場合があるので、レイヤ名とジオメトリタイプの組み合わせになっています。（ラインの道路、ポリゴンの道路とか）
- リストからレイヤを取り出すと、sfクラスになっているので、簡単にプロットやデータ処理ができます。
- 複数タイルにまたがる場合は、単純にコンバインしてあります。
- なので、タイル境界のデータは、オーバーラップしてて、データ処理に使いたい場合は、ディゾルブする必要があります。
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)
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
sourceCpp("pbf2geojsonWrapper.cpp")
タイルのx,y,zoomからタイルのWEBメルカトル座標と解像度を返す関数です。
tileinfo<-function(x,y,zoom){
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){
myurl<-paste0(baseurl, paste(zoom, x, y, sep = "/"), ext)
info<-tileinfo(x,y,zoom)
tileX<-info$bbox[1]
tileY<-info$bbox[4]
tileSpanX<-info$res*256
tileSpanY<-info$res*256*-1
print(myurl)
con <- gzcon(curl(myurl))
mvt <- readBin(con, raw(), n=1e6, endian = "little")
close(con)
mvt<-paste(mvt,collapse = "")
mvtjson<-decodePBF(zoom,x,y,tileX,tileY,tileSpanX,tileSpanY,mvt)
Encoding(mvtjson)<-"UTF-8"
mvtjson<-gsub("_row","__row",mvtjson)
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
}
}
tmp_fet<-unlist(unname(features), recursive=FALSE)
features <- split(setNames(tmp_fet, rep(names(features), lengths(features))), names(tmp_fet))
if(combine==TRUE){
features<-features %>% map(bind_rows.sf)
}
return(features)
}
sfクラスに変換してリスト形式に格納している各タイルを一つのsfクラスにバインドする部分です。
bind_rows.sf<-function(x){
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)
}
だいたい、こんな感じです。
