GRASS GIS で dem(tiff)データに対して測量を行う(富士山の体積を測ろう！)

# 概要
[SENSYN Robotics(センシンロボティクス)](https://www.sensyn-robotics.com/)の衛です。

この投稿は、GRASS GISを使った測量の方法を説明することを目的とする。

最近、会社では、お客様の要望より、ドローンで撮った写真に基づいて[DSM](https://en.wikipedia.org/wiki/Digital_surface_model)[(tif)](https://en.wikipedia.org/wiki/GeoTIFF)データを
作成し、地面にある物体の体積を測量するタスクがあった。

手順は大体以下の通り。
1.測量用、ドローンによるエリア撮影を実施
2.オルソ画像、DSM(tif)データ作成
3.DSMデータに対してGRASS GISを使って測量を行う
　＊この投稿は3.のみを対象にする。
　＊この投稿はターミナルのコマンドを対象にしているが、全てのコマンドは基本、[GUI](https://gis-oer.github.io/gitbook/book/materials/GRASS/GRASS.html)でも同じ目的達成できる。

使ったのは[GRASS GIS](https://grass.osgeo.org/)である。
[GRASS GISインストールする必要がある。](http://grassmac.wikidot.com/downloads)

# 測量原理
一枚写真からでは、被写体の立体的な表現（情報記録）ができませんが、
ドローンが一定的なオーバラップ率で撮った、同じ地面の場所の複数写真を、
[立体視](https://ja.wikipedia.org/wiki/%E7%AB%8B%E4%BD%93%E8%A6%96)あるいは[ステレオマッチング](https://en.wikipedia.org/wiki/Computer_stereo_vision)の原理使って変換をかけると、その場所の地面高度が含まれたデータが作成できます。そのうちの一つはDSMデータです。

DSMデータには、２d（例：経度と緯度）の各セル（タイル）に、該当場所の高さが保存されている。
（それはGRASS GISのRasterデータに相当する。）

富士山を例としする。下記の図の０がある暗いセルは地面（あるいは海抜0m）としする。1~4のセルは高さが徐々に上昇する富士山の部分を示している。
<img width="492" alt="Screen Shot 2020-06-05 at 16.03.54.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/d57bf95c-6fd4-cbd5-f5c4-2dfc91ffe47e.png">

シンプルに、測りたいエリア（ポリゴン）の高さの平均を算出して、
エリア（ポリゴン）の面積x高さの平均で、測量対象の体積がわかる。

また、各セルの体積が計算でき、それを加算すると全体の体積となる。

斜面においてある物体に対しては、「基準平面」を斜面にそって作れば良い。
また、基準平面を一律同じ高さ（例えばDEMの全体の最低点）を設定するなど、
いろいろ測量のニーズによって、工夫できるが、
今回はGRASS GISの使用方法を説明するため、海抜0mより上の部分を全部体積とカウントする。

# GRASS GIS セットアップ
GRASS GISのデータ構造は、location＝＞mapset＝＞各種[raster](https://gis-oer.github.io/gitbook/book/materials/00/00.html#%E3%83%A9%E3%82%B9%E3%82%BF%E3%83%87%E3%83%BC%E3%82%BF)と[vector](https://gis-oer.github.io/gitbook/book/materials/00/00.html#%E3%83%99%E3%82%AF%E3%83%88%E3%83%AB%E3%83%87%E3%83%BC%E3%82%BF)データとなる。それを見るmonitorや各種ツールもある。
基本[GRASS GIS](http://grassmac.wikidot.com/downloads)をダウンロードしたら、インストールしてから開いたら、
ウィザードが開かれる。
![Screen Shot 2020-09-01 at 16.47.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/90629fc4-1540-0593-d742-32bb7ae5a7eb.png)
「2. Select GRASS Location」では、Newを選んで、下記のように入力。
![Screen Shot 2020-09-01 at 16.48.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/dc85e4fb-1497-9882-e61e-e8ba25fad81c.png)
「Create a generic Cartesian system (xy)」を選んで次々とする。
![Screen Shot 2020-09-01 at 16.49.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/2930c3ff-68b6-25f9-0d41-438fb83d3cfd.png)
「Start GRASS session」をクリックする。
![Screen Shot 2020-09-01 at 16.52.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/2011414c-7eab-05e0-b747-d8aefe2506a4.png)


# GRASS GIS を使った測量の詳細

* 日本のelevation tifをダウンロードする。
https://www.gsi.go.jp/kankyochiri/gm_japan_e.html
ダウンロードしたファイル名：gm-jpn-el_u_1_1.zip
tifファイル名：el.tif

* DEM/DSM(TIF)をGRASS GISにインポートする方法：
単純に、下記でインポートできる。
`r.import input=/Users/d-wei/Downloads/gm-jpn-el_u_1_1/jpn/el11.tif output=japan_el -o --overwrite`
![Screen Shot 2020-09-01 at 14.36.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/1a644771-e077-7797-a3cd-12142b898e87.png)

* DEM/DSM(TIF)から作られた高さ含まれたデータがrasterにあるかどうかを確認：
`g.list raster | grep japan`
![Screen Shot 2020-09-01 at 14.47.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/07798ce6-db55-e2f3-e94e-15ce54a354b3.png)  
`r.info japan_el`
![Screen Shot 2020-09-01 at 16.41.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/89af2e72-6776-523a-4160-60237b3d2549.png)


* ハマったこと：regionの設定
インポートしたrasterが、今の`calculate region`外に存在する場合、そのrasterに対するクエリや計算ができなので、インポートしたrasterに合わせてregionを設定し直す必要がある。
`g.region raster=japan_el`

* モニターを一個起動する：
`d.mon wx0`
<img width="1019" alt="Screen Shot 2020-06-05 at 17.06.04.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/e7478dea-2190-0b3c-534d-4d3bca9b8410.png">

* japan_elをモニターに表示する：
`d.rast japan_el`
![Screen Shot 2020-09-01 at 14.50.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/05771184-e376-195e-0ecb-a3e55b44a16b.png)
表示しない場合、ズームの問題の可能性があるので、monitorの下記ボタンをクリックすること。
![Screen Shot 2020-09-01 at 16.55.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/830cd176-f8fd-188a-e1a6-febdfd5698f2.png)

* 東京湾の高さを確認する：
`r.what map=japan_el coordinates=139.916382,35.561278`
見事に0.
139.916382| 35.561278||0
* 富士山の高さを確認する：
`r.what map=japan_el coordinates=138.729172,35.361616`
結果：あれ、低い。。。
138.729172|35.361616||154
大丈夫、標高3,776までに補正してあげる。
`r.mapcalc.simple --overwrite expression=A*3776/154 a=japan_el output=japan_el_full`
それから確認すると
`r.what map=japan_el_full coordinates=138.729172,35.361616`
もちろん、正しくなっている。
138.729172|35.361616||3776


* 下記ツールを使って、富士山周りを適当にgeojsonとして作る、vectorとしてインポート
https://geoman.io/geojson-editor
![Screen Shot 2020-09-01 at 14.55.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/dac0eea6-4a4a-c2b7-6f8c-9306262e31af.png)
ダウンロードしたgeoman.jsonをvectorとしてインポート
`v.import input=/Users/d-wei/Downloads/geoman.json output=mount_fj -o --overwrite`
完了後には確認できる。
`v.info mount_fj`
![Screen Shot 2020-09-01 at 16.42.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/492d6ef8-31d8-8ce0-f7d4-b1b56c05efdc.png)
説明：一つのポリゴンで囲んだものなので、boundariesが1, areasが１、
　　centroidsはバウンダリーの内部か外部かを決める物で、もちろん富士山エリアの内部である。

* mount_fjをモニターに表示する：
`d.mon wx1`
`d.vect mount_fj`
![Screen Shot 2020-09-01 at 16.57.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/accc1567-d900-b816-dc33-1ef341976dc0.png)

* 色を変えたり、レイアを重ねて表示させると
![Screen Shot 2020-09-01 at 19.25.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/9092cfa9-56c2-ed77-4eda-19a1fb84b4c6.png)

* maskを作る
maskを作ると、マスク区域以外の部分は計算しないようになるので、今計算したいvectorであるmount_fjに対してmaskを作る。
`r.mask vector=mount_fj`
GUIを使って見ると、切り取られた様子は下記通り
![Screen Shot 2020-09-01 at 19.55.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/7f10e582-e91b-114c-60db-e2472164decc.png)


* そのまま海抜から体積を計算する場合：
`r.volume input=japan_el_full`
![Screen Shot 2020-09-01 at 19.57.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547017/87c19700-ff60-615c-9b6b-a44966e8294b.png)
最終結果が40.03（この座標系は適当な度数なので）
もちろん、海抜から上の部分を体積とするケースは、実運用中には少ないのだが、今回はあくまでも例として。

* 実際測量する場合：
下記の記事の通りに、経度緯度の座標系から、メーター単位の座標系に変換する必要が最初にある。
https://qiita.com/ericofjp/items/8800eef383819c96ef53

* ここからは適当な話
富士山エリア、１経度*１緯度の面積は1m^2ではなく、大体100km*100km=10,000,000,000m^2
なので最終的には40.03*10,000,000,000m^2=400,000,000,000m^2.
400km^3

検証としては、測量したものが円錐と考えよう。
その半径は、v.info mount_fj表示したProjectionからおおよそ推測すると、
(（35.41843-35.31008）/2 + (138.81397 -138.65948)/2)/2 = 0.06571
高さはもちろん3776m
なので体積は3.14*0.06571^2*3776/3=17。
170km^3
まあ結果の40.03と同じ桁なので、OKとしよう。

[googleで富士山の体積を検索](https://www.google.com/search?ei=8CtOX4StIJDr-QakhobwAw&q=%E5%AF%8C%E5%A3%AB%E5%B1%B1%E3%81%AE%E4%BD%93%E7%A9%8D&oq=%E5%AF%8C%E5%A3%AB%E5%B1%B1%E3%81%AE%E4%BD%93%E7%A9%8D&gs_lcp=CgZwc3ktYWIQAzICCAAyAggAOgcIABCxAxAEOgoIABCxAxCDARAEOggIABCxAxCxAzoGCAAQBBADOgQIABAEOgQIABBDOgcIABCxAxBDOgkIABBDEEYQ-wE6CggAELEDELEDEEM6DggAELEDEAQQJRBGEPsBOgYIABAEECVQt9EKWJrqCmCG8ApoAHAAeAGAAe4CiAH9F5IBCDAuMTQuMi4xmAEAoAEBqgEHZ3dzLXdpesABAQ&sclient=psy-ab&ved=0ahUKEwiE5dvv6MfrAhWQdd4KHSSDAT4Q4dUDCA0&uact=5)すると、大体そんなもんだと。

* todo: JGD2011(EPSG:6677 )と、WGS84/UTMzone5４N（EPSG:32653）の解説(https://epsg.io/ )
* tifファイルのメタデータ、gisのprojectionの解説（https://www.metadata2go.com/ ）
* 無料dataの入手先：landsat, SRTM(https://www.usgs.gov/ ), ASTER GDEM(https://ssl.jspacesystems.or.jp/ersdac/GDEM/J/ )

国土地理院(https://www.gsi.go.jp/top.html )