etopo1を利用してQGIS3で海底地形図を描く

前回 はＱＧＩＳで海岸線を描きましたが、

今回は海底地形図を描きます。

（１）海底地形データの入手

（２）ＱＧＩＳでラスタデータとして取り込み

の順です。

（１）海底地形データの入手

今回はＮＯＡＡのetopo1を使用します。etopo1は１分（だいたい1.8km）刻みの全地球の標高値のデータです。

以下のページから入手可能です。

https://www.ngdc.noaa.gov/mgg/global/global.html

データは４種類ありまして、

・南極と北極の標高値を氷の表面にするか、氷の下の地面にするか

・データの場所を、グリッドの真ん中にするか、セルの真ん中にするか

の組み合わせになります。下のグリッドかセルかについては

こちらに解説があります。

https://www.ngdc.noaa.gov/mgg/global/gridregistration.html

とりあえず今回は、氷の表面（ETOTP1 ice Surface) で、

セルの真ん中(cell registrated）を選択しておきます。

ＱＧＩＳで使用するので、georeferenced tiff形式をダウンロードします。

数１００メガありますのでご注意ください。

出てきたファイルを解凍して、どこかのフォルダにまとめて収納します。

（２）QGISでの描画

(2-1)「レイヤ」→「レイヤの追加」→「ラスタレイヤの追加」から

出てきた画面で、先ほどのetopo1（拡張子tiff)を選択します。

 

そうすると、画面左側のレイヤのところに、上で選択したetopo1のファイル名のレイヤが、

画面中に地形図が表示されます。

etopo1のファイルのプロパティで、測地系を確認してみます。

画面左側レイヤの表示されていところで、etopo1のファイル名のところを右クリックして

プロパティを表示させ、出てきた画面の左上の「情報」を選択すると測地系などが確認できます。

EPSG 4326のWGS84になっていると思います。

（２）海岸線を重ねて描いてみる

前回 描いた海岸線を重ね書きしてみます。

前回の海岸線データは「ベクタ形式」（今回のetopo1はラスタ形式）なので

前回同様「レイヤ」→「レイヤの追加」→「ベクタレイヤの追加」で

海岸線データを選択するとご覧のように海岸線が上書きされます。

出てこないよ！という場合は、レイヤの書き順がずれている可能性があります。

画面左側レイヤの欄内を確認して、海底地形のレイヤ(etopo1）よりも

下に海岸線のレイヤ（line)が表示されている場合は、海岸線を一番奥に描いてから

海底地形を上に重ね書するということが行われます。

この海岸線のレイヤをドラッグして海底地形のレイヤの上に置き換えると

海底地形を一番奥に書いてから、海岸線を描くという順になりますので

海底地形が表示されると思います。

海岸線のレイヤのプロパティで確認できますが、海岸線のファイルも

測地系がEPSG4326 WGS84となっており、etopo1と一致しているため

海岸線の場所も正しく表示されていると思います。

QGISで海岸線を描画する（GMTのpscoast的なこと）

  題名の通りなのですが、QGISで海岸線を描いてみます。

GMTでいうとpscoast的なことをします。

＊はじめに

(0-1)QGIS のインストール

まず qgis を web からダウンロードし、インストールします。

アドレスはこちらhttps://www.qgis.org/ja/site/forusers/download.html

長期安定版 (2020/6現在 Ver3.10) を選択してください。

(0-2)海岸線のデータについて

国内データでしたら、

たとえば国土交通省のＧＩＳサイトから

http://nlftp.mlit.go.jp/index.html

国土数値情報ダウンロードサービスへ入って

http://nlftp.mlit.go.jp/ksj/index.html

海岸線データのダウンロードコーナーへ入ると

http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C23.html

入手可能です。今回は全世界のデータを使うので別サイトから行きます。

＊海岸線を描こう

（１）海岸線データの入手（２）ＱＧＩＳで海岸線を描いてみる

の順に進めます。

（１）海岸線データの入手

全世界データの海岸線を探したところ、OpenStreetMapプロジェクトの派生で

海岸データを抽出してくれているサイトがありました。

https://osmdata.openstreetmap.de/

上のサイトからcoastlineコーナーへ入ると

https://osmdata.openstreetmap.de/data/coastlines.html

shapefile形式（ＱＧＩＳでいうところのベクタ形式）でダウンロードできます。

投影図法がWGS84かメルカトル図法かを選ぶことができますが、

WGS84が便利だと思います。QGIS内で「再投影」という機能を

使うと相互変換可能ですので、あまり気にしないで大丈夫です。

数100メガバイトありますので、ご注意ください。ZIP形式ですので

解凍して、すべてのファイルを一緒のフォルダに保存しておきます。

（２）ＱＧＩＳで海岸線を描いてみる

(2-1)保存したcoastlineデータをベクタとして取り込む

「レイヤ」→「レイヤの追加」→「ベクタレイヤの追加」と進みます。

 

出てきた画面の「ソース」で、先ほど保存した海岸線データを選択します

拡張子は.shpです。ダウンロード時にファイル名を変えていなかったらline.shpだと思います。

選択したら「追加」をクリック

（２－２）表示を見ます

画面に戻ると、こんな具合に海岸線が表示されていると思います。

（２－３）測地系を確認する

今回（１）でＷＧＳ８４をダウンロードしましたが、その測地系を確認してみます。

画面左側のレイヤーエリアの、現在表示しているレイヤの箇所で右クリックして

プロパティを選択

出てきた画面で「ＣＲＳ」のところに情報が載っています。

今回のファイルは、「EPSG:4326 - WGS 84 - 地理学 」となっています。

基本的はこのように、ダウンロードしたファイルセット内に測地系の情報が書き込まれていて

ＱＧＩＳに取り込んだ際に自動的に読み込まれます。

（２－４）投影方法を変えてみる（メルカトル図法を試してみる）

（１）のところで、メルカトル図法と相互変換できます、と記載しましたので

そちらの操作を試してみます。

画面左側レイヤーのところで、海岸線データを選択した状態で

「ベクタ」→「データ管理ツール」→「レイヤの再投影」と進みます。

出てきた画面で「変換先ＣＲＳ」の右側になるボタン（ＣＲＳの選択）をクリックします

、

出てきた画面の上側の「フィルター」欄に[3857]と、おまじないを打ち込むと

世界の座標参照系欄にMercator WGS84 Pseudo-Mercator EPSG:3857

がでてくると思いますので、そちらを選択して「OK」をクリック

戻った画面の「変換先ＣＲＳ」に先ほど選んだEPSG:3857 WGS84 Pseude-Mercator

が入っていることを確認してから「実行」をクリック

「実行」の上のバーが作業の進行とともに伸びて行って、１分ぐらいすると

完了します。完了すると画面左側レイヤーのところに

「再投影された」というレイヤーができています。

このレイヤーのところで、もともとあったレイヤー(ここではline）のチェックボックスを外して

「再投影された」レイヤーのチェックだけのコスト、メルカトル図法の図が出てきます。

少し上の、変換先のＣＲＳを変更することで、いろいろな投影方法にすることができます。

検索でつかったおまじないの数字はEPSGコードといわれるものです。

使っているうちに、いろいろ出てくるとおもいますが、

４３２６ぐらいがよく出てきます。詳しくは検索してく調べてください。

QGISで利用するベースマップの利用規約

ＱＧＩＳですが、立ち上げると真っ白な画面が表示されるところからスタートします。

立ち上げ後、通常「ベースマップ」と呼ばれる、下地になる図を表示させてから

作業開始すると便利です。このベースマップとして、

たとえばgoogle mapなどが使用できすが、

実際に使用する場合は、著作権に注意する必要があります。

使用許諾についてわかりやすいベースマップとして

（１）国土地理院地図と（２）Open Street Mapがあります。

（１）国土地理院地図

利用規約はこちら

https://www.gsi.go.jp/kikakuchousei/kikakuchousei40182.html

下記の「国土地理院の地図の利用手続」というページに詳しく書かれています。

https://www.gsi.go.jp/LAW/2930-index.html

こちらのページ内の

「1．国土地理院の地図の利用手続 （測量成果の複製・使用申請フロー）」

の箇所で、フローチャート式で使用方法の説明がありますので一読すると便利です。

（２）OpenStreetMap

ホームページはこちら

https://openstreetmap.jp/

利用規約についてはこちらが便利です

https://www.openstreetmap.org/copyright

XYZ形式の海底地形図を3Dプリンタでプリントする

４年ほどご無沙汰しておりました。
２年ほど前、研究室に３Dプリンタが来て、立体物がプリント出来るようになってきています。
今回は３Dプリンタで、海底地形図を３Dプリントしよう！、と色々ためした備忘録です。

元のデータは、１行毎に緯度、経度、水深が１行ったでーたが
何万行かつらなるｃｓｖデータ（いわゆるｘｙｚ形式）です。

プリントルートとして３つためしました。
（１）ｘｙｚ形式を点群データと見なして
-------------------------------------------------------
（１－１）MeshLabなど点群処理ソフトで海底地形の薄皮一枚（いわゆるサーフェス）のデータをつくる
（１－２）Fusion360など3D cadで薄皮に厚みをつけて、実体化（いわゆるソリッド）のデータをつくる
（１－３）ｓｔｌ形式で保存して３Dプリンタで印刷
-------------------------------------------------------

（２）XYZ形式をGMTで等高線化して
-------------------------------------------------------
（２－１）Fusion360などの3D CADで押し出して、実体化（いわゆるソリッド）のデータをつくる
（２－２）STL形式で保存して３Dプリンタで印刷
-------------------------------------------------------

の２つをまず考えたのですが、どちらもソリッド化がうまくできませんでした。
（２－２）は気合いで通せますが、等高線一つに１押出し、という操作がシンドイです。

（３）FOSS4G系のソフトであるをQGISを使用して
　XYZ形式をDEM形式であると見なして
-------------------------------------------------------
（３－１）点群データを読み込んでレイヤを作成し
（３－２）QGISの機能でレイヤ内のデータをデータ補間して、３Dっぽい表示にする
（３－３）QGISのアドイン（DEMto3D）でｓｔｌ形式で出力
-------------------------------------------------------
としたら一応プリントできました。

（３－２）’からの別ルートとして
QGISのアドインのqgis2threejsを使用するとSTL形式まで出力できる
という情報もあったのですが、私のQGISのバージョン(3.4)では
qgis2threejsビューアーでの表示はできても、STL形式で出力する方法がみつかりませんでした。

以下（３）の手法を記載します。
-------------------------------------------------------
ｘｙｚについて、ｘ、ｙは小数点形式の緯度経度
（133.423423度、32.32423度みたいな形）
ｚは水深のｍからスタートしました。

(3-0)まずqgisをwebからダウンロードし、インストールします。アドレスはこちら
https://www.qgis.org/ja/site/forusers/download.html
今回長期安定版の3.4を選択しました

(3-1)早速qgisを立ち上げ、
１）「プロジェクト」→「新規作成」で新しいプロジェクトをつくって
２）「レイヤ」→「レイヤの追加」→「デリミティッドテキストレイヤ」の追加を選択し出てきた画面で
・[フィアイル名]のところで、XYZ形式の元データを選択
・[レイヤ名]を適当につける。何もしないとファイルの張子名になるのかな？
・[ファイル形式]は適切なものを選択（今回はCSV形式でやりました）
・[レコードとフィールドのオプション]で無視するヘッダの行数を選択（今回は０でした）
・[ジオメトリ定義]でｘフィールド（東西経度の列）、Yフィールド（南北緯度の列）を選択
　[ジオメトリのCRS]でｘｙｚデータの測地系を選ぶ（今回はWGS84　EPSG:4236でした）
・[レイヤ設定]はデフォルトのまま
で「追加」をクリックすると
画面中央に、ｘｙｚファイル内の点がｘｙプロットされると思います（深度情報は無し）

(3-2)データ補間します。
1）画面左側「レイヤ」エリア内に２）で指定したレイヤ名が出ていると思いますが
　その左側のチェックボックスにチェックが入った状態で、レイヤ名をクリックすると
　チェックボックスごと反転し選択状態になりますので、その状態で
　「プロセッシング」→「ツールボックス」と選択すると画面右側にツールボックスがでてきます
2）ツールボックスから「データ補間」の右向き矢印をクリックすると下向き矢印に代わり
　「IDW補間」「TIN補間」「ヒートマップ・・・」がでてきます。、
　補間ですがIDWが調子良さそうです
3）「IDW補間」をダブルクリックすると、新しいウインドウがでています。パラメータタブで
・[ベクタレイヤ]を(3-1)でつくったレイヤを指定
・[補間する属性]で深さ方向の列を選択
・１行下の枠の右上に出ている＋マークをクリックすると、
　ベクタレイヤ、属性、タイプというウインドウの中に、上で指定したレイヤと深さ方向の列がでてきます。
　その状態でタイプは「ポイント」を選択
・[距離係数]はあとで出てきた図で適当に調整しますが、３～１０ぐらいかなぁ、
・領域は右側の...の場所をクリックすると、「キャンパス領域を使用」「キャンパス上で領域を指定する」「レイヤの領域を使う」
　が選べますので、プリントしたい箇所を適当にえらんでください。
　ちなみに「キャンパス上で領域を指定する」を選択すると、2)でプロットされた画面にいけますので
　ここでプリントしたい場所をマウスで指定して話すと、領域に数値が入ります。
・[出力クラスタサイズ]で地形データの出力する細かさを指定します。
　例えば行に数字をいれると、上で指定した領域の面積に合わせて、自動的に列が決ります。
以上設定が終わりましたので、「実行」をクリックすると、地形図が計算されて、画面中央に（おそらく白黒で）表示
されるとともに、画面左側のレイヤに「補完された」というレイヤがでてきます。
画面中央の図がきれいにかけるまで、適当にパラメータを調整してください。

（３－３）最後にSTL形式の出力にかかります。
1)STLの出力にはプラグインのDEMto3Dを使用しますのでインストールします
・”インターネットに接続した状態で”
　「プラグイン」→「プラグインの管理とインストール」を選択すると
　プラグインがずらっと出てきます。
・「すべてのプラグイン」をクリックすると名前順に一覧がでますのでその中からDEMto3Dにチェックをいれて
　「プラグインをインストール」をクリックします。これでプラグインがインストールできました
2)DEMto3Dを使ってSTL形式で出力させます。
・「ラスタ」→「DEMto3D」→「DEMto3Dprinting」をクリック
・[Layer to print]で５）にて描いた画像のレイヤを選択
・[print extent]で○っぽいマークをクリックすると、上のレイヤの緯度経度が自動で入ります
・[model size]でspcingは推奨の0.2を
　sizeはプリントしたい大きさをmmで入力します
　scaleはそれにあわせて自動で1/1000とか表示してくれます
・[exaggeationterrain]でｚ軸の深さ方向の拡大量を入力、どうも１０が最大のようです。。。
・[height base]で地形データをプリントするときのゲタを何m分はかせるか入力（0でもいいです）
以上で「EXPORT　to STL」をクリックして、
成功すると、STLをはき出すファイル名が聞かれますので「保存」をクリック
失敗すると、「attention fill the data correctly」が出ますので、きちんとデータを入力したり
レイヤデータを正しく選択するなどの処理をして、成功できるよう頑張りましょう。

このデータを３Dプリンタのスライサーにもっていけばプリントできますが
出てきた海底地形データは、深さ方向が逆向き、深さ方向の拡大が中途半端
な状態だと思います。きれいにプリントするためには
ｘｙｚのデータに戻って、
・深さデータを、どこかの海底面からの高さデータに変換
・高さデータを数倍程度拡大しておく
の処理をした後で、描画するときれいにできると思います。

今回高知県沖合の海底データを、ダビンチの2.0 Duoで印刷したらこんな感じになりました。
また時間があれば画像付で書きます。