pos ファイルからLat Lon Altを計算するスクリプト - 自然環境保全のための周辺技術で作成したcsvファイルは緯度、経度なので、それを他の座標系に変換するスクリプトを作成しました。Photoscanへの入力時に平面直角座標系にしておきたい場合に使います。
pos ファイルからLat Lon Altを計算するスクリプト
QGIS Bezier Editing Plugin を作りました。
QGISでもIllustratorみたいにベジエ曲線で地物を描きたい！ということでBezier Editing というプラグインを作成し、公式プラグインに登録してみました。
プラグインでは、こんなことができます。
- アンカーとハンドルでベジエ曲線を描いてフィーチャーを作成。フィーチャーはアンカー間を10点で補間する仕様になっています。
- フィーチャーをベジエ曲線に逆変換して、再編集。ただし、プラグインで作成したフィーチャーだけが再編集可能です。
- ポイント、ライン、ポリゴンに対応。ポリゴンは、最初と最後のアンカーがスナップしている必要があります。
- フリーハンドの曲線もベジエ曲線に変換できます。
- 分割と結合も可能。
インストールは、QGISのメニューからプラグイン > プラグインの管理とインストール で 「Bezier Editing」 と検索してください。
操作のイメージは以下をどうぞ
編集ツール
フリーハンドツール
使い方のドキュメントは、こちらにあります。
github.com
バグやPull Requestがある場合は、こちらにお願いします。
github.com
QGIS3にプラグインを移植する方法
qgis2用に作ったプラグインをqgis3に対応させるためのメモです。
qgis2to3の準備
OSGeo4Wをインストール
（QGISのスタンドアローンインストーラーでインストールされるものだとpy3_envでパスがうまく指定されないのでpipがC:\OSGeo4Wのpythonを読み込もうとするので、OSGeo4Wネットワークインストーラーを使います。）
OSGeo4W shell を管理者で実行
py3_env
pip install qgis2to3
qgis2to3の実行
OSGeo4W shellを普通に実行
py3_env
cd C:\Users\mizutani\.qgis2\python\plugins
mkdir myplugin_qgis3
python "c:\OSGeo4W64\apps\Python37\Scripts\qgis2to3" -w -n -o myplugin_qgis3 myplugin
importの修正
importを以下に変更
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtWidgets import *
from qgis.PyQt.QtGui import *
from qgis.core import *
from qgis.gui import *
qgis3のプラグイン置き場に移動
mv myplugin_qgis3 %USERPROFILE%\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\myplugin
リソースのコンパイル
cd %USERPROFILE%\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\myplugin
qt5_env
pyrcc5 -o resources.py resources.qrc
QGIS3の起動とエラーの修正
qgis3を起動してプラグインでエラーになる箇所を修正
ドローンで撮影した動画と飛行軌跡を同期させOpenlayersで表示する方法
（とりあえず忘れないように殴り書き。あとで追記します。）
ドローン（Phantom4pro）で撮影した動画が、どこからどちらの方向を向いて撮影したのかOpenlayersで表示させる方法です。
動画と地図が同期しているので、地図上で場所を指定するとその場所を撮影した動画が表示され、逆に動画の再生位置を変更すると、地図上のアイコンが撮影場所に移動します。
フライトログと動画を同期する
フライトログの取得
DJI Assistant 2でDATファイルを取得します。
https://www.dji.com/jp/phantom-4-pro-v2/info#downloads
フライトログの変換
DatConでDATファイルをCSVファイルに変換します。
https://datfile.net/
変換されたCSVのヘッダーの意味はこちらで確認できます。
https://datfile.net/DatCon/fieldsV3.html
CSVの列の編集
CSVをEXCELで編集します。
以下の列を残して、不要な列を削除します。
| IMU_ATTI(0):Longitude | IMU_ATTI(0):Latitude | IMU_ATTI(0):yaw360 | flightTime |
列名を以下に変更します。
| lon | lat | yaw | sec |
sec列の単位はミリ秒となっているので、1000で割って小数以下を切り捨て単位を秒にしておきます。
CSVの行の編集（動画との同期）
フライトログと動画を同期させるためにCSVの行を削除します。
あらかじめ動画は必要な範囲を切り取るなどの編集をしておきます。
（切り貼りして不連続に編集した動画は同期できません）
CSVをQGISに読み込んで、フライトログを表示させます。
（x,yには、lat,lonを指定します）
動画の開始位置がCSVの開始行となるように、それ以前の不要なポイントを削除します。
同様に、動画の終了位置がCSVの最終行となるように、それ以降の不要なポイントを削除します。
QGISからCSVに書き出します。
CSVをEXCELで開き、sec列の先頭が0（秒）から始まるように、sec列全体からsec列の先頭の値を引きます。（最終行の秒数は、動画の終了秒とほぼ同じになるはず）
秒数が重複しているデータを削除します。
CSVをJSONに変換
CSV2JSONでCSVをJSONに変換します。Minifyにチェックしておきます。
https://www.csvjson.com/csv2json
Rでおっぱい山を分析する
これはFOSS4G Advent Calendar 2017 24日目の記事です。
「おっぱい山」とは、おっぱいの形をした山（英：breast-shaped hill）のことです。（出典： Wikipedia）
今回は統計ソフト「R」を使って各地のおっぱい山を分析してみたいと思います。
なぜ、そんなことするかだって？
そこにおっぱい山があるからさ
（ジョージ・マロリー）
まずは、おっぱい山を探す所からです。地理院地図の検索機能を利用してみましょう。
検索窓に"おっぱい山"と入力してみます。
しかしながら、0件、見つからず。東大の協力も力及ばずと言った結果となりました...
仕方ないのでググってデータを集めました。潜在的なおっぱい山はもっとあるはずですが、なかなか表に出てこないようです。
データ：https://gist.github.com/tmizu23/29c24df1662bbf976e2407f843f43c15
| area | name | lat | lon |
|---|---|---|---|
| 北海道上士幌町 | 西クマネシリダケ、ピリベツ岳 | 43.524422 | 143.229961 |
| 岩手県八幡平市 | 七時雨山 | 40.070009 | 141.107197 |
| 福島県郡山市 | 安達太良山 | 37.621130 | 140.287864 |
| 福島県天栄村 | 二岐山 | 37.247941 | 139.965477 |
| 東京都小笠原村母島 | 乳房山 | 26.659560 | 142.161412 |
| 東京都小笠原村父島 | 乳頭山 | 27.096348 | 142.210046 |
| 埼玉県小川町 | 笠山 | 36.018509 | 139.189873 |
| 山梨県甲府市 | 要害山 | 35.703142 | 138.599333 |
| 広島県安芸太田町 | 574ピーク | 34.583811 | 132.211576 |
| 沖縄県宮古島市 | 69.1ピーク | 24.750397 | 125.405974 |
おっぱい山の可視化
データ分析の基本は可視化です。
これらのおっぱい山をRと地理院タイルを利用して、可視化してみましょう。
地理院タイルをRで読み込めるようにベクトルタイル Advent Calendar 2017 16日目の記事で作成したコードを改造しました。指定した範囲の画像タイル、PNG標高タイルをダウンロードし、Rasterクラスに格納するので、表示・分析に利用できます。
使い方
最初に、以下の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
おっぱい山の可視化
前置きが長くなりましたが、可視化の結果はこちらです。
日本を代表するおっぱい山です。両ピークに名前が付いているのが特徴です。
なかなかの良型です。おっぱい山として文句はありません。
日本百名山からのエントリーです。片乳ながら、乳頭のリアリティを追求しています。
デコルテが激しく何かを主張しています。もっとリスペクトされても良いかもしれません。
母島のおっぱい山です。あるべくしてある。必然の結果です。
父島のおっぱい山です。父にも乳はあるんです。
地元では名が知られているが垢抜けない。伸びしろに期待です。
史跡と一体となった妖麗なおっぱい山です。誘われたら断れません。
一見、変哲もなくつまらない。それでいいんです。
言わんとすることは分かるんです。でも危険と裏腹な恐ろしさがあります。
まとめ
地理院タイルとRで、おっぱい山の可視化が簡単にできるようになりました。
しかしながら、我々が為すべきことは山積みです。
どのようにして、おっぱい山ができたのか？
なぜ、おっぱい山はそこにあるのか？
これから、おっぱい山はどうなってゆくのか？
おっぱい山の頂き目指し、これからも一歩一歩登っていきたいと思います。
それでは、良いクリスマスを！
バイナリベクトルタイルをRで表示する
これはベクトルタイル Advent Calendar 2017 16日目の記事です。
前回、QGISのVector Tiles Reader Pluginでバイナリベクトルタイルを表示してみましたので、今回は、そのコードを参考にRでバイナリベクトルタイルを表示できるようにしてみます。
できあがりイメージ
OpenStreetMapデータの場合
（hfuさんのOSM日本ベクトルタイルデータを使わせてもらいました。）
植生データの場合
使い方
とりあえず使い方です。
コード全部を取ってきて、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
