無料で使える航空写真とOpenDroneMapを利用して、あなたの街の3Dデータ（点群）を作成する！ #QiitaEngineerFesta2022

（今回作っていく広域点群データの完成図です）

3Dデータを有効に可視化したい

PLATEAUのデータを利用することで、地図 + 地形データの上に詳細な建物データを配置するなど3Dデータを活用したリッチな表現が比較的簡単に出来るようになりました。

が、3Dモデルはテクスチャ（ビルの外観など）を貼り付けるとWeb地図上ではかなーり動作が重たくなってしまう上、モデルが整備されていない・もしくは整備されているが、詳細な形状がモデリングされておらず四角い箱になっていたり、テクスチャが存在しなかったりする都市だと、こんな感じでちょっと寂しい絵面になってしまいます。

そこに、詳細な色付き点群データを利用して乗っけると一気にリッチな感じになりますね！

詳細な広域3Dメッシュモデルがテクスチャ付きで整備されていれば、それを利用すれば良い話ではあるのですが、全国分整備するにはお金も時間もかかります。

なので、今回は無料で利用できる素材・ツールを使って札幌市の点群データを作ってみたいと思います！

航空写真のダウンロード

地図・空中写真閲覧サービスを利用して札幌市のデータをダウンロードしてみましょう！

リンク先に遷移した後、「検索方法」のタブを空中写真に変更し、西暦2020年（最新の航空写真）を指定しましょう。　
複数年に渡るデータを利用してしまうと、建物が増えていたりなくなっていたり…太陽の位置が異なるため色味が変わってしまったりとフォトグラメトリをする上で不具合が出る（というか、エラーで処理が止まる）ので、必ず単年のデータを利用するようにしましょう。

ダウンロードした画像はこんな感じの空から撮った写真になります。

ど真ん中にあるのが札幌駅ですね。

WebODMをインストールする

以下の記事に全ての詳細が書いています（ありがてぇ）

WebODMを無料で簡単インストールする全工程

ここではdockerとgitがお使いのマシンにインストールされている前提でコマンドだけ紹介します！
以下のコマンドでGitHubからコードをダウンロードして、リポジトリのディレクトリに潜り、shを起動させるだけ！
これだけでWebODMが起動します！

git clone https://github.com/OpenDroneMap/WebODM --config core.autocrlf=input --depth 1
cd WebODM
./webodm.sh start

ちなみに公式サイトのインストール方法はこちらに書かれており、推奨環境は以下となっております。

Latest Generation CPU
100 GB of disk space
16 GB RAM

そのほか、利用可能なメモリの3/4割り当てろとか、CPUは半分くらい割り当てろとか要求されますので、そのようにDocker desktop等から利用メモリとCPUの設定を変更しておきましょう。

WebODMで画像を読み込ませてGCPを入れる

今回は札幌駅周辺で画像を集めてみましたが、11枚ほどダウンロードしてみました。

東京などの大都市だと、このくらいの範囲で40枚とかそのくらい集まるんですが、地方都市では最低限のレベルでしか撮影されていないようです…
データは少ないので、期待する精度が出るかは微妙なところなのですが、とりあえずやってみましょう！

./webodm.sh startでWebODMが起動している状態であれば、ブラウザでhttp://localhost:8000/に接続してみましょう。
するとこんな画面が開くかなと思います。

もしログイン画面が出ている場合は適当に設定してアカウントを作りましょう！
今後WebODMを利用するのに必要になるので、忘れないようにしてくださいね！

上記の画像のような画面になったら、画像を読み込ませてGCPを打っていきましょう。

GCPとは

GCPとは「（正確な）座標がわかっているポイント」のことです。

WebODMはOpenDroneMapという名の通り、ドローンで空撮した個別の画像から、写真測量（フォトグラメトリ）を経て大きな画像を作り、地図上に落とし込むことができるシステムですが、正確な座標がわかっているポイントが存在しないと、フォトグラメトリ時に（当然）歪みが出ます。

これを補正する目的で、空撮の前にいくつかのポイントを事前に測量して正確な座標を割り出し、空撮時にはそのポイントにマーカー（1m × 1mの正方形のものがよく利用される）を置いておくことで、画像に正確な座標を割り当てることができます。

詳しくはグラウンドコントロールポイント (GCP) の重要性についてをご参考ください。

今回は、自分以外が撮影した航空写真を利用するので、GCPを設置させることはできませんが、画像の中で緯度経度が特定できる箇所に関しては、フォトグラメトリの前に割り当ててあげることで、ある程度精度を向上させることが可能です。

国土地理院が公開している地理院地図を利用することで、ある程度正確な座標を入手することができますので、これを利用して緯度経度と標高を入力していきましょう。

ちなみに、注意点として、地理院地図に記載の「標高」は地表面の標高ですので、ビルの屋上などをGCPにしてしまうと、垂直方向が大きくずれてしまいます。
例えば、日本で一番高い建物である「あべのハルカス」の屋上にフォーカスを合わせているのに、標高は16.4mと表示されていますね。

逆に、富士山頂付近にフォーカスを当てると3500m超の値が出てきます。
あくまで、「地表面の標高」の数値であることを頭に入れ、GCPは建物や木々が存在しない地表が見えている箇所を選定するのが良いでしょう。

GCPを設置する

ということで設置していきましょう！

サイドバーに「GCPインターフェース」というものがあるのでこちらをクリックします。

するとこんな画面に遷移するので「Choose images / drag here」と書いている箇所に航空写真をドラッグして画像を読み込ませましょう。

右側の地図画面から地図っぽいアイコンをクリックし、「MAP PROVIDER」の「Custom」の箇所にhttps://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpgと入力し、Applyを押して国土地理院提供の航空写真を読み込ませましょう。

画像をクリックすると画面左側に画像が開くのと同時に、真ん中にオレンジ色のマーカーが表示されますが、Deleteで削除しましょう。
同時に、地図画面に表示されたマーカーも削除しましょう。

特徴のあるわかりやすい箇所を見つけたら、画面中央の+ボタンをクリックし…

その箇所をクリックしてポイントを打ちます。

その後、画面右側の地図画面でも同じことをしましょう。

これだけでは画面右側の地図と左側の画像のどのポイントが対応しているのか分からないので、対応させたいポイント同士を1回ずつクリックし、対応を取ります。
認識されるとポイントが緑色になり、どちらかのポイントにカーソルを合わせるだけで、対応しているポイント同士が大きく表示されるようになります。

これを5箇所程度設置し、なおかつ設置したGCPの箇所が他の画像にも存在している場合、流用してあげましょう。

大変ですががんばりましょう。
（とはいえ、元の航空写真の画質の問題、かつ後から目視でGCPを設定している関係で、どう頑張っても水平に60cm以上、垂直に数メートル誤差が発生するので、気楽にやりましょう。）

最終的に、この箇所にGCPを設置しました。

設置が終わったら画面右上の「EXPORT FILE」をクリックしましょう。

するとこのように設置したGCPの座標と、対応する画像ファイルの名称が記載されたテキストファイルが表示されますので、Saveを押して保存しましょう。

ただし、左から3列目のZ値（高さ方向の値）が全て0になっていますね。
これを与えてあげましょう。

左から1、2列目の値はUTM座標のゾーン54として記載された緯度経度とは全く違うXY座標ですので、これを緯度経度に直してみましょう。
こちらのサイトなどを利用し、座標を入力、変換ボタンをクリックするとこのように緯度経度で座標が表示されますね！

ここから得た座標を利用して地理院地図の左上の検索窓にスペース区切りで座標を入力してエンターを押すと、その箇所に飛ぶことができます。

ここで得られた標高をダウンロードしたテキストファイルに記載していきます。

最終的なテキストファイルは以下のようになりました。

+proj=utm +zone=54 +ellps=WGS84 +datum=WGS84 +units=m +no_defs 
532461.85 4770022.73 10.2 2319.78 2401.14 CHO20201-C6-22.jpg
532461.85 4770022.73 10.2 2318.14 725.93 CHO20201-C6-23.jpg
532461.85 4770022.73 10.2 6560.35 2277.88 CHO20201-C7-22.jpg
534066.31 4769360.44 6.6 5014.07 1816.27 CHO20201-C6-23.jpg
533747.86 4767979.73 9.9 4496.07 2465.46 CHO20201-C6-24.jpg
535119.84 4766065.20 11.5 6773.99 3991.45 CHO20201-C6-25.jpg
531415.46 4770334.12 11.6 4817.95 1745.92 CHO20201-C7-22.jpg
530384.32 4767981.50 17.3 3091.94 4001.24 CHO20201-C7-23.jpg
530384.32 4767981.50 17.3 3083.79 2326.81 CHO20201-C7-24.jpg
529563.96 4766623.87 22.1 1688.12 2902.67 CHO20201-C7-25.jpg
527075.79 4769328.49 12.3 1788.22 2494.89 CHO20201-C8-22.jpg
527075.79 4769328.49 12.3 1787.29 813.75 CHO20201-C8-23.jpg
528310.43 4768042.16 19.1 3857.94 2963.77 CHO20201-C8-23.jpg
528310.43 4768042.16 19.1 3863.39 1253.11 CHO20201-C8-24.jpg

フォトグラメトリを行う

ここまで長かったかと思いますが、やっとフォトグラメトリを行なって点群データを作ることができます。

ダッシュボードに戻り、「画像とGCPを選択」をクリックします。

ダウンロードした画像と編集後のtxtファイルを全て選択しましょう。

その後、設定を変更していきます。

お手持ちのPCスペックによって設定は様々なのですが、僕は以下のあたりだけ変更しました。

3d-tiles: 有効
auto-boundary: 無効
debug: 有効
dsm: 無効
feature-qyality: high
pc-quality: high

「画像のサイズを変更」は「いいえ」に設定した上で、「レビュー」をクリック、そのまま処理開始をクリックしましょう。

おそらく10分程度は処理に時間がかかると思いますので、気長に待ちましょう。
また、スペックが足りないとエラーで処理が停止するので、qualityを下げるなど、工夫をしましょう。

待つこと16分…タスクが完了したようです！

データはいろいろな種類でダウンロードできるのですが、今回は「Point Cloud（3D Tiles）」を選択しましょう。
するとzipファイルがダウンロードされると思います。

点群データを地図上に表示してみる

以下の記事で少し詳しく紹介しているのですが、Cesium.jsという地図上に3Dデータを配置することができるJavaScriptのライブラリがあります。
今回はこれを利用して先ほどダウンロードした3DTilesを表示していきましょう！

ヘッダーが壊れているLASファイルを修正して、XY座標を入れ替えた上でpy3dtilesで3DTilesを作成する

ダウンロードしたzipファイル解凍した上で、解凍後のディレクトリと同じ場所に以下のようなindex.htmlを配置します。

<!DOCTYPE html>
<html>
　<head>
    　　<meta charset="UTF-8">
    　　<title>Cesium</title>
    　　<script src="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Cesium.js"></script>
    　<link href="https://cesium.com/downloads/cesiumjs/releases/1.72/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
    <style>

        #cesiumContainer {
            position: absolute;
            top: 0;
            left: 0;
            height: 100%;
            width: 100%;
            margin: 0;
            overflow: hidden;
            padding: 0;
            font-family: sans-serif;
        }

        html {
            height: 100%;
        }

        body {
            padding: 0;
            margin: 0;
            overflow: hidden;
            heght: 100%;
        }

    </style>
    　</head>

　<body>
<div id="cesiumContainer"></div>
<script>
    const raster_tile = new Cesium.OpenStreetMapImageryProvider({
        url: 'https://cyberjapandata.gsi.go.jp/xyz/pale/',
        credit: new Cesium.Credit('地理院タイル', '', 'https://maps.gsi.go.jp/development/ichiran.html'),
        maximumLevel: 18,
    })
    const viewer = new Cesium.Viewer('cesiumContainer', {
        imageryProvider: raster_tile,
        baseLayerPicker: false,
        geocoder: false,
        homeButton: false
    })

    viewer.camera.setView({
        destination : Cesium.Cartesian3.fromDegrees(140.00, 36.14, 20000000.0)
    })

    const tileset = viewer.scene.primitives.add(
        new Cesium.Cesium3DTileset({
            url: './3d_tiles_pointcloud/tileset.json'
        })
    )
    tileset.style = new Cesium.Cesium3DTileStyle({
        pointSize: 2
    })

    viewer.zoomTo(tileset)

</script>
</body>
</html>

ターミナルを起動し、index.htmlが存在しているディレクトリと同じ階層に移動して、以下のコマンドを実行します。

python -m http.server 3000

その後、ブラウザでhttp://localhost:3000/に接続しましょう。

すると、以下のように期待した通りの場所に点群データが表示されているかと思います！！！

札幌駅周辺を上から見てみると結構いい感じです。

横から見ると建物の側面などが結構残念な感じではありますが、まぁ無料で公開されているあまり画質のよくないデータですし、画像の枚数が圧倒的に足りていないため、こんなもんでしょう…！

終わりに

ということで今回は無料で公開されている航空写真を利用して11枚の画像にGCPを設置し、WebODMでフォトグラメトリ + Cesiumに表示するところをやっていきました！

正直11枚の画像（僕がドローンで空撮していた時代は200枚超の画像を使っていました）で綺麗な点群データが出来上がるのか心配でしたが、案外ちゃんとできてよかったです！

みなさんもぜひやってみてくださいね！