Cloud Optimized Point Cloud(COPC)の作成と表示方法

ドローン画像のフォトグラメトリで作成した点群データをブラウザで快適に表示できるCOPC形式の点群データの作成と表示をやってみました。

copc.io

las、lazファイルの用意

地理座標系(緯度、経度)以外で書き出したlasファイルもしくはlazファイルの点群データを用意します。高度は楕円体高です。(標高になっている場合は補足で説明している方法で変換します。)

COPCへの変換

変換方法は2通りあります。

その1

QGIS3.26以降を用意します。las、lazファイルを画面にドラッグするとxxx.copc.lazのようなファイルが出来上がります。今のところファイル名に日本語は使えません。

その2

PDALを使用します。WindowsならQGIS3.26をインストールした際に入るOSGeo4Wで利用できます。
OSGeo4Wで以下のコマンドを実行しするとlazをcopcに変換できます。

pdal translate -i xxx.laz -o xxx.copc.laz --writers.copc.forward=all

COPCファイルのアップロード

どこかのWEBサーバーにCOPCに変換したファイルをアップロードします。WEBサーバーではCORSを受け付けるように.htaccessを以下のように設定しておきます。

Header set Access-Control-Allow-Origin "https://viewer.copc.io"

COPCデータの表示

https://viewer.copc.io
にアクセスしてデータボタンからAdd point cloudを押してcopcファイルのURLを貼り付けます。
もしくは以下のように直接URLを指定してアクセスします。
https://viewer.copc.io/?copc=https://www.example.com/xxx.copc.laz

  • Ctrlを押しながらドラッグで回転できます。
  • 色はデフォルトでElevationになっているのでパネルからRGBに変更すると実際の色で表示できます。
  • 点のサイズを変えたり、表示数(Point budget)を変えたりすると見た目の詳細さが変わります。

COPC Viewerについて

以下によるとオープンソースではないようです。見た目はなんとなくCesiumっぽいので、そのCOPC対応版なのかな?と思ってみたり。
https://hobu.co/copc-viewer.html

補足 lazが標高になっている場合

手元にある点群データが平面直角座標系10で高さは標高値となっている場合の変換方法です。COPC Viewerの地形データは楕円体高となっているため高さを合わせるためにはcopcのデータも楕円体高に変換しておく必要があります。(点群データの高さは一般的には楕円体高なんですかね?知りません)

PDALのパイプライン用のファイルを以下のように作成します。

  • xxx.lazは読み込むファイル
  • 最初のin_srsは平面直角座標系10(epsg:6678)のproj指定にgeoidgridを追加したもの。geoidgridで指定しているgtxファイルはこちらで説明しているものを使用。
  • 最初のout_srsはEPSG:6667(緯度、経度、楕円体高)。ここで一旦、平面直角座標系10の標高を緯度経度の楕円体高に変換する。
  • 次のout_srsはEPSG:6678(平面直角座標系10)。これで楕円体高の平面直角座標系10になる。ただ、EPSG:6678は2D用なので3D用には別のものがある気がするけど分かりません。(EPSG:6678+224みたいな感じ??)
  • 最後にwriters.copcでcopcに変換します。

pipeline.json

[
{
"type" : "readers.las",
"filename" : "xxx.laz"
},
{
"type":"filters.reprojection",
"in_srs":"+proj=tmerc +lat_0=40 +lon_0=140.833333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +geoidgrids=D:/PDAL/jp_gsi_gsigeo2011.gtx",
"out_srs":"EPSG:6667"
},
{
"type":"filters.reprojection",
"in_srs":"EPSG:6667",
"out_srs":"EPSG:6678"
},
{
"type" : "writers.copc",
"filename" : "xxx.copc.laz"
}
]

pdal pipeline pipeline.json で変換を実行します。

Phantom4 RTKの飛行ログを取得する方法

1. プロポとPCをUSBで接続します。
2. GS RTK アプリを起動して、画面を上からスワイプします。
3. USB ストレージを接続する ボタンを押します。
4. PCのエクスプローラでD:\DJI\dji.prof.mg.gsp\LOG\tcp_769にアクセスします。
5. フォルダ内の飛行ログをコピーします。
6. DJI Flight Log Viewer | Phantom Help このサイトに飛行ログをアップロードしてKMZファイルに変換します。

標高demデータから楕円体高のデータをgdalで作成する方法

なんだか楕円体高のdemデータが欲しくなるときがあります。
そんな時のために取得方法を紹介します。


2022.8.19追追記
Phantom4 RTKの地形認識モード使用時にPPKの場合(D-RTK 2を使用しない)は、高度は気圧計で算出されるため、標高やその日の天気で誤差が大きくなってしまうようです。で、やっぱり楕円体高でDEMを用意する必要があるようです。

https://forum.dji.com/forum.php?mod=redirect&goto=findpost&ptid=264115&pid=2734061

ジオイド高の分だけ高く飛んだというのは、偶然それぐらいの誤差が出ていただけだった可能性が高いので取り消し

2021.10.19追記
Phantom4 RTKの地形認識モードのDSMは、普通に標高で良いのかもしれません....実際に飛ばしたら、ジオイド高の分だけ高く飛んでいました...

1.

QGIS 3.20をインストールして最新のgdal(とproj)を使えるようにします。

QGISのダウンロード

2.

C:\Program Files\QGIS 3.20.3\share\proj に入っている jp_gsi_gsigeo2011.tifを作業フォルダにコピーします。

3.

OSGeo4W shellを起動して作業フォルダに移動します。

4.

日本のgeoidモデルのtifをgtx形式に変換します。

gdal_translate -ot Float32 jp_gsi_gsigeo2011.tif jp_gsi_gsigeo2011.gtx

ジオイドモデルの元データは、国土地理院で提供されています。使用に関する許可申請等は国土地理院のサイトおよびprojファイルに入っているjp_gsi_README.txtを確認ください。

5.

jp_gsi_gsigeo2011.gtxをC:\Program Files\QGIS 3.20.3\share\projのフォルダにコピーします。

6.

標高demデータを緯度経度(epsg:6668)で用意します。

例えばこれで作成
基盤地図情報 標高DEMデータ変換ツール | コンテンツ | 株式会社エコリス

7.

gdalwarpで楕円体高になるようにパラメータを指定して変換します。

gdalwarp -s_srs "+init=epsg:6668 +geoidgrids=jp_gsi_gsigeo2011.gtx" -t_srs epsg:6667 merge.tif merge_ellipsoidal.tif

※ Warning 1: +init=epsg:XXXX syntax is deprecated. と出ますが気にしません。

補足

  • こんなことをしなくても gdalwarpで -s_srs をepsg:6697 (標高)-t_srsをepsg:6667(楕円体高)と指定すれば変換できると思ったのですが駄目でした。なんで駄目かを誰か教えてください!
  • -co TFW=YES を追加すると.tfwも出力できます。
  • これはDEMなので、DSMとして利用する場合は、建物とか木の高さを考慮して飛行高度を設定してください。
  • 水域などnodataの場所がある場合は気をつけてください!(nodataを0とかにしてるとジオイド高が値に入って、それっぽくなってしまいます。そのままドローンを飛ばすと水域にダイブすることになります。)

Phantom4RTKのPPK処理

Phantom4RTKで撮影したデータをPPK(Post Processing Kinematic)するためのメモです。

処理の概要

  • Phantom4RTKではGNSSの2周波の観測データを取得できます。対応しているGNSSは、GPS(米)、GLONASS(露)、BeiDou(中)、Galileo(欧)です(みちびき(日)は未対応)。
  • 観測データおよび航行データは、0.2秒間隔(5Hz)で取得して、RTCM3.2 MSM5 形式のファイル「PPKRAW.bin」に保存されます。また、観測データをRINEX形式に変換したファイル「Rinex.obs」も保存されます。(PPKRAW.binを、FormatにRTCM3を指定しrtkconvで変換すると、Rinex.obsと同様の観測データのほかに、航行データが取得できます。)
  • また、各画像の撮影時刻と、観測データを受信したアンテナの位置と画像を撮影したCMOSセンサーの位置の差は「Timestamp.MRK」に保存されます。
  • 観測データをRTKLIBで処理することで正確な位置を算出し、CMOSセンサーとの位置の差を考慮して撮影画像に正確な位置を付与します。
  • その際、観測データの位置情報は0.2秒間隔で、撮影位置とは一致していないので、観測データの取得時刻と撮影時刻の関係から線形補間します。

RTKLIBで正確な位置を算出

移動局(Rover)のデータにはRinex.obsを使用し、基準局(Base)と航行データ(.nav)は、近くの電子基準点のデータをダウンロードして使用します。RTKLIBの処理の詳細は、以下を参照してください。

ポイントとしては、次の処理のためにRTKLIBののOptionsのOutputの設定でHeaderをOFFにして、Time formatをww ssss GPST(1980年1月6日からの週数と週始めからの秒数)にして、Field separatorを「,」にしておきます。また、HeightはGeodeticにして出力を標高値にしておきます。(こうすることで、後述のMetashapeで処理した後のdemデータが、標高値になります(未実施なので要確認)。←おそらく問題ないと思いますが、楕円体高で処理した場合とで精度に違いが出ないか少し気になってます。)

Phantom 4 RTK - PPK Processing Workflow | Drone Data Processing
https://swest.toppers.jp/SWEST21/program/pdfs/s5a_public.pdf
https://www.naro.affrc.go.jp/publicity_report/publication/files/drone_gnss.pdf

撮影位置を補間

RTKLIBの処理で作成した.posファイル(0.2秒間隔の位置情報)と撮影時刻を記録した「Timestamp.MRK」から撮影位置を算出します。
算出には、以下からダウンロードできるEXCELファイルを利用します。
使用方法は、以下のサイトとファイルに書いてある通りですが、position fileの貼り付けは、別のシートで.posファイルをカンマ区切りで読み込んでおいて貼り付けると簡単です。(数式の列番号が間違っているところがあるので、自分で修正が必要です。)

Phantom 4 RTK - PPK Processing Workflow | Drone Data Processing

撮影位置を画像に付与

EXCELから出力されたcsvファイル(画像のファイル名と位置情報の対応情報)をMetashapeでインポートすれば、画像に位置情報を付与できます。
以下の資料にもあるように解の精度(Fix,Float,Single)によって位置精度を変更するように改造すると良いかもしれません。(以下資料のTimestamp Converterは公開されていない?ようです)

以下の12ページ目参照
Phantom4 RTKで、後処理キネマティック(PPK)をする方法


精度の目安は、fixで水平0.03m 垂直0.05m floatで水平、垂直0.5mにしようと思います。あと、画像がサブフォルダに分かれていて同じファイル名が複数あるとcsvと対応づけできないので、Metashapeに読み込む前にファイル名をフォルダ名+ファイル名に変換しておきます。

参考バッチファイル
フォルダ名をファイル名の頭に追加するバッチ| OKWAVE

RTKGPS+(カスタムバージョン)の更新

以前ビルドしたものがpixel3で起動できなくなってたので修正しました。

修正したapkはこちらからダウンロードできます。

コードはこちら
github.com

以下作業メモ

以前修正したコードを取ってくる。
git clone -b Customized https://github.com/tmizu23/RtkGps.git

RTKLIBの取得(最新版に変更)

cd RtkGps
git submodule deinit -f jni/RTKLIB
git rm -rf jni/RTKLIB
rmdir /s .git\modules\jni\RTKLIB
git submodule add -b rtklib_2.4.3 https://github.com/tomojitakasu/RTKLIB.git jni/RTKLIB

jni/RTKLIB-unixsocket.patchのパッチをあてる
android studioVCS -> Apply patch... でファイルを指定する

以前の記事に従いソースを修正

simonlynen_android_libsの取得

git submodule update -i

build.gradle

コメントアウト
//ndkVersion '21.3.6528147'

変更
classpath 'com.android.tools.build:gradle:4.0.1'

local.properties

追加
ndk.dir=C\:\\Users\\mizutani\\AppData\\Local\\Android\\Sdk\\ndk\\21.3.6528147

以下のメッセージに従う

Minimum supported Gradle version is 6.1.1. Current version is 4.4.

Please fix the project's Gradle settings.
Fix Gradle wrapper and re-import project

AndroidManifest.xml

削除
android:minSdkVersion="18"
android:targetSdkVersion="18"

基盤地図情報 標高DEMデータ変換ツールのコンパイル方法

基盤地図情報 標高DEMデータ変換ツール Version1.7.0コンパイル方法です。
(QGIS3.10以降だと、以前のバージョンで変換したdemのCRSが不明となってしまったので、gdalとprojを更新して対応しました。ついでに32bit OS での動作は非対応にしました。その際の自分用メモです。

1. Visual Studio Community 2017をインストール

https://visualstudio.microsoft.com/ja/vs/older-downloads/
ここからダウンロードしてインストールして、x64 Native Tools コマンドプロンプトを起動します。

3.gdalをコンパイル

gdalの3.0.4をダウンロードします。
https://github.com/OSGeo/gdal/releases/download/v3.0.4/gdal-3.0.4.tar.gz

nmake.optを変更します。
・nmake.optのSETARGVを変更。
SETARGV = "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.16.27023\lib\x64\setargv.obj"

・nmake.optのPROJを変更
PROJ_INCLUDE = -IC:\Users\mizutani\Desktop\proj-6.3.2\build\distro\include
PROJ_LIBRARY = "C:\Users\mizutani\Desktop\proj-6.3.2\build\distro\lib\proj.lib"

・nmake.optのssqlite3を変更
SQLITE_INC=-IC:\Users\mizutani\Desktop\proj-6.3.2\build\sqlite3
SQLITE_LIB=C:\Users\mizutani\Desktop\proj-6.3.2\build\sqlite3\sqlite3.lib

gdalを以下のコマンドでビルドします。
nmake -f makefile.vc MSVC_VER=1910 WIN64=yes
nmake -f makefile.vc devinstall MSVC_VER=1910 WIN64=yes

c:\warmerda 以下にgdal一式がインストールされます。

4. dem.cppをコンパイル

基盤地図情報 標高DEMデータ変換ツールの中に移動して以下を実行します。
cl -Ic:\warmerda\bld\include c:\warmerda\bld\lib\gdal_i.lib dem.cpp

5. gdalのexe、dll、dataをコピー

C:\warmerda\bld\binから基盤地図情報 標高DEMデータ変換ツールの中にコピーします。
・gdal300.dll
・gdalbuildvrt.exe
・gdalwarp.exe
・gdaldem.exe

C:\warmerda\bld\dataフォルダを基盤地図情報 標高DEMデータ変換ツールの中にコピーします。

6. projをコピー

C:\Users\mizutani\Desktop\proj-6.3.2\build\distro\share\projフォルダを基盤地図情報 標高DEMデータ変換ツールの中にコピーします。

以上です。

補足

変換結合.vbsの中でgdalを使用するための 環境変数を設定してあります。

tempEnv.Item("GDAL_DATA") = "data"
tempEnv.Item("GDAL_FILENAME_IS_UTF8") = "YES"
tempEnv.Item("PROJ_LIB") = "proj"