概要

SENSYN Robotics(センシンロボティクス)の衛です。

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

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

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

使ったのはGRASS GISである。
GRASS GISインストールする必要がある。

測量原理

一枚写真からでは、被写体の立体的な表現（情報記録）ができませんが、
ドローンが一定的なオーバラップ率で撮った、同じ地面の場所の複数写真を、
立体視あるいはステレオマッチングの原理使って変換をかけると、その場所の地面高度が含まれたデータが作成できます。そのうちの一つはDSMデータです。

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

富士山を例としする。下記の図の０がある暗いセルは地面（あるいは海抜0m）としする。1~4のセルは高さが徐々に上昇する富士山の部分を示している。

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

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

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

GRASS GIS セットアップ

GRASS GISのデータ構造は、location＝＞mapset＝＞各種rasterとvectorデータとなる。それを見るmonitorや各種ツールもある。
基本GRASS GISをダウンロードしたら、インストールしてから開いたら、
ウィザードが開かれる。

「2. Select GRASS Location」では、Newを選んで、下記のように入力。

「Create a generic Cartesian system (xy)」を選んで次々とする。

「Start GRASS session」をクリックする。

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
DEM/DSM(TIF)から作られた高さ含まれたデータがrasterにあるかどうかを確認：
g.list raster | grep japan

r.info japan_el
ハマったこと：regionの設定
インポートしたrasterが、今のcalculate region外に存在する場合、そのrasterに対するクエリや計算ができなので、インポートしたrasterに合わせてregionを設定し直す必要がある。
g.region raster=japan_el
モニターを一個起動する：
d.mon wx0

japan_elをモニターに表示する：
d.rast japan_el

表示しない場合、ズームの問題の可能性があるので、monitorの下記ボタンをクリックすること。
東京湾の高さを確認する：
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

ダウンロードしたgeoman.jsonをvectorとしてインポート
v.import input=/Users/d-wei/Downloads/geoman.json output=mount_fj -o --overwrite
完了後には確認できる。
v.info mount_fj

説明：一つのポリゴンで囲んだものなので、boundariesが1, areasが１、
　　centroidsはバウンダリーの内部か外部かを決める物で、もちろん富士山エリアの内部である。
mount_fjをモニターに表示する：
d.mon wx1
d.vect mount_fj
色を変えたり、レイアを重ねて表示させると
maskを作る
maskを作ると、マスク区域以外の部分は計算しないようになるので、今計算したいvectorであるmount_fjに対してmaskを作る。
r.mask vector=mount_fj
GUIを使って見ると、切り取られた様子は下記通り
そのまま海抜から体積を計算する場合：
r.volume input=japan_el_full

最終結果が40.03（この座標系は適当な度数なので）
もちろん、海抜から上の部分を体積とするケースは、実運用中には少ないのだが、今回はあくまでも例として。
実際測量する場合：
下記の記事の通りに、経度緯度の座標系から、メーター単位の座標系に変換する必要が最初にある。
https://qiita.com/ericofjp/items/8800eef383819c96ef53
ここからは適当な話
富士山エリア、１経度１緯度の面積は1m^2ではなく、大体100km100km=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.140.06571^23776/3=17。
170km^3
まあ結果の40.03と同じ桁なので、OKとしよう。

googleで富士山の体積を検索すると、大体そんなもんだと。

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 )