1.はじめに
国土地理院では、「基盤地図情報サイト」で5mメッシュと10mメッシュの標高データを提供しています。
このうち5mメッシュは航空レーザー測量から作成されたものなので、レーザー測量が苦手とする湖沼などの水面部は欠測となっています。
一方、10mメッシュは地形図からの読み取りのため、このような欠測はありません。
このようなことから、できれば詳細なデータである5mメッシュを使いたいけれど、欠測箇所があるのでそのままでは使えない。そこで何とかして、欠測箇所を推定値または10mメッシュの値で補完したい、という意見を見かけることがあります。
私が出没する「QGIS初心者質問グループ」でも話題になったことがあって、10mメッシュから補完するのが容易だろうということで落ち着き、作業手順も投稿されています(当時の投稿はこちら)。
ただQGISの仕様が少し変わっていて、全く同じ手順ではうまく行かないようです。そこで作業手順を見直して、Qiitaへ投稿することにしました。
2.GeoTiffへの変換
国土地理院で提供しているデータはxml形式ですが、QGISで利用するためにはGeoTiffへ変換する必要があります。
主たる変換方法は2つ有り、1つはエコリスさんが提供している「基盤地図情報標高DEM変換ツール」を使う方法です。
このツールは10年以上前から公開されていますので、QGISユーザーのほとんどはお世話になっていると思います。関連ページのURLを記載します。
https://www.ecoris.co.jp/contents/demtool.html
もう1つは2021年からMIERUNEさんが提供している「QuickDEM4JP」というプラグインです。こちらは公式リボジトリーからインストールすることができます。
説明記事がQiita内にありますので、URLを記載します。
https://qiita.com/nokonoko_1203/items/b99aa733cb215305f8aa?fbclid=IwAR22EYldb8vqP90NLG1fi1OkfL0IfIQBt5pl4f1bzYmnAYbv3RNINfjiZM4
2つとも使いやすい優れたツールですが、若干の違いがありますので、次表に掲載します。
基盤地図情報 標高DEMデータ変換ツール |
QuickDEM4JP | |
---|---|---|
形 態 | スタンドアローン | QGISプラグイン |
出力CRS | 経度緯度、UTM、平面直角座標 (JGD2000又はJGD2011) |
QGISで扱うことのできるCRSを 任意に選択できる |
No Dataの扱い | 「0」または「-9999」 | 「-9999」のみ |
その他 | - | Terrain RGBの出力が可能 |
次からツールの使用方法について説明しますが、基盤地図情報からダウンロードしたファイルは、解凍して、変換元になるxmlファイルを1つのフォルダにまとめておく必要があります。
2-1.基盤地図情報 標高DEM変換ツール
起動は「変換結合.vbs」を実行します。DEM変換を良く行う方は、ショートカットを作っておくと良いと思います。
起動したら、メッセージが出ますので、そのまま「OK」ボタンを押します。
次の操作では、参照座標系を「緯度・経度」、「UTM」、「平面直角座標」のうちいずれかを選びます。傾斜などを計算する場合に単位が「m」である投影座標が都合良いこと、投影座標のうち「平面直角座標」の方が精度が高いことから「平面直角座標」を選ぶことが多いように思います。
投影座標を選んだ場合、「UTM」ではゾーン番号、「平面直角座標」では系番号を入力します。
陰影起伏図の作成については、必要であれば「はい」、不要の場合は「いいえ」を選択します。
xmlファイルをまとめておいたフォルダを選択し、「OK」ボタンを押します。
海域の標高の選択画面が出ます。海域という言葉を使っていますが、正確には「nodata」となっている場所にどのような値を付与するかということです。
よって「10mメッシュ」では言葉どおり「海域」となりますが、「5mメッシュ」には湖沼や河川で「nodata」となっている箇所もありますので留意してください。
ここで「はい→0」を選ぶと「nodata」と「標高0m」が同じ値になって見分けがつかなくなりますので、「いいえ→-9999」を推奨します。
「OK」ボタンを押して、処理を実行します。
処理が終了するとメッセージが出ますので「OK」ボタンを押します。
指定したフォルダ内に出力された「merge.tif」が、目的のDEMファイルです。
2-2.QuickDEM4JP
QuickDEM4JPは、QGISのメニューの「プラグイン - Quick_DEM_for_JP - Quick_DEM_for_JP」で呼び出します。
「入力設定」の「形式」では「'xml'または'xml'を含む'zip'」または「'xml'を含むフォルダ」を選ぶようになっていますが、この作業例ではxmlファイルをフォルダにまとめていますので「'xml'を含むフォルダ」とし、次の「DEM」の設定でxmlファイルをまとめたフォルダを指定します。
「出力」の「形式」では「GeoTiff」にチェックを入れ、次の「出力先」でパスを含めたファイル名を指定し、「CRS」で参照座標系を指定します。
なお「出力」の「形式」では「Terrain RGB」にもチェックを入れて出力することもできますが、入門者にはあまり必要無いものだと思います。QGISを使い慣れてから、必要に応じて出力するというスタンスで良いと思います。
処理が終了するとメッセージが出ますので「OK」ボタンを押します。
3.補完作業
3-1.処理前の準備
変換して得たGeoTiffファイルですが、作業説明をわかりやすくするため、それぞれ05m_DEM.tif、10m_DEM.tifと名前を変更して作業を行います。
処理前の準備としては、nodataを値で埋めるという作業を行います。これはnodataとして設定されているセルがラスタ計算機の計算対象から外れてしまうため、これを回避するため処理です。
以前はレイヤプロパティの透過性タブでnodata値に入っているチェックを外すだけで回避できたのですが、現在はチェックを外しても回避できない仕様になっています。
nosataを値で埋めるためのツールは、プロセッシングツールの「ラスタツール - nodataセルを値で埋める」です。
埋める値を-9999として実行するとnodataとして扱われなくなるので、ラスタ計算機で扱えるようになります。
これを5mメッシュ、10mメッシュの両方で行い、出力ファイル名の後ろに「_filled」を付けて出力します。
3-2.ラスタ計算機でデータ補完
5mDEMで-9999の値を持つセルに、10mDEMの値で補完します。
ラスタ計算機では、等号(=)や不等号(!=)で記述した場合、真の場合に「1」、偽の場合に「0」を返します。それを利用して次のような式にしました。
("05m_DEM_filled@1" != -9999) * "05m_DEM_filled@1"+("05m_DEM_filled@1" = -9999) * "10m_DEM_filled@1"
出力ファイルは05m_DEM_mix.tifとし、空間範囲は「選択レイヤの領域を使用」で05m_DEM_filledと一致させています。ラスタ計算機の設定は次のとおりです。
3-3.nodataを値で埋める
出力のレイヤプロパティをみると、nodata値が設定されています。
この部分も含めて最後の処理でnodataセルを設定したいので、プロセッシングツールの「ラスタツール - nodataセルを値で埋める」を使って-9999で埋めて、nodataとしたいセルの値を揃えます。なお、出力は一時レイヤで差し支えないので「一時ファイルに保存」という設定にしています。このファイルを保存したい場合は、適宜、出力先を設定してください。
3-4.便宜的にClip機能を使って、nodata値を設定する
ラスタの値からnodataを指定するツールがあれば良いのですが、残念ながら見つけることができませんでした。そこでClip機能にnodata値を設定するオプションがあるので、それを便宜的に利用してnodata値を設定することにします。
さきほどの出力はレイヤパネルに「出力ラスタ」と表示されます。このレイヤを対象にメニューの「ラスタ - 抽出 -範囲を指定して切り抜く」を使います。
ポイントとしては「切り抜き範囲」を右端にある逆三角をクリックして展開し「レイヤから計算 -出力ラスタ」で設定することと、「出力バンドに指定nodata値を割り当てる[オプション]」を「-9999」に設定することです。
また、この出力は最終目標ですので、「出力ファイルの設定」はファイルとしてきちんと保存する設定としてください。
これで作業は終了となります。