Hatena::ブログ(Diary)

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

2017-12-24

Rでおっぱい山を分析する

13:26 | Rでおっぱい山を分析するを含むブックマーク

これはFOSS4G Advent Calendar 2017 24日目の記事です。


「おっぱい山」とは、おっぱいの形をした山(英:breast-shaped hill)のことです。(出典: Wikipedia)

今回は統計ソフト「R」を使って各地のおっぱい山を分析してみたいと思います。


なぜ、そんなことするかだって?

そこにおっぱい山があるからさ

(ジョージ・マロリー)


f:id:tmizu23:20171224121249p:image

Rと地理院タイルで作成した北アルプスCS立体図


まずは、おっぱい山を探す所からです。地理院地図の検索機能を利用してみましょう。


検索窓に"おっぱい山"と入力してみます。

しかしながら、0件、見つからず。東大の協力も力及ばずと言った結果となりました...

f:id:tmizu23:20171222224906p:image



仕方ないのでググってデータを集めました。潜在的なおっぱい山はもっとあるはずですが、なかなか表に出てこないようです。

データ:https://gist.github.com/tmizu23/29c24df1662bbf976e2407f843f43c15


areanamelatlon
北海道上士幌町西クマネシリダケ、ピリベツ岳43.524422143.229961
岩手県八幡平市七時雨山40.070009141.107197
福島県郡山市安達太良山37.621130140.287864
福島県天栄村二岐山37.247941139.965477
東京都小笠原村母島乳房山26.659560142.161412
東京都小笠原村父島乳頭山27.096348142.210046
埼玉県小川町笠山36.018509139.189873
山梨県甲府市要害山35.703142138.599333
広島県安芸太田町574ピーク34.583811132.211576
沖縄県宮古島市69.1ピーク24.750397125.405974

おっぱい山の可視化

データ分析の基本は可視化です。

これらのおっぱい山をRと地理院タイルを利用して、可視化してみましょう。


地理院タイルをRで読み込めるようにベクトルタイル Advent Calendar 2017 16日目の記事で作成したコードを改造しました。指定した範囲の画像タイル、PNG標高タイルをダウンロードし、Rasterクラスに格納するので、表示・分析に利用できます。


コード:https://github.com/tmizu23/gsitiles.R


 使い方


最初に、以下のRライブラリをインストールしてください。


install.packages(c("raster", "curl", "png", "jpeg", "spatialEco"),dependencies=T)

関数をsource("gsitiles.R")で読み込んでから、こんな感じで使用します。

#使用例
source("gsitiles.R")

#標高データを取得・表示・投影変換・書き出し
dem<-getTile(137.665215,36.294582,137.7,36.314582,15,type="dem5a")
plot(dem)
demUTM53<-projectRaster(dem, crs=CRS("+init=epsg:3099"))
writeRaster(demUTM53,"demUTM53.tif")

#標準地図とCS立体図を乗算合成・表示・geoTIFF書き出し
cs<-makeCS(dem,shaded=TRUE)
stdmap<-getTile(137.665215,36.294582,137.7,36.314582,15,type="std")
stdmap_cs<-overlayColor(stdmap,cs,type="multiply")
plotRGB(stdmap_cs)
writeRaster(stdmap_cs,"stdmap_cs.tif",datatype="INT1U") #datatypeを指定

#陰影起伏図の作成・表示
slope<-terrain(dem,"slope")
aspect<-terrain(dem,"aspect")
shade<-hillShade(slope,aspect)
plot(shade,col=grey(0:100/100),legend=F)

#航空写真と陰影起伏を透過表示
photo<-getTile(137.665215,36.294582,137.7,36.314582,15,type="photo")
plot(shade,col=grey(0:100/100),legend=F)
plotRGB(photo,alpha=180,add=T)

  • getTile関数で取得したい範囲の左下と右上の座標、ズームレベル、タイルの種類を指定します。(左下と右上の座標が同じ場合は、その座標を含む1タイルを取得します。)
  • タイルの種類は、標準地図:"std"、シームレスフォト:"photo"、標高PNG:"dem10b","dem5a","dem5b"を選択できます。提供されているズームレベルと範囲は、地理院タイルのページを参照ください。
  • crop=TRUEを指定すると、指定範囲で切り抜きます。デフォルト(FALSE)では、指定範囲を含むタイルの範囲になります。

この他に以下の関数を用意しました。

  • makeCS(dem,shaded=TRUE):

   標高データからCS立体図を作成。shaded=TRUEで陰影起伏図とオーバレイ合成したCS立体図を作成

  • overlayColor(stackA,stackB,type="overlay | multiply"):

   RGB値の入ったstackクラス同士を画像合成。種類は、オーバレイ合成:"overlay"、乗算合成:"multiply"。(stackAとstackBは同じ範囲、ズームレベルの必要があります)


CS立体図については、こちらの記事を参考にさせて頂きました。

http://waigani.hatenablog.jp/entry/2017/12/03/225834

https://qiita.com/frogcat/items/7e91d3070a7a8d3e2c94


おっぱい山の可視化

前置きが長くなりましたが、可視化の結果はこちらです。

f:id:tmizu23:20171224090542p:image

日本を代表するおっぱい山です。両ピークに名前が付いているのが特徴です。

f:id:tmizu23:20171224090538p:image

なかなかの良型です。おっぱい山として文句はありません。

f:id:tmizu23:20171224090533p:image

日本百名山からのエントリーです。片乳ながら、乳頭のリアリティを追求しています。

f:id:tmizu23:20171224090528p:image

デコルテが激しく何かを主張しています。もっとリスペクトされても良いかもしれません。

f:id:tmizu23:20171224090523p:image

母島のおっぱい山です。あるべくしてある。必然の結果です。

f:id:tmizu23:20171224090518p:image

父島のおっぱい山です。父にも乳はあるんです。

f:id:tmizu23:20171224090513p:image

地元では名が知られているが垢抜けない。伸びしろに期待です。

f:id:tmizu23:20171224090508p:image

史跡と一体となった妖麗なおっぱい山です。誘われたら断れません。

f:id:tmizu23:20171224090504p:image

一見、変哲もなくつまらない。それでいいんです。

f:id:tmizu23:20171224090456p:image

言わんとすることは分かるんです。でも危険と裏腹な恐ろしさがあります。


まとめ

地理院タイルとRで、おっぱい山の可視化が簡単にできるようになりました。

しかしながら、我々が為すべきことは山積みです。


どのようにして、おっぱい山ができたのか?

なぜ、おっぱい山はそこにあるのか?

これから、おっぱい山はどうなってゆくのか?


おっぱい山の頂き目指し、これからも一歩一歩登っていきたいと思います。

それでは、良いクリスマスを!

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

2017-12-09

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

11:46 | バイナリベクトルタイルをQGISで表示するを含むブックマーク

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

hfuさんの記事でベクトルタイルの理解が深まってきたので、このへんで実際にQGISを使ってベクトルタイルを見てみましょう。

ベクトルタイルはWebGIS界隈で整備が進んできましたが、ここになってデスクトップGISへの逆輸入も始まっているようです。


ArcGIS Pro:

https://pro.arcgis.com/ja/pro-app/help/data/services/use-vector-tiled-layers.htm

QGIS:

https://github.com/geometalab/Vector-Tiles-Reader-QGIS-Plugin


ということで、今回はバイナリベクトルタイルをQGISで表示してみます。


f:id:tmizu23:20171201184735p:image:w600

バイナリベクトルタイルをQGISで表示


QGISでバイナリベクトルタイルを表示してみる

QGISでバイナリベクトルタイルを表示するには、Vector Tiles Reader Pluginを使用します。QGIS 2.18以上で動作します。


1. QGIS 2.18にVector Tiles Reader Pluginを入れます。

プラグインの管理の設定画面で、「実験的プラグインも表示する」にチェックを入れるとリストに表示されます。


2. メニューから ベクタ > Vector Tiles Layer > Add Vector Tiles Layer..と選びます。

デフォルトで、以下のOベクトルタイルの配信データがセットされています。

  • Mapzen.com
  • OpenMapTiles.com
  • OpenMapTiles.com(with custom key)

これはOSMのデータをベクトルタイルに変換したもので、Mapzenバージョン、OpenMapTilesバージョンが選択できます。

何が違うのかはhfuさんのこちらの記事が参考になります。

https://qiita.com/hfu/items/a761c5875d73c164bcd4

https://qiita.com/hfu/items/afbafacb70197711ee29


3. 「Connect」を押してレイヤ情報を取得します。

レイヤの一覧が表示され、表示オプションを設定できます。

表示オプションは、「Base Map defaults」、「Analysys Defaults」、「Inspection Defaults」、「Manual」がセットされています。

  • Base Map defaults:表示用。JSONのスタイルを適用。32タイルを読み込み
  • Analysys Defaults:分析用。タイルの境界バッファ部分を切り取り、すべてのタイルをマージ。16タイル読み込み
  • Inspection Defaults:確認用。1タイルのみ読み込み。
  • Manual:自分でオプションを設定

f:id:tmizu23:20171201184743p:image:w600


4. 「Add」ボタンを押して、表示します。

f:id:tmizu23:20171201184706p:image:w600



たったこれだけで、全国編集可能なOSMのデータが取得できてしまいます。ベクトルタイル便利ですねー

※ちょっと前に試した時は、デフォルトのスタイルをいい感じで付けてくれてたけど、今日、プラグインをアップデートしたら、いまいちなスタイルになってました。まあ、そのへんは開発途中ということで。


こうなってくるとベクトルタイルも自分で作って、QGISで表示させたくなります。

ということで、、、

shpファイルからバイナリベクトルタイルを作る

自分のshpファイルからバイナリベクトルファイルを作ってみましょう。


1. shpファイルを用意します。

今回は、植生図のデータを利用しました。

http://gis.biodic.go.jp/webgis/sc-023.html


2. shpファイルをgeojsonに変換します。

QGISでgeojson形式で書き出せばOKです。


3. tippecanoeを用意します。

Windowsでtippecanoeを利用するために、Bash on Ubuntu on Windowsを入れておきます。インストール方法はググってください。

tiooecanoeのソースコードをダウンロードします。

https://github.com/mapbox/tippecanoe


bashコンソールでgit clone すればダウンロードできます。

git clone https://github.com/mapbox/tippecanoe.git

tippecanoeをコンパイルします。コンパイルに必要なソフトをapt-getでインストールしてからmakeします。

cd tippecanoe
sudo apt-get install make g++ sqlite3 libsqlite3-dev zlib1g-dev 
make
sudo make install

6. tippecanoeでgeojsonをpbf(mvt)に変換します。

以下のコマンドでgeojsonをpbfファイルに変換します。ファイルはvegフォルダに1/1/1.pbfのように作られます。変換オプションで、ズームレベルや簡素化なども設定できます。

 tippecanoe -e veg veg.geojson

以下のコマンドだと、mbtiles形式になります。

 tippecanoe -o veg.mbtiles veg.geojson

tippecanoeで変換すると拡張子はpbfとなるようですが、vector-tile-specを読むとmvtとすべきみたいですね。

https://github.com/madefor/vector-tile-spec/blob/master/2.1/README.md



バイナリベクトルタイルをネットで配信してQGISで表示する

バイナリベクトルタイルができたのでQGISで読み込めるように配信してみましょう。


1. レイヤ情報のtile.jsonを作成します。

Vector Tiles Reader Pluginにレイヤの情報を教えるために、tilejson形式のtile.jsonを作成します。

https://github.com/mapbox/tilejson-spec/tree/master/2.2.0


tilejsonのspecはリンクの通りですが、プラグインでは最低限の情報があれば大丈夫そうです。

tilesのURLは、配信場所のURLに合わせます。boundsやcenterはtippecanoeで変換されたvegフォルダの中のmetadata.jsonに書かれています。


こんな感じです。

{
  "tilejson": "2.2.0",
  "tiles": ["https://map.ecoris.info/mvt-tiles/veg/{z}/{x}/{y}.pbf"],
  "name": "veg mvt",
  "attribution": "Biodiversity Center of Japan",
  "maxzoom": 14,
  "minzoom": 0,
  "bounds": [140.746596,37.836344,140.871589,37.919672],
  "center": [140.767822,37.848832,14],
  "vector_layers": [{
    "id": "veg",
    "description": "vegetation data (Biodiversity Center of Japan)",
    "maxzoom": 14,
    "minzoom": 0
  }]
}



2. データをサーバーにアップロードします。

tile.jsonをvegフォルダの中に入れて、ウェブサーバーにアップロードします。


3. プラグインにタイル情報を登録します。

QGISのメニューからベクタ > Vector Tiles Layer > Add Vector Tiles Layer..で「New」を選択して、Nameとtile.jsonのURLを登録します。今のところNameに日本語は入れられません。


f:id:tmizu23:20171201184740p:image:w450


4. ConnectしてAddして表示します。

これで自分のベクトルタイルが表示されます。ベクトルタイルなので、自分の好きなように着色したり編集できます。読み込みオプションで結合していない場合は、境界部分のバッファーに重なりがあります。


f:id:tmizu23:20171209114121p:image:w600



おまけ github上のベクトルタイルをQGISで読み込む

hfuさんがgithubでベクトルタイルを公開されているので、それをQGISで読み込んでみます。

試しに、このベクトルタイルを表示してみます。

https://github.com/hfu/chome-vt


自分のgistでレイヤ情報を作成します。

https://gist.github.com/tmizu23/5d427e657e6553c361881a3ec16683ed


RAWボタンを押して、そのURLをプラグインのタイル情報のTile Json URLに設定します。


もちろん、属性のラベルも表示できます。

f:id:tmizu23:20171209113540p:image:w600



おまけ mbtilesをQGISで読み込む

tippecanoeでmbtilesに変換した場合も、Vector-Tiles-Readerで表示できます。

QGISのメニューからベクタ > Vector Tiles Layer > Add Vector Tiles Layer..で「MBTiles」タブを選択して、ファイルを選べばOKです。


※pbfも「Directoris」タブで読み込めるはずですが、属性に日本語があるとmetadata.jsonの読み込みでエラーになるようです。(こっちはmetadata.jsonを参照するようです。)


まとめ

バイナリベクトルタイルもQGISでデータを取得、表示できるようになってきました。ベクトルタイルは、データそのものなので、編集も着色も自由にできます。WebGIS用と思われているベクトルタイルですが、その恩恵はデスクトップGISにもありそうです。バイナリベクトルタイルの普及が進んで、県ごとのshpファイルを別々にダウンロードして、解凍して、結合してみたいな面倒なことをしなくて済むようになると良いですね;-)



(そういえば、WFSってのが昔あった気がしますが、どうなったんですかねー?)


来週は、Vector-Tiles-Readerのソースコードを移植して、Rでバイナリベクトルタイルを読み込めるようにしてみたいと思います。


では、また。

2017-11-21 QGISでWMTSのズームレベルを制限して印刷する このエントリーを含むブックマーク

QGISのWMTSで地理院地図を表示させて、印刷しようとすると、ズームレベル17のデザインで印刷したいのに、ズームレベル18の詳細な地図が印刷されてしまうことがあります。

それを防ぐ方法を紹介します。

1. 地理院のGithubからWMTSの定義ファイルをダウンロードします。

https://github.com/gsi-cyberjapan/experimental_wmts

2. 標準地図の定義でz2to18となっているところを、z2to17と変更します。

<Layer>

<ows:Title>標準地図</ows:Title>

<ows:Identifier>std</ows:Identifier>

...

...

...

<TileMatrixSet>z2to17</TileMatrixSet>

...

...

...

</Layer>

3. z2to18のTileMatrixSetの定義の部分をコピーして、下に貼り付け、z2to18をz2to17に変更し、ズーム18のTileMatrixを削除します。

<TileMatrixSet>

<ows:Identifier>z2to18</ows:Identifier>

...

...

...

<TileMatrix>

<ows:Identifier>18</ows:Identifier>

...

...

</TileMatrix>

</TileMatrixSet>

<TileMatrixSet>

<ows:Identifier>z2to17</ows:Identifier>

...

...

...

<TileMatrix>

<ows:Identifier>17</ows:Identifier>

...

...

</TileMatrix>

</TileMatrixSet>

4. 定義ファイルを自分のwebサーバーに置いて、QGISのWMTS機能で呼び出します。

これで、ズームレベル17までの地図しか表示しないので、思った通りの印刷ができます。

補足

TileLayer Pluginの定義の指定でも同様のことができるのかもしれませんが、試していません。

2017-03-28 写真とGPSをリンクさせてQGISで写真を表示する方法 このエントリーを含むブックマーク

写真の撮影時間とGPSログをリンクさせてQGISで撮影位置と写真を表示する方法を紹介します。

前提として、GarminのGPSでトラックログを取りながら、GPS非対応のデジカメで写真を撮っていることを想定しています。デジカメの時間を、GPSの時間とズレないように正しい時間にセットしておきます。

1. GPSログを利用して写真に位置情報を埋め込む

ExifToolをダウンロードします。http://www.sno.phy.queensu.ca/~phil/exiftool/

GPSのトラックログをGPXファイルで保存しておきます。

写真をpicフォルダに入れておきます。

以下のコマンドをコマンドプロンプトで実行すると、picフォルダの中の写真全部に位置情報が埋め込まれます。

exiftool(-k).exe -geotag=Current.gpx pic

カメラの時間がズレている場合は、以下のコマンドで調整できます。(GPSの時間から1分20秒引いて同期させる場合)

exiftool(-k).exe -geosync=-1:20 -geotag=Current.gpx pic

2. QGISにPhoto2shapeプラグインをインストール

「森林土木memo」さんの情報をもとにPhoto2shapeプラグインをインストールします。

http://koutochas.seesaa.net/article/417307012.html

3. Photo2shapeの実行

picフォルダと出力するシェープファイルを指定、"Append to existing file"のチェックを外して実行します。

写真に埋め込まれた位置情報をもとにポイントが作成されます。

4. 写真を表示

追加されたレイヤのプロパティ→ディスプレイを開いてHTMLに以下を追加します。属性テーブルに写真のパスと日付が入っているので、それを表示させます。

<img src="[% filepath %]" width="400" height="300"><br/>
[% filename %]<br/>
[% img_date %]

picフォルダをqgsプロジェクトファイルとともに移動させたい場合は以下のように指定します。

<img src="[% @project_folder %]\pic\[% filname %]" width="400" height="300"><br/>
[% filename %]<br/>
[% img_date %]

ツールバーのマップチップスのアイコンを選択して、ポイントにマウスを合わせると写真が表示できます。

5. 写真を表示2

eVisプラグインでも写真を表示できます。

eVisプラグインはデフォルトでは有効になっていないので、「プラグインの管理」でチェックを入れます。

Photo2Shapeで読み込んだレイヤを選択して、メニューのデータベースからeVisイベントブラウザを起動します。

オプション→ファイルパスで、filepathを選択して「記憶する」にチェックを入れて保存します。

メニューのデータベースからeVisイベントIdツールを選択して、Photo2Shapeで読み込んだポイントをクリックすると写真が表示できます。

オプションの設定で、属性のfilenameを利用すれば、相対パスにも設定できます。

参考:

http://www.sno.phy.queensu.ca/~phil/exiftool/geotag.html

https://staff.aist.go.jp/t-yoshikawa/Geomap/QGIS_memo.html

http://koutochas.seesaa.net/article/417307012.html

Connection: close