データ無しの水面の標高を基盤地図情報数値標高モデルのGeotiffより推察して埋めたい

[arai...@gmail.com]

基盤地図情報基本項目「水涯線」を取り入れて考えていましたが、自力で答えを出せませんでした。

目的は海、河川、池などがデータ無しの場合が有るので海抜0として扱うのは「3D」的に無理が有るためです。

当初、水涯線をラスター化すれば何とかなるのかと考えましたが、そのラスター化すら出来ずエラーを見るだけでどうしようもありませんでした。

調べていて、そもそも水涯線に高さ情報なんて乗っているはずが無いから、簡単に標高を埋めることは難しいのではないか？と考えるに至り、この質問となりました。

おそらく、水涯線の形状を一定のグリッド等でサンプルするポイントを決め、ラスタレイヤの近く標高を参照しながら埋めていくことになるのか？と考えたのですが、今の知識では

どうしようもなく、またしてもお世話になるしかなく。。。。

さて、表題のことをするのは可能なんでしょうか？

できるなら、どんな手続きをとるのでしょうか？

また、ラスター化のときのエラーをどうすればよいのか、ご教示ください。

情けないのですが、よろしくお願いいたします。

[福岡]

こんばんは
突破口になればというレベルの投稿ですが、よろしくお願いいたします。

「no data」に関しては周囲のDEMからデータ補完してみてはどうでしょうか？
「ラスタ － 解析 －欠損値の補完」で付近のDEMからデータを補完することができます。

ただ、これだと問題が１つあって、内水面に関しては良いと思うのですが、海岸付近の海面も付近のDEMからデータ補完されてしまいます。
これに対処する方法としては、陸地部分のマスクレイヤを作って「ラスタ － 抽出 －クリッパー」でくり抜く方法があると思います。
これで海の部分が「no data」でくり抜かれます。

最後に、私だったら海の標高は「0m」にしたいと思いますので、データの入れ替えを行います。
データの入れ替えはプロセッシングツールの「SAGA － Raster tools － Reclassify values」 で「replace no data values」にチェックを入れて値を「0」にするとできると思います。

[福岡]

ステップが多く煩わしい手順ですが、それなりのDEMを得ることができました。
もっと良い手順があるかもしれませんので、参考程度にしてください。

昨日投稿してから、あの手順でうまく行くか試してみることとしました。
海と欠損値が入っているようなところを選んだところ、海岸線が複雑なのは良いとしても島などがあり、Shapeファイルをフィルタとして切り抜くのには無理がありそうだと感じました。

ここで方針転換することとしました。海と陸地を分類するラスタファイルを作成して、ラスタ計算機で処理することを目標としました。
まず、上のラスタをラスタ計算機で処理します。欠損値として－9999が指定しているとうまく動作しないので、プロパティの透過性タブにある「データなしとする値」のチェックを外します。
外すと次のように表示されます。

ラスタ計算機で処理します。ラスタ計算機で使っている「!=」は不等号です。
（　）内の条件が真なら「1」、偽なら「0」を返します。式は－9999（海など）の値でないのなら1×9999ということで9999。
海などの場合は0×－9999という計算をします。

得られたラスタです。

このラスタをShapeファイルに変換します。使うツールはプロセッシングツールにあるGRASSのr.to.vectです。

得られたShapeを海の部分だけ残して、後は削除したものです。

今度はShapeファイルをラスタに戻してやります。使うツールはプロセッシングツールにあるGRASSのv.to.rast.attributeです。
セルサイズのところに入っている数字は、他のラスタのメタデータを確認して揃えてあります。

得られたラスタです。陸となっている部分は「no data」となっているので、もう一手間掛けます。

プロセッシングツールにある SAGA － Raster Tools － Reclassify values を使います。
上のラスタを対象として、下の方にある「replace no data values」を「1」、「replace other values」を「0」にして実行します。

得られたラスタです。海が「1」、陸地が「0」で分けることができました。

欠損値のデータ補間を行います。最初に示したラスタへの処理です。処理前に、透過性タブにある「データなしとする値」にチェックを入れておきます。
ラスタ － 解析 － 欠損値の補間　を使って得たラスタです。昨日に投稿で書いたように海岸線付近が余計に補正されています。

このラスタを陸と海を2値で分けたラスタを使って修正します。

最終的に得られたラスタです。海の部分の値が「0」となって、陸地部分に欠損値はありません。

長々と書きましたが、一部でも参考になる手順があれば幸いです。

[arai...@gmail.com]

福岡様

いつもありがとうございます。

とりあえずは最初のアドバイスを実行しているのですが、欠損値の補完はできてると思います。

ですが、クリッパーがどうしてもエラーでクリアできず困っていたのですが、tiff保存を上書きから別名保存ですとクリア。

QGISでは上書き保存はご法度でしょうか？以前の質問でクリッパーのクリアができてなかったので放置していました。

ちゃんと整理してみます。

やった事

1：行政区域をDLし、shpに変換、QGISに取り込み、マスクとして利用

2：細かく分かれたポリゴンを選択、編集モードから選択地物の結合、CRSの統一、クリッパーで抽出＜-今ここです。

で、直前のアドバイスは今の私がすんなり理解するには高度すぎるので、クリッパーが整理できれば最初のアドバイスの方針を試そうと思います。

「SAGA － Raster tools － Reclassify values」は試せてませんのでどうなるか分かりませんが、コツコツと勉強してみてから高度なプロセスも

試してみようと思います。最初の画像は取り込んだばかりの生のGeotiffでよろしいですよね？

「ラスタ － 解析 －欠損値の補完」で検索距離だけ50にしてやってるのですが、デフォルトの100とかでは結果は埋まる距離が変わるだけと理解してよいのでしょうか？

エイヤで埋めるだけ埋めて海だけ0を入れたら済むのかなー？と考えたりしています。河口付近が気持ち悪いのかもしれませんが。

[福岡]

こんにちは

ラスタをクリッパーで抽出した後のデータがどのようになるのかを資料で確認してみました。
ラスタにさわり始めた頃、参考にしていた資料です。

QGIS はじめてのラスタ解析
https://www.slideshare.net/Mayumit/qgis-38312837

24～26ページにクリッパーの取り扱いが掲載されています。
そして26ページには「データの無い値」の設定で、どのような出力となるかについて解説されています。
今回の場合、一番下の「データがない場合：チェックなし」とすると良いように思います。

検索距離ですが、検索範囲が拡がると当然欠損値が埋まる範囲が広くなると思います。
ただ、それだけではなく、補完値への影響があるものと想像していました。
イメージとしては、河川で欠損値があってこれを埋める場合、検索範囲が狭いと河道付近の河原の値に近くなり、
検索範囲が広くなると土手を越えた範囲も計算に含まれるというものです。
実際の仕事などで欠損値の補完は使ったことが無いので、値の検証はしていないのですが。。。

なお最初の画像は、基盤地図情報からダウンロードし、Geotiffに変換したものからクリッパーで切り抜いたものです。
欠損値があるのは航空レーザー測量の5mメッシュと思ってダウンロードしましたが、未整備の区画があって、そのままだと使えなかったからです。

最後に、既にご存じかもしれませんが、プラグインの「Value Tool」をインストールしておくことをお勧めします。
このツールはラスターの値を確認するのに便利で、マウスポインタ位置で値を確認することができます。

[arai...@gmail.com]

福岡様

ありがとうございます。

クリッパー

やはり、私の環境（仕事場）ではGeotiff上書き保存ではエラーですね。誰がやっても同じで、できないなら保存時に許されないアラートを出すべきだと思いました。

他のソフトではありえないですからね。書き換えるぞ！と言いながらNGは。。。

検索距離の大小

海岸付近、河口等で両側から合わさる部分が微妙な幅が有る箇所では、補間値が変わるようです。PhotoShopの差の絶対値で確認できました。

両方向の当たりが完全にオーバーラップするところは一致してると思います。今のところ補間値がシビアでなければいけないことがありませんので、深く考えずに

行けそうです。

SAGA － Raster tools － Reclassify values

どこにあるのか迷いましたが、海部分を0で埋めることができたようです。

最終のPNGまで試していませんが、最初の手法で何とかなるように思いました。

ラスタ計算機手法

プロセスの考え方はなんとなく理解できそうですが、やはりすらすらと思い浮かばない初心者にはきついです。。。はい。

最後の合体プロセスがなんか失敗している状況ですが、

やっていての不明点

1：河川が有り、途中で切れるケースで海のポリゴン部分だけ残したが、ラスタ化データが河川の切れた部分までしか生成されなかった。（エリアが小さい）

2：ラスタ化の際に海部分に1セル？のゴミで囲まれる？あとあと問題にならない？

参考資料

一度見ていました。が、その時はこんな事をするとは思わず、真剣さが無い読み方でした。もう一度読み直します。

必要に駆られていないと関係なさそうな箇所は覚えてないもんですね。覚えるところは覚えたとは言えませんが。。。

Value Tool

レイヤごとに合わさった確認ができて便利そうですね。それをどう使うと助かるのか、まだ経験がありません。。。

高さはどこかに表示されないかなぁ。とは思っていたのでそれは勉強できました。

[福岡]

こんばんは

投稿が遅れて申し訳ありません。
ラスターを使った処理がうまくいかない原因について考えていました。

ポリゴンからラスタへの変換で得たGeoTiffに不具合があるということですので、１つだけ確認をお願いいたします。
この変換で得られるGeoTiffは、他のGeotiffと重ねて使用しますので、縦横のピクセル数が一致していなければなりません。
もし違っていましたら、他のGeotiffのレイヤプロパティのメタデータタブの一番下にあるプロパティで「ピクセルの大きさ」を確認してメモしてください。

また、私の試したデータでは北と西が海になっていたので、「領域の範囲」を特に定める必要がなかったことに気がつきました。
もし、変換したGeoTiffを他のGeoTiffと重ねて狭い範囲だった場合には、変換の際に「領域の範囲」を指定して変換する必要があります。
他のGeotiffのレイヤプロパティのメタデータタブの一番下にあるプロパティで「レイヤ領域（レイヤのオリジナルプロジェクション)」を確認してメモしてください。
これらを「v.to.rast.attribute」でポリゴンをラスタに変換する際に「GRASS GIS 7 領域の範囲」と「GRASS GIS 7 領域のセルサイズ」にメモしておいた数値を入力してから、変換を実行してください。

縦横のピクセル数が一致していた場合には、申し訳ありませんがもう一度考えさせてください。m(__)m

[arai...@gmail.com]

福岡様

ありがとうございます。

「GRASS GIS 7 領域の範囲」

数値の順番やコピペの仕方で四苦八苦し、なんとか白黒2値のデータを得られましたが、なぜか最後のラスタ計算機では補間データと同じデータになりました。

何か間違っているかもしれませんが、上手く行ったとしても海から河川の上流まで海扱いになりますよね？

仮に大阪湾から琵琶湖まで続いてるデータになっているとしたら、3D的にそれは困ることになりますので、陸地近辺を埋めるだけ埋め、補間しておいて行政区域で抜き

海はデータ無しで置いておいても、最終のPNGでは黒になって、結果同じようになりますよね？

山深いところまで河川が上がっていた場合、3D的には海抜0まで河川が地割れ状態になります。これはまずいことなんです。適当に埋まっているのがベターです。

>Shapeファイルをフィルタとして切り抜くのには無理がありそうだと感じました。

簡単な最初の考え方ではどう無理が有るのでしょうか？

[福岡]

混乱させてしまって申し訳ありません。m(__)m

私が試してみたデータには島が多く、陸部分に関して細かいポリゴンが多数生じました。
そのShapeファイルをマスクレイヤに使うとエラーが発生するため、ラスターでの処理に切り替えたものです。
そういうことですので、Shapeファイルで切り抜けるのであれば何も問題ありません。
余計なお時間を取らせてしまって、申し訳ありませんでした。

[arai...@gmail.com]

あ、そういうことでしたか、分かりました。

ちなみにサンプルの場所はどこでしょうか？

一度簡単な手法を試してみたいと思います。

混乱と言うよりも何か違うのだろうと考えていたので、違う考え方の手法で試すのも後々の考え方の参考にもなっているので気にしないでください。

何か理解しようとしてやらないと頭には残りません。。。

[福岡]

おはようございます。

サンプルの場所は男鹿半島の北西部です。
生のデータが最初からあった方が良いと思いますので、後ほどメールの方でダウンロードするurlをお伝えします。
一連のGeoTiff、Shapeファイルを圧縮ファイルにしましたので、よろしくお願いいたします。

[arai...@gmail.com]

福岡様

データまでご用意いただきありがとうございます。

行政区域で抜く方法を試したのですが、なんかエラーで失敗してます。

手順が間違っているのか、上手く行った例と比べて見ようと思います。

島が多いからではないように感じた気もするので、じっくりやってみます。

まずはお礼申し上げます。

[福岡]

おはようございます。
もっと簡便な方法がないかを考えていたのですが、発想を変えると単純な方法があることに気がつきました。

データに欠損値があるのはレーザー測量で水面の標高を計測できないためで、海や内水面の一部に欠損値が発生していました。
また、国内全てに5mメッシュデータは整備されていません。
一方、10mメッシュは国内全てが整備されていて、内水面にも欠損値がありません。

そこで5mメッシュをベースにして、欠損値のある箇所について10mメッシュデータから補完すると良いのではないかと考えました。
これなら作業手順も単純です。

作業手順としては、同じ2次メッシュ範囲で5mメッシュと10mメッシュ（10B（地形図の等高線）のデータをダウンロードします。
エコリスさんの変換ツールでDEMに変換するのですが、5mメッシュでは海の標高を－9999とします。これで海と内水面の欠損値のある場所が－9999になります。
一方10mメッシュでは海の標高を0で変換します。これで海の標高が0mで内水面全てが標高を持っているDEMとなります。

これをラスタ計算機で処理します。DEMファイルの名称をわかりやすいように05mDEMと10mDEMとしました。
数式は次のとおりです。
( "05mDEM@1" != -9999 ) * "05mDEM@1" + ( "05mDEM@1" = -9999 ) * "10mDEM@1"

( "05mDEM@1" != -9999 ) は真のときに「1」、偽のときに「0」を返してきます。「!=」は不等号です。
よって「( "05mDEM@1" != -9999 ) * "05mDEM@1" 」は、5mDEMが欠損値で無いときは5mDEMの値、欠損値の場合は「0」となります。
そして 「 ( "05mDEM@1" = -9999 ) * "10mDEM@1"」は、5mDEMが欠損値であったときは10DEMの値、欠損値であった場合は「0」となります。

ラスタ計算機を使うのが1回だけですので、作業が楽だと思います。
なお5mDEMですが、レイヤプロパティの透過性タブで「データなしの値」のチェックを外しておいてください。
こうしないと、ラスタ計算機で値を-9999と認識しませんので、海や内水面が「no data」で出力されてしまいます。

[arai...@gmail.com]

福岡様

ありがとうございます。

これはドンピシャな発想ですね！まだ手が付けれないのでやってないのですが、私でも簡単そうなのが予想できます。

まあ、これまでラスタ計算機等を勉強したせいも有りますが。。。w

QGISと違うところ（3D）も光が見えてきてQGISとの合わせ技が実りそうな予感。。。

かまって頂いて、ありがたい限りです。m(__)m