地理院タイルから三次元地形モデルを作ってみる

## 概要
このエントリは、[FOSS4G Advent Calendar 2020](https://qiita.com/advent-calendar/2020/foss4g) 12日目の記事です。
<a href="https://twitter.com/t_mat" targe="_blank">@t_mat</a>と申します。
本職は測量設計会社の環境調査部門の人で、いわゆる「非IT系エンジニア」ですが、プログラミングに興味があり、夜な夜な、その時々に思いついたものを作っています。

今回は、国土地理院タイルデータから広域の三次元地形モデルを作成するツールを作ってみました。
ブラタモリやたまたま視聴した[日本地形連合](http://www.jgu.jp/)の公開シンポの影響で地理・地形の面白さを知り、最近ちょこちょこモデルを作って眺めています。
![EoperFRVoAQ8qlc.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/16dbf4c8-0a16-e83d-c762-0ce56ae98248.jpeg)
いまどきは様々なツールで三次元地形を閲覧することができるので、需要があるかはハテナ？ですが、モデルをファイルに出力できれば、色々と活用できるかもと思い、投稿することにしました。

## 開発環境
開発環境は以下のとおりです。実行にはJRE1.8以降が必要です。
　開発言語：Java1.8
　プロジェクト管理：maven
　依存ライブラリ：JTS、TINFOUR、GSON
　IDE：Eclipse
なお、作成したモデルは、Windows標準のPrint3Dなどでも閲覧できますが、[MeshLab](https://www.meshlab.net/)を使いました。

## 作ったもの
今回、個人的に作っていたクラスを整理し、以下のツールを作成しました。

- TileImageGetter.jar：GeoJsonで定義した範囲の地理院タイルを取得・結合するツール</li>
- PlyCreator.jar：テクスチャ画像とDEM画像からPLY形式の三次元地形モデルを生成するツール</li>
- DEMInterpoler.jar：DEM画像を任意解像度に補間するツール
- TerrainImageProcesser.jar：DEM画像からCS立体画像等を生成するツール
- ImageMul.jar：複数の画像の乗算画像を生成するツール

ソースコードはGithubにアップしました(12/7)。 https://github.com/termat/geoply
また、実行形式JarファイルとサンプルデータをGoogleDriveにおいておきます。
https://drive.google.com/file/d/1QLCffoc1tbwUDIczGlIQ26wvnezcP7P3/view?usp=sharing
 geoply
 ├── TileImageGetter.jar
 ├── PlyCreator.jar
 ├── DEMInterpoler.jar 
 ├── TerrainImageProcesser.jar
 ├── ImageMul.jar
 ├── aso　　　 　(サンプルデータ)
 ├── izu　　 　 　(サンプルデータ)
 ├── toyooka　 　(サンプルデータ)
 └── takao　　 　(サンプルデータ)

## 三次元地形モデルの生成手順
三次元地形モデルの生成手順は以下のとおりです。
QGISでモデル領域のGeojsonを作成し、その領域のDEM画像(標高PNG)と航空写真等を地理院タイルから取得してピクセル情報をXYZRGBデータとし、TINで面を構成してサーフェイスモデルの三次元地形モデルを作成します。
![img01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/8486f172-a1a0-44c1-ed17-430e3afb5ec2.jpeg)

## 各ツールの説明
以下に各ツールの使用方法を説明します。
なお、コマンドラインの引数は、オプション等を付けて指定する形にすべきですが、今回は省略しています。
また、大きなモデルを作る時は-Xmxコマンドなどで、JVMに多めのメモリを割り当てる必要があります。
#### TileImageGetter.jar
TileImageGetterは、地理画像タイルを取得・結合し、画像データを出力するツールです。
以下のコマンドで、設定ファイル（json）に定義した条件の画像が取得できます。

```shell
java -jar TileImageGetter.jar [設定ファイルのパス]
```
設定ファイルの内容は以下のとおりです。
サンプルをお使いになる場合は環境に合わせてパス等を変更する必要があります。

```json
{
  "num":2,   //平面直角座標系のNoを指定
  "geojson":"C:/Data/geoply/aso/aso.geojson",    //取得する領域のgeojson
  "zoom":15,  //取得する画像のズームレベル
  "url":["https://cyberjapandata.gsi.go.jp/xyz/dem5a_png/",  //取得するタイルのurlのリスト
	"https://cyberjapandata.gsi.go.jp/xyz/dem5b_png/",　 //上位URLのタイル画像が無い場合
	"https://cyberjapandata.gsi.go.jp/xyz/dem5c_png/"],  //は下位URLを検索する
  "out":"C:/Data/geoply/aso/",   //出力先
  "name":"dem",  //出力ファイル名のprefix
  "ext":"png"    //取得するタイル画像の拡張子、航空写真の場合はjpg
}
```
本ツールで行っている処理は、下図のとおりで、実際の処理はnet.termat.geo.TileReader.java
こうした処理を行うため、Geojsonを読み込むクラス（GeojsonData.java）、緯度経度と平面直角座標系を変換するクラス（LotLatXY.java）、ラスタ用座標ファイルの処理クラス（WorldFile.java）を作成しています。
なお、座標変換は国土地理院の「[Gauss-Krüger 投影における経緯度座標及び平面直角座標相互間の座標換算についてのより簡明な計算方法](https://www.gsi.go.jp/common/000061216.pdf)」を、ラスタ用座標ファイル（ワールドファイル）は農研機構の「[ラスタデータ用の座標ファイル（ワールドファイル）について](http://cse.naro.affrc.go.jp/takahasi/gisdata/avg/worldraster.htm)」を参考にしました。
![img02.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/a18a69da-ec38-4c20-9154-a25d3a2f47a2.jpeg)
#### PlyCreator.jar
DEM画像とテクスチャとして使用する画像から、PLYファイルを生成するツールです。
モデル化の考え方は以下のとおりで、DEM画像のXYZ座標とテクスチャ画像のRGBをVertexとし、TINにより面情報を作成して、サーフェイスモデルの三次元地形モデルを生成します。
以下のコマンドで、PLY形式のモデルが出力されます。

```shell
java -jar PlyCreator.jar [テクスチャ画像のパス] [DEM画像のパス] [出力ファイルのパス]
```
![img03.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/86911186-ea6c-c32b-c759-44ee9cda3443.jpeg)

サンプルでは、以下のコマンドで、画像の様な三次元モデルが生成されます。

```shell
# DEM画像の取得
C:\Data\geoply> java -jar TileImageGetter.jar ./aso/dem.json
# 航空写真の取得
C:\Data\geoply> java -jar TileImageGetter.jar ./aso/photo.json
# 三次元モデルの生成
C:\Data\geoply> java -jar PlyCreator.jar ./aso/photo.png ./aso/dem.png ./aso/photo.ply
```
![img00.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/0deee9ce-5eb3-49ae-4f78-9a63940f9005.jpeg)

#### DEMInterpoler.jar
DEM画像を任意の解像度にリサイズするツールです。
例えば、航空写真はズームレベル16、DEM画像はズームレベル15など、ズームレベルが異なる場合に、補間によりズームレベル16相当のDEM画像を生成するのに使用します。
DEM（標高値）の補間はTIN補間を使用しています。
以下のコマンドで、参照画像と同じ解像度のDEM画像が生成されます。なお、DEM画像はモデル領域周辺のDEMを複数指定することができます。

```shell
java -jar DEMInterpoler.jar [参照画像のパス] [出力ファイルのパス] [DEM画像のパス]･･･
```
![img04.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/658a5b6a-927b-e44a-b593-4c4a2c92394b.jpeg)

#### TerrainImageProcesser.jar
写真画像だけでモデルを作るのも味気ないので、DEM画像からCS立体図等を生成するツールをついでに作りました。net.termat.geo.image配下の画像処理クラス、net.teramt.geo.tool配下のグラジエント処理クラスを使用して作成しています。
以下のコマンドで、傾斜量図、曲率図、CS立体図が生成されます。なお、[コマンド]はcurve、slope、csが指定できます。

```shell
java -jar TerrainImageProcesser.jar  [コマンド] [DEM画像のパス] [出力ファイルのパス]
```
![img05.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/b54c9a80-9fb7-4078-740f-77775a9eb412.jpeg)

#### ImageMul.jar
GIMPを使えば良いのですが、手っ取り早く乗算画像を作りたかったので作成したものです。
以下のコマンドで指定したファイルの乗算画像が生成されます。なお、乗算する画像は2枚以上指定する必要があります。

```shell
java -jar ImageMul.jar  [出力ファイルのパス] [乗算する画像ファイル]・・・
```
![img06.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/86da345d-1c9e-ab73-c72a-9a23f83f44d9.jpeg)

## モデルの作成例
#### 伊豆半島（サンプル：izuフォルダ）
サンプルファイルでは、以下のコマンドで、画像の様な三次元モデルが生成されます。
なお、植生図は、エコリスさんが公開されている[植生図タイル](https://map.ecoris.info/)を使用させて頂きました。

```shell
# DEM画像の取得
C:\Data\geoply> java -jar TileImageGetter.jar ./izu/dem.json
# 植生図画像の取得
C:\Data\geoply> java -jar TileImageGetter.jar ./izu/vega.json
# DEM画像の補間（※DEM標高値の欠測地域があっため）
C:\Data\geoply> java -jar DEMInterpoler.jar ./izu/std.png ./izu/dem2.png ./izu/dem.png
# 植生図三次元モデルの生成
C:\Data\geoply> java -jar PlyCreator.jar ./izu/vega.png ./izu/dem2.png ./izu/vega.ply
```
以下が生成された伊豆半島の三次元地形モデルです。
竹林（黄色）は東向き斜面に多いようですが、何か要因があるのか気になります。
![izu01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/703dc56a-dc1f-dac7-880d-565769baacf0.jpeg)

#### 高尾山（サンプル：takaoフォルダ）
以下のコマンドで、画像の様な三次元モデルが生成されます。
なお、高尾山では、DEMはズームレベル15、航空写真と標準地図はズームレベル16で取得しているので、DEMInterpolerでズームレベル16相当のDEM画像を作成しています。

```shell
# DEM画像取得
C:\Data\geoply> java -jar TileImageGetter.jar ./takao/dem.json
# 航空写真取得
C:\Data\geoply> java -jar TileImageGetter.jar ./takao/photo.json
# 標準地図取得
C:\Data\geoply> java -jar TileImageGetter.jar ./takao/std.json
# ズームレベル16相当のDEM画像生成
C:\Data\geoply> java -jar DEMInterpoler.jar ./takao/std.png ./takao/dem2.png ./takao/dem.png
# DEM画像からCS立体図を生成
C:\Data\geoply> java -jar TerrainImageProcesser.jar cs ./takao/dem2.png ./takao/cs2.png
# 航空写真とCS立体図の乗算画像を生成
C:\Data\geoply> java -jar ImageMul.jar ./takao/csphoto.png ./takao/cs2.png ./takao/photo.png
# 乗算画像の三次元地形モデルを生成
C:\Data\geoply> java -jar PlyCreator.jar ./takao/csphoto.png ./takao/dem2.png ./takao/csphoto.ply
# 標準地図の三次元地形モデルを生成
C:\Data\geoply> java -jar PlyCreator.jar ./takao/std.png ./takao/dem2.png ./takao/std.ply
```
![takataka00.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/f79dfc00-2f2b-414d-a5b2-1bf8cd6bd1a6.jpeg)
![takataka01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/55a0a16a-fbb7-9ac3-c3ff-7a754b0f8d29.jpeg)
その他、エコリスさんが公開されている[植生図タイル](https://map.ecoris.info/)などを活用すれば、立体植生図的なものも作成できます。
![003.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/9bb5fcf1-06d1-56fc-eb4b-6e1529a867e9.jpeg)

## 最後に
わずか数分で地形モデルが作れるのはオープンソースの文化とオープンデータの取組のおかげです。感謝。
こうしたモデルは、見て楽しむ？だけでなく、二次元不等流計算等の数値解析に使えそうなので、そのうち、試してみたいと思っています。
あと、個人的に気になっているのが、Windows標準のPrind3Dで「見積原価」って表示が出てくるので、ひょっとすると、3Dプリンタで地形モデルを出力できるのかな？ということです。
3Dプリンタを所有していないので、試したことはないのですが、出力できたら面白そうだと思っています。

![003.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/4166e491-dae7-3f64-c0cd-61463a2ae360.png)

## 2020/12/13 追記
12/12放送のブラタモリ（城崎・豊岡：玄武岩）を見て、地質図を取得したいと思い、コマンドを追加しました。
※12/13にGithubとサンプルを更新しました。

#### GeologyGetter.jar
産総研の[20万分の1日本シームレス地質図V2 Web API](https://gbank.gsj.jp/seamless/v2/api/1.2.1/)から（元データが20万分の1ですが）地質図画像を取得します。
以下のコマンドで参照画像と同じ領域の地質図画像を出力します。
なお、[平面座標系No]は、平面直角座標系のNoです（豊岡は兵庫県なので5系）。
また、Web APIで凡例情報が取得できる場合は、出力ファイルと同名のjsonファイルで出力します。
※凡例情報が取得できない場合は、「null」のjsonファイルになります。

```shell
java -jar GeologyGetter.jar  [平面座標系No] [参照画像のパス] [出力ファイルのパス]
```
サンプルでは、以下のコマンドで城崎～豊岡にかけてのモデルが出力されます。
なお、この区域では、凡例情報の取得には失敗するようです。

```shell
# DEM画像取得
C:\Data\geoply> java -jar TileImageGetter.jar toyooka/dem.json
# 航空写真取得
C:\Data\geoply> java -jar TileImageGetter.jar toyooka/photo.json
# DEMに欠測箇所が多いので補間
C:\Data\geoply> java -jar DemInterpoler.jar toyooka/photo.png toyooka/dem2.png toyooka/dem.png
# シームレス地質図画像を取得
C:\Data\geoply> java -jar GeologyGetter.jar 5 toyooka/dem.png toyooka/geology.png
# 航空写真と地質図の乗算画像を生成
C:\Data\geoply> java -jar ImageMul.jar toyooka/photo2.png toyooka/photo.png toyooka/geology.png
# 乗算画像の三次元地形モデルを生成
C:\Data\geoply> java -jar PlyCreator.jar toyooka/photo2.png toyooka/dem2.png toyooka/photo2.ply
```
こんな感じのモデルが出力されます。
手前の紫色の範囲が玄武岩が分布する区域です。なお、玄武洞付近はDEMがないため、内挿で標高を設定しています。
![002.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/5ec09e27-2c28-96a1-a126-06026873f8e5.jpeg)