Hatena::ブログ(Diary)

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

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

2017-01-25 360°パノラマ画像の作り方 このエントリーを含むブックマーク

追記 2017/2/15

huginのパノラマ結合はいまいちだったので、Image Composite Editorで「sphere」を指定してパノラマを書き出して、それをペイントとかで縦横比1:2のキャンバスに貼り付ける方法の方がいいかもしれません。(水平線が中央になるように貼り付ける。空と地面を撮影していない場合は上下は空白になります。)


ドローンで撮影した画像を合成して、360°パノラマ画像をブラウザで表示する方法を紹介します。

※これまではMicrosoftのImage Composite EditorPhotosynthで出来たのですが、残念ながら2017年の2月でPhotosynthがサービス終了となるので、その代替手段です

パノラマ画像の作成

写真を合成してパノラマ画像を作成するためにHuginというソフトを利用します。

http://hugin.sourceforge.net/

インストールして、起動したらアシスタントメニューに従って作業をするだけですが、「パノラマを作成...」の前に、別のウインドウのメニューにあるステッチャーの作成で以下のように設定します。

  • 投影法をEquirectangular(正距円筒図法)にする。
  • 画角を水平方向:360、鉛直方向:180にする
  • 切り抜きをキャンパスサイズに合わせる

出力ファイル名は、panorama_fused.tifとしておきます。

パノラマ画像をタイル化

合成したパノラマ画像をタイル化します。

以下のプログラムを保存します。

https://github.com/mpetroff/pannellum/blob/master/utils/multires/generate.py

pythonのコンソールで、以下のように実行します。

python generate.py -n "C:\Program Files\Hugin\bin\nona.exe" panorama_fused.tif

outputフォルダにタイル化された画像ができます。

html作成

ブラウザで表示するために、以下のライブラリを使って、htmlを作成します。

https://pannellum.org/

こんな感じです。

<!DOCTYPE HTML>
<html>
<head>
    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <title>仙台市 井土浦(2016年3月4日ドローンにより撮影)</title>
    <link rel="stylesheet" href="https://cdn.pannellum.org/2.3/pannellum.css"/>
    <script type="text/javascript" src="https://cdn.pannellum.org/2.3/pannellum.js"></script>
    <style>
    #panorama {
        width: 100%;
        height: 500px;
    }
    h1 {
       text-align:center;
       font-size:14pt;
    }
    </style>
</head>
<body>
<H1>仙台市 井土浦(2016年3月4日ドローンにより撮影)</H1>
<div id="panorama"></div>
<script>
pannellum.viewer('panorama', {
    "type": "multires",
    "multiRes": {
        "basePath": "output",
        "path": "/%l/%s%y_%x",
        "fallbackPath": "/fallback/%s",
        "extension": "jpg",
        "tileResolution": 512,
        "maxLevel": 4,
        "cubeResolution": 3176
    },
    "autoLoad": true,
    "compass": true,
    "northOffset": 0,
    "pitch": -17
});
</script>

</body>
</html>

この部分の設定は、outputフォルダのconfig.jsonにあるのでコピペすればOKです。

        "tileResolution": 512,
        "maxLevel": 4,
        "cubeResolution": 3176

できあがり

http://www.ecoris.co.jp/introduction/drone/panorama.html

ツアーバージョン

http://www.ecoris.co.jp/introduction/drone/natori/panorama.html

2016-12-20

おっぱい山をFOSS4Gでアーカイブする方法

11:57 | おっぱい山をFOSS4Gでアーカイブする方法を含むブックマーク

この記事はFOSS4G Advent Calendar 2016の20日目の記事です。


今年のFOSS4G Tokyoの特別セッションは「アーカイブ×FOSS4G」というテーマで、大変刺激を受けました。

ということで私もそれにのっかって、おっぱい山をFOSS4Gでアーカイブする方法を紹介したいと思います。

f:id:tmizu23:20161219114950p:image
仙台市上空よりおっぱい山を望む(巨乳化処理済み)

1. おっぱい山を作る

公園の砂場に行き、おっぱい山を作ります。おっさん一人で砂場で遊んでいると通報されますので、自分の子か親戚の子と作業してください。

f:id:tmizu23:20161204110446j:image:w300f:id:tmizu23:20161204111206j:image:w300
作業の様子おっぱい山

2. 写真を撮る

おっぱい山が出来上がったら、UAV(ドローン)を用いた公共測量マニュアル(案)に従って写真を撮ります。ドローンはありませんので、スマホを手に砂場の周りをぐるっと回ってシャッターを押してください。注意点は、写真に位置情報が記録されるよう設定しておくこと、7割以上オーバーラップさせて撮影すること、公園にいるママさん達に不審がられないようにすることです。

f:id:tmizu23:20161204110120j:image:w300f:id:tmizu23:20161214230147p:image:w300
撮影の様子。自分の影が写らないように。撮影したおっぱい山(80枚)


3. OpenDroneMapで3Dモデルとオルソ画像を作る

写真が準備できたら、OpenDroneMapで3Dモデルとオルソ画像を作成します。OpenDroneMapは、ドローンで撮影した空中写真から地形の3Dモデルとオルソ画像を作ることができるオープンソースのソフトウェアです。*1ドローンだけでなく、スマホで撮影した写真からでも作成できます。


3-1. Bash on Windowsの導入

OpenDroneMapをWindowsで使うためにBash on Windowsを導入します。Bash on Windowsとは、Windowsで動作するBashシェルです。これによってOpenDroneMapのようなLinux用のソフトもWindowsでコンパイル&実行できるようになります。導入方法はここでは説明しないので、適当にググってください。


3-2. OpenDroneMapのコンパイル

Bash on Windowsを起動して、以下のコマンドでGitHubからOpenDroneMapのソースをダウンロードします。

git clone https://github.com/OpenDroneMap/OpenDroneMap.git

~/.bashrcにパスを追加して、sourceコマンドで設定を読み直します。/your/path/の部分はソースのある場所に置き換えてください。

export PYTHONPATH=$PYTHONPATH:/your/path/OpenDroneMap/SuperBuild/install/lib/python2.7/dist-packages
export PYTHONPATH=$PYTHONPATH:/your/path/OpenDroneMap/SuperBuild/src/opensfm
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/your/path/OpenDroneMap/SuperBuild/install/lib
source ~/.bashrc

OpenDroneMapのフォルダに移動してコンパイルします。依存するライブラリなど諸々インストールされるので結構時間がかかります。(1時間ほど)

cd OpenDroneMap
bash configure.sh

3-3. GCP(グラウンドコントロールポイント)の作成

OpenDroneMapを実行する前に、写真の位置あわせをするためのGCPを作成しておきます。スマホの写真には位置情報が記録されているのでGCPがなくてもオルソ画像は出来ますが、精度をあげるために作成します。GCPファイルは、gcp_list.txtという名前で以下のフォーマットで作成します。

coordinate system description
x1 y1 z1 pixelx1 pixely1 imagename1
x2 y2 z2 pixelx2 pixely2 imagename2
x3 y3 z3 pixelx3 pixely3 imagename3

coordinate system descriptionは位置の投影法、x1,y1,z1は位置座標、pixelx1,pixely1は写真のピクセル座標、imagenameは写真ファイル名です。1枚の写真に対して最低5点ぐらいはあると良いようです。具体的には、こんな感じです。

WGS84 UTM 54N
494399.2 4232697.1 5 333 1590 IMG_6785.jpg
494399.2 4232697.1 5 758 1698 IMG_6786.jpg
494399.2 4232697.1 5 877 1942 IMG_6787.jpg

手動でこのファイルを作ろうとすると結構手間ですが、支援ツールもあるようなので、利用してみると良いかもしれません。*2


3-4. OpenDroneMapを実行する

プログラムを実行する準備として、OpenDroneMapフォルダの直下にoppaiyamaフォルダを作成します。さらにoppaiyamaフォルダの中にimagesフォルダを作って、そこにスマホで撮影したおっぱい山の画像を入れておきます。gcp_list.txtはoppaiyamaフォルダの直下に入れます。

その後、以下のコマンドでプログラムを実行します。写真の枚数や設定にもよりますが、80枚で30分ぐらいはかかるので、気長に待ちます。

sh run.sh --project-path oppaiyama --odm_georeferencing-useGcp  --odm_orthophoto-resolution 1000

実行パラメータは、色々ありますが、今回は、gcpを使用して、オルソ画像の解像度を1000pixcel/mと設定しておきます。

パラメータの詳細は、ここを参照してください。https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters


3-5. データを確認する

プログラムが終了したら、oppaiyamaフォルダの中に3Dモデルとオルソ画像が出来上がっています。3Dモデルは、odm_texturingフォルダの中に「odm_textured_model_geo.obj」という名前で入っているので、MeshLabなどで読み込んで確認してみましょう。また、オルソ画像は、odm_orthophotoフォルダの中に「odm_orthophoto.tif」という名前で入っているので、QGISに読み込んで確認してみましょう。

f:id:tmizu23:20161216090006p:image:w300f:id:tmizu23:20161216090007p:image:w300
MeshLabで3Dモデルを表示QGISでオルソ画像を表示


4. OpenDroneMapのデータをCesium用に変換する

OpenDroneMapで作成した地形データをCesiumで表示するためデータを変換します。Cesiumとは、地図をブラウザで3D表示できるオープンソースのjavascriptライブラリです。Cesiumの使い方は、このあと説明しますが、とりあえず読み込むためのデータを作成しておきます。

CesiumにOpenDroneMapのデータを読み込む方法は二通りあります。一つは、objファイルをgltf形式の3Dモデルに変換して読み込む方法です。もう一つは、lasファイル(点群データ)をCesium用の標高タイルに変換して読み込む方法です。


4-1.(方法その1) objファイルをgltf形式の3Dモデルに変換する

おっぱい山を3Dモデルとして読み込むために、まずmeshlabでobjファイルをcollada(dae)形式に変換します。

odm_texturingフォルダの中のodm_textured_model_geo.objをMeshLabに読みこんでcollada(dae)形式でエクスポートします。その際、MeshLabでいらない部分を消すなど編集しておくと良いかもしれません。ファイル名は「oppaiyama.dae」としておきます。

次に、collada(dae)形式をgltf形式に変換するためcollada2gltfというソフトをGitHubからダウンロードしてコンパイルします。*3

git clone https://github.com/KhronosGroup/COLLADA2GLTF.git
cd COLLADA2GLTF
cmake .
make
sudo make install

以下のコマンドで、oppaiyama.daeをoppaiyama.gltfに変換します。-eはテクスチャ等も同梱するオプションです。

collada2gltf -e -f oppaiyama.dae -o oppaiyama.gltf

f:id:tmizu23:20161218223226p:image

collada2gltfが簡単に利用できるのもBash on Windowsのおかげです。slコマンドも動くよ!


4-2.(方法その2)lasファイルを標高タイルに変換する

おっぱい山を標高タイルとして読み込むために、LastoolsでlasファイルをDEMデータに変換します。まず、Lastoolsを以下からダウンロードます。

http://www.cs.unc.edu/~isenburg/lastools/

Lastoolsの一部のツールは残念ながらフリーウェアではありません。今回使用するlas2dem.exeも商用使用だと有料になりますが、個人的に使用する分にはサイズ制限はありますが使用できるようです。*4ソフトを解凍したら以下のコマンドで変換します。ただし、las2dem.exeはWindows用の実行ファイルなので、Bash on Windowsではなく、別途、PowerShellを起動して実行します。変換するlasファイルは、odm_georeferencingフォルダに入っているodm_georeferenced_model.ply.lasを使用します。

cd Lastools\bin
las2dem.exe -i odm_georeferenced_model.ply.las -o dem.tif -elevation -step 0.001 -utm 54N -nodata 0

オプションの設定は、高さ方向0.001mの解像度で、UTM54の投影法、nodataは0の標高データをdem.tifで出力するとなっています。

次に、DEMデータをCesium用の標高タイル(heightmap-1.0 terrain format)*5に変換するため、Cesium Terrain Builderをコンパイルします。Cesium Terrain Builderのコンパイルには、gdal2.0.0以上が必要なので、以下のコマンドでgdalをアップデートしておきます。(Bash On Windowsでの作業に戻っているので間違えないように)

sudo add-apt-repository ppa:ubuntugis/ppa && sudo apt-get update
sudo apt-get install libgeos-c1v5 gdal-bin

gdalをアップデートできたら、Cesium Terrain Builderのソースをダウンロードしてコンパイルします。

git clone https://github.com/geo-data/cesium-terrain-builder.git
cd cesium-terrain-builder
mkdir build && cd build && cmake ..
sudo make install

ソフトは準備できましたが、今回のデータは対象範囲が狭すぎるためそのままでは上手く変換できません。そのため、Cesium Terrin Builderで変換する前にgdal_translateで巨乳化処理をしておきます。

gdal_translateの-a_ullrオプションで範囲を1000倍に拡大し、-scaleオプションで高さも1000倍にして最低標高を調整しておきます。ついでに-projwinオプションで必要な範囲にクリップします。(nodataを範囲に含むと変換結果がおかしくなる?ので切り取る)

gdal_translate -projwin 494398.737 4232696.642 494400.030 4232694.956 -scale 4.740 5.160 -80 340 -a_nodata -9999 -a_ullr 494398.737 4232696.642 495893.03 4231008.956 dem.tif  dem_scale.tif

Cesium Terrain Builderのctb-tileコマンドで、DEMデータを標高タイルに変換します。標高タイルはterrainフォルダに保存します。

mkdir terrain
ctb-tile -o terrain dem_scale.tif

標高タイルには以下の定義ファイルも必要なので、「layer.json」という名前でterrainフォルダの中に入れておきます。

{
  "tilejson": "2.1.0",
  "format": "heightmap-1.0",
  "version": "1.0.0",
  "scheme": "tms",
  "tiles": ["{z}/{x}/{y}.terrain"]
}

また、0/1/0.terrainと0/0/0.terrainは必須(?)のようなので、0/0/0.terrainを、そのまま0/1/0.terrainとしてコピーしておきます。(今回のデータ範囲は、0/0/0の領域ではないので作成されていないため)

オルソ画像も地図タイルにしておきます。dem.tifと同じように範囲を1000倍に拡大して、png形式の地図タイルに変換しておきます。なお、ctb-tileで作成される地図タイルはtmsなので原点は南になります。

gdal_translate -projwin 494398.737 4232696.642 494400.030 4232694.956 -a_ullr 494398.737 4232696.642 495893.03 4231008.956 odm_orthophoto.tif ortho_scale.tif
mkdir ortho
ctb-tile -p mercator -f png -o ortho ortho_scale.tif

ctb-tileのオプションの詳細は、こちらにあります。https://github.com/geo-data/cesium-terrain-builder


5. Cesiumでおっぱい山を表示する

おっぱい山のデータが出来上がったので、これをCesiumで表示してみます。

5-1. Cesiumを用意する

f:id:tmizu23:20161216140657p:image:w360:right

Cesiumのサイトからライブラリをダウンロードします。

https://cesiumjs.org/

解凍するとファイルとフォルダが色々入っていますが、使用するのはAppsフォルダの中のHelloWorld.htmlとBuildフォルダの中のCesiumフォルダ一式になります。

ひとまず動作確認のために、Bash on Windowsで以下のコマンドを打って、簡易ウェブサーバーを起動します。その後、ブラウザで http://localhost:8000/ にアクセスして「Hello World」のリンクをクリックしてみましょう。Ctrl+Cでサーバーを終了できます。

cd Cesium-1.28
python -m SimpleHTTPServer

5-2. gltfの読み込み

gltfを読み込む場合は、以下のコードをHelloWorld.htmlに追加します。3Dモデルを正しい方向で読み込むには、中心を軸に-90度回転させる必要があります。高さの中心は、今回は見た目で微調整しました。Cesiumの高さの基準は楕円体面(海面0mではない)なので、そのへんの調整も必要です。

    var center = Cesium.Cartesian3.fromDegrees(140.9360018,  38.2422495, -5.5);//地形を読み込む場合は42.5にする
    var heading = 0.0;
    var pitch = 0.0;
    var roll = -Math.PI / 2.0;
    var hpr = new Cesium.HeadingPitchRoll(heading, pitch, roll);
    var transform = Cesium.Transforms.headingPitchRollToFixedFrame(center, hpr);
    var model = viewer.scene.primitives.add(Cesium.Model.fromGltf({
         url:"./oppaiyama.gltf",modelMatrix : transform, scale : 1.0})
    );

5-3. 標高タイルとオルソ画像の読み込み

標高タイルとオルソ画像の地図タイルを読み込む場合は、以下のコードをHelloWorld.htmlに追加します。オルソ画像の地図タイルは、南が原点となっているので{reverseY}と指定します。

var terrainProvider = new Cesium.CesiumTerrainProvider({
    	url : './terrain'
});
viewer.terrainProvider = terrainProvider;

var ortho = new Cesium.UrlTemplateImageryProvider({ 
    url: './ortho/{z}/{x}/{reverseY}.png'
});
viewer.scene.imageryLayers.addImageryProvider(ortho);

5-4. おっぱい山の表示

ついにCesiumで表示できるところまで来ました。5-1.と同様に簡易ウェブサーバーを起動して確認したいところですが、gzip形式で圧縮されている標高タイルは、そのままではSimpleHTTPServerでは配信できません。そこで、元のSimpleHTTPServer.pyをコピーして以下のコードを加えます。こうすると拡張子が.terrainの標高タイルはgzip形式として配信されます。

cp /usr/lib/python2.7/SimpleHTTPServer.py ./Cesium-1.28
self.send_response(200)
self.send_header("Content-type", ctype)
## SimpleHTTPServer.pyの95行目以降に以下の3行加える(send_head関数内)##
base, ext = posixpath.splitext(path)
if ext == ".terrain":
   self.send_header("Content-Encoding", "gzip")
###################################################

なお、通常のウェブサーバーの場合は.htaccessを以下のように記述して設置しておきます。

AddType application/octet-stream .terrain
AddEncoding x-gzip .terrain

SimpleHTTPServer.pyを変更したら簡易サーバーを起動しましょう。5-1.の時と違って、-m オプションは付けずに.pyも含むので間違えないように。

python SimpleHTTPServer.py

http://localhost:8000/Apps/HelloWorld.html にアクセスしてみましょう。(簡易サーバーだと、読み込みに少し時間がかかります。また、時々gltfの読み込みに失敗するので、表示されない場合はShift+リロードしてみます。)

f:id:tmizu23:20161217155138p:image:w600

OpenDroneMapで作成した3DモデルをCesiumで表示

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld1.html

f:id:tmizu23:20161217210020p:image:w600

OpenDroneMapで作成した点群データを標高タイルに変換してCesiumで表示

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld2.html


6. アーカイブをVRで体験

ここまでの作業で、おっぱい山をアーカイブして、Cesiumで表示することができました。Cesiumでは、さらにヴァーチャルリアルティ(VR)表示も可能なので、やってみましょう。以下のようにhtmlを変更すると、画面右下にVRボタンが表示されるようになります。

var viewer = new Cesium.Viewer('cesiumContainer', {
    vrButton : true
});

VRボタンを押すと2画面表示になり、Google Cardboardのようなヘッドセットで見ると臨場感あふれるバーチャルリアリティ体験ができます。

※ VR表示はiPhoneには対応していませんが、PCでVR表示させた画面をリモートデスクトップでスマホに表示すれば可能になります。

f:id:tmizu23:20161218144335j:image:w600

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld.html

まとめ

子供と作った砂場のおっぱい山は、夕日とともに崩れてなくなってしまいますが、このようにアーカイブしておけば、思い出と一緒に永遠に残しておくことができます。

プライスレスなこの活動、FOSS4Gがあれば誰でもできます。皆さんもぜひお試しください!それでは、良いクリスマスを。


追伸:

Google Cardboadは、このあと別のおっぱい山の閲覧に有効活用させていただきました。VRすげーよ。ムフフ。



参考リンク

この記事に使用した画像、コード

https://github.com/tmizu23/cesium_oppaiyama


おっぱい山アーカイブ

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld.html

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld1.html

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld2.html


OpenDroneMap関連

https://github.com/OpenDroneMap/OpenDroneMap

https://github.com/OpenDroneMap/OpenDroneMap/wiki/Running-OpenDroneMap

https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters

http://lists.osgeo.org/pipermail//opendronemap-users/2016-June/000314.html

https://github.com/wolkstein/OpenDroneMap-GCP_LIST.TXT-generator

データ変換関連

http://www.cs.unc.edu/~isenburg/lastools/

https://www.cs.unc.edu/~isenburg/lastools/download/las2dem_README.txt

https://github.com/KhronosGroup/COLLADA2GLTF

https://github.com/geo-data/cesium-terrain-builder

https://github.com/AnalyticalGraphicsInc/obj2gltf

http://www.sarasafavi.com/installing-gdalogr-on-ubuntu.html

http://blog.mastermaps.com/2014/10/3d-terrains-with-cesium.html

http://blog.mastermaps.com/2016/09/creating-tins-in-saga-gis.html

http://blog.mastermaps.com/2016/09/creating-tin-from-raster-dem.html

https://cesiumjs.org/data-and-assets/terrain/formats/heightmap-1.0.html

http://cesiumjs.org/data-and-assets/terrain/formats/quantized-mesh-1.0.html

http://www.pdal.io/

Cesium関連

https://cesiumjs.org/

http://www.slideshare.net/makinux7/cesium3

http://www.slideshare.net/yuki1225/foss4g-2016-hokkaidocesium

http://www.slideshare.net/Ihizaki/cesium-63857910

http://stackoverflow.com/questions/38709743/cesiums-height-value-is-altitude-or-above-sea-level

そのほか

http://meshlab.sourceforge.net/

https://vr.google.com/intl/ja_jp/cardboard/

http://psgsv2.gsi.go.jp/koukyou/public/uav/

http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q10144883875

http://www.osgeo.jp/events/foss4g-2016/foss4g-2016-tokyo/foss4g-2016-tokyo-coreday

*1:商用だとPhotoScanとかPix4Dなどがありますが、これらはウン十万円もする高価なソフトで、個人ではとても手がでません。そのようなツールと同等なものを使えるのは、ほんとFOSS4Gのおかげです。ただ、まだまだ開発途中で商用のものに及ばない部分は沢山ありますが。

*2https://github.com/wolkstein/OpenDroneMap-GCP_LIST.TXT-generator

*3:obje形式から直接gltf形式に変換できるobj2gltfというソフトもありますが、テクスチャが表示できない不具合(?)があったので見送りました。

*4:PDALとかOpenDroneMapでDEMに変換できると良いのですが

*5:Cesium用の標高タイルはheightmap-1.0 形式のほかに、quantized-mesh-1.0 形式というものがあります。こちらの形式は標高データをTINにして、それをタイルに切り分ける方法のようですが、まだ使いやすい変換ソフトが無いようなので試せていません

katkat 2017/04/05 12:45 はじめまして。
いつも参考にさせていただいております。

全くの初心者なのですが、
別記事の「2016-03-18 OpenDroneMapの使い方」では
win10pro のPCでDocker Tool-Boxを使用してなんとか結果表示はできたのですが、
こちらの記事中とおり「Bash on Windows」上でOpenDroneMapで実行させると、
「sh run.sh --project-path oppaiyama --odm_georeferencing-useGcp --odm_orthophoto-resolution 1000」の行で、「sh」のコマンドが不明とエラーが出て先に進めませんでした。
何か当方の設定が間違っているのでしょうか?

以前のDocker Tool-Boxはアンインストールしてしまい、ODMを使用した各種データの取得が不可能になってしまいました。

多少windowsの知識がある程度で、linux等ほぼ解らず手を付けて泥沼化しております。

何か有用な知恵があればとコメントしました。

tmizu23tmizu23 2017/04/07 20:03 コメントありがとうございます。
GitHub上のソースは頻繁に更新されていて、上手くコンパイルできないコミットもあるようです。
(私も最新のコミットでは、上手くコンパイルできませんでした。)
このコミットだと大丈夫でした。
https://github.com/OpenDroneMap/OpenDroneMap/tree/bd23610cec516029de92596ed18397401136d0c5

あと、bash on windowsでもapt updateとapt upgradeをしておいた方が良いと思います。

tmizu23tmizu23 2017/04/08 11:09 もう少し検証したところ、bash configure.sh を2回実行しないと上手くコンパイルできないかもしれません。
あと、実行時のコマンドは「sh run.sh --project-path . oppaiyama」 と指定するように変更になったみたいです。
--project-path にプロジェクトのパス(例だとカレントパス「.」) その後にプロジェクト名(例だと「oppaiyama」)を指定します。