ラスタデータの切り抜き及びその部分の補間方法について

[ねこ]

お世話になっております。

主に2つ質問があるのですが、既出でしたら申し訳ありません。

使用しているQGISのバージョンは3.12.2、GRASSGISのバージョンは7.8.2です。

ある地域一帯のラスタレイヤをレイヤＡ（基盤地図情報ダウンロードサービスの5ｍメッシュセット1つ分位だと考えてください）、

小さい成層火山1つの範囲を示したポリゴンレイヤをレイヤＢとしたとき、

レイヤＡからレイヤＢの範囲を取り除いたレイヤＣを作成するにはどうすればよいのでしょうか。

ベクタレイヤ同士ならば差分（difference）を使えばよいと考えられますが、ラスタレイヤはマスクレイヤによる切り抜きしか無かったので…

また、作成したレイヤＣにできたレイヤＢの範囲の穴をTINなりIDWで補間するにはどのツールを使えばよいのでしょうか。

もしその山が該当地域になかった場合の推定を行いたいと考えております。

もし方法をご存じの方がいらっしゃればご教示いただきたく、どうかよろしくお願いいたします。

国土数値情報ダウンロードサービスが1週間ほど前に大幅にリニューアルされて戸惑っている日々です。

[福岡]

ねこ　様

こんにちは

実際に手を動かしていないのでイメージだけで申し訳ありませんが、ラスタ計算機でDEMを合成するという方法はいかがでしょうか？

以前、5mメッシュに内水面の標高が含まれていないのを10mメッシュの値で補間する方法について投稿したことがあります。

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

https://groups.google.com/forum/#!topic/qgisshitumon01/mTnnGJGk0-U

最後の方で、ラスタ計算機を使う方法について投稿しています。

この方法を使えば、TINから生成したDEMをマスクレイヤでくり抜き、基盤地図情報からダウンロードしたDEMと合成するというように応用できそうに思います。

ご要望に添う内容であれば、ご覧くださいますようお願いいたします。

[ねこ]

福岡　さま

お礼が遅くなり申し訳ありません。

返信していただきありがとうございます。

まずTINから生成したDEMをマスクレイヤでくり抜くことですが、

ラスタ計算機で値を-9999にしたラスタレイヤを作成することはでき、

そこから5ｍメッシュデータとの組み合わせ方の理解が至らなかったため

最終的にラスタの結合を利用してくり抜いたような形にできました。

このラスタでくり抜き値が-9999となったエリアを、基盤地図情報からダウンロードしたDEMと合成するのではなく、

周りの5ｍメッシュの地形データから補間したダミーデータを作り埋めるにはどのようにすればよいのでしょうか。

難しい質問ではありますが、ご教示いただけますと幸いです。

[ありた]

こんにちは

（手元でテストできてませんが。。。）

ひとつめの目的については解決されたようですが、他の方法として、まず対象のラスタAの

サイズ以上の大きさのベクタを作成し、そのベクタと除外したいベクタBの差分を取ります。

つまりベクタBを反転させたベクタを作り、そのレイヤでマスクする方法もあるかと思います。

ふたつめですが、説明文を見る感じだと GRASS の r.fill.stats や r.fillnulls が利用できそうです。

[ねこ]

ありた様

お礼が遅くなり申し訳ありません。

福岡様の手法と共に参考にさせていただきました。

無事に切り取り方法については理解して慣れることができました。

しかし二つ目の目的である、作成したレイヤＣにできたレイヤＢの範囲の穴をTINなりIDWで補間する手段について、

GRASS の r.fill.stats や r.fillnullsを利用してみたのですが、エラーで表示されない・全く変わらないのどちらかでレイヤの穴を埋めることはできませんでした。

やはり5ｍメッシュのデータで1㎞四方の穴を埋めるのは無理があるのでしょうか。

度々の質問で申し訳ありませんが、ご教示いただけますと幸いです。

[ありた]

こんにちは

エラーが発生した場合は、エラーメッセージを全文転記してください。

「全く変わらない」というのはどういう意味でしょうか。入力ラスタと同じ内容の

ラスタデータが出力されるということでしょうか。

また、エラーが発生する場合と全く変わらない場合で再現性（この条件ならば。。。）が

あればちゃんと書いてください。

あとは nodata に対し補間するものですが、 r.null あたりで nodata の設定はちゃんと

できていますか？　（レイヤのプロパティ設定からでも可能かも）

[ねこ]

ありた様

 

説明不足で大変申し訳ありませんでした。

nodata の設定ですが、0のデータを追加のnodata値としてレイヤのプロパティで設定しました。

エラーが発生する場合と全く変わらない場合の件ですが、r.fillnullsを利用した際にはエラーが発生し、r.fill.statsを利用した際には変わらないレイヤが表示される状況になっています。

全く変わらないという件ですが、非常に似ているレイヤが出力されたというほうが正しかったのかもしれません。

まずGRASS の r.fill.stats と r.fillnullsを利用する際の画面です。値はすべてデフォルト状態です。

r.fill.statsを実行した際のログを添付いたします。

QGIS version: 3.12.2-București

QGIS code revision: 8a1fb33634

Qt version: 5.11.2

GDAL version: 3.0.4

GEOS version: 3.8.1-CAPI-1.13.3

PROJ version: Rel. 6.3.1, February 10th, 2020

プロセシングアルゴリズム...

アルゴリズム 'r.fill.stats' を開始しています...

入力パラメータ:

{ '-k' : False, '-m' : False, 'GRASS_RASTER_FORMAT_META' : '', 'GRASS_RASTER_FORMAT_OPT' : '', 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : '38140.64,47623.1745,-138036.8719,-129837.2135 [EPSG:6676]', 'cells' : 8, 'distance' : 3, 'input' : 'C:/Users/AAA/Desktop/a/nuki.tif', 'maximum' : None, 'minimum' : None, 'mode' : 0, 'output' : 'TEMPORARY_OUTPUT', 'power' : 2, 'uncertainty' : 'TEMPORARY_OUTPUT' }

g.proj -c proj4="+proj=tmerc +lat_0=36 +lon_0=138.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

r.in.gdal input="C:\Users\AAA\Desktop\a\nuki.tif" band=1 output="rast_5ed704bfe07874" --overwrite -o

g.region n=-129837.2135 s=-138036.8719 e=47623.1745 w=38140.64 res=26.56172128851541

r.fill.stats input=rast_5ed704bfe07874 mode="wmean" distance=3 power=2 cells=8 output=output676edb4ec1c744bfbae8680ec5156e44 uncertainty=uncertainty676edb4ec1c744bfbae8680ec5156e44 --overwrite

g.region raster=output676edb4ec1c744bfbae8680ec5156e44

r.out.gdal -t -m input="output676edb4ec1c744bfbae8680ec5156e44" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ad1109496a5e4dcdae0a48bd41c286d2\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

g.region raster=uncertainty676edb4ec1c744bfbae8680ec5156e44

r.out.gdal -t -m input="uncertainty676edb4ec1c744bfbae8680ec5156e44" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ca8beef361e14a739e8681a9f0d1d6c5\uncertainty.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

Starting GRASS GIS...

警告: マップセットの平行ロックはウィンドウズではサポートされていません。

一時ファイルを削除しています...

Executing <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> ...

C:\Users\AAA\Documents>chcp 932 1>NUL

C:\Users\AAA\Documents>g.proj -c proj4="+proj=tmerc +lat_0=36 +lon_0=138.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

警告: データム <Unknown_based_on_GRS80_ellipsoid> はGRASSに認識されていません そしてパラメータが見つかりません

Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default

Projection information updated

C:\Users\AAA\Documents>r.in.gdal input="C:\Users\AAA\Desktop\a\nuki.tif" band=1 output="rast_5ed704bfe07874" --overwrite -o

警告: データム <Japanese_Geodetic_Datum_2011> はGRASSに認識されていません そしてパラメータが見つかりません

Over-riding projection check

Importing raster map <rast_5ed704bfe07874>...

0..3..6..9..12..15..18..21..24..27..30..33..36..39..42..45..48..51..54..57..60..63..66..69..72..75..78..81..84..87..90..93..96..99..100

C:\Users\AAA\Documents>g.region n=-129837.2135 s=-138036.8719 e=47623.1745 w=38140.64 res=26.56172128851541

C:\Users\AAA\Documents>r.fill.stats input=rast_5ed704bfe07874 mode="wmean" distance=3 power=2 cells=8 output=output676edb4ec1c744bfbae8680ec5156e44 uncertainty=uncertainty676edb4ec1c744bfbae8680ec5156e44 --overwrite

W-E size of neighborhood is 7 cells.

S-N size of neighborhood is 7 cells.

Input data range is 0.000000 to 1079.000000.

Input data type is 'single' (4 bytes) and output data type is 'single' (4 bytes).

Minimal estimated memory usage is 0.026 MB.

Interpolating:

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

r.fill.stats 完了. Processing time was 0h0m1s.

C:\Users\AAA\Documents>g.region raster=output676edb4ec1c744bfbae8680ec5156e44

C:\Users\AAA\Documents>r.out.gdal -t -m input="output676edb4ec1c744bfbae8680ec5156e44" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ad1109496a5e4dcdae0a48bd41c286d2\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.

Checking GDAL data type and nodata value...

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

Using GDAL data type <Float32>

Exporting raster data to GTiff format...

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

r.out.gdal 完了. File <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ad1109496a5e4dcdae0a48bd41c286d2\output.tif> created.

C:\Users\AAA\Documents>g.region raster=uncertainty676edb4ec1c744bfbae8680ec5156e44

C:\Users\AAA\Documents>r.out.gdal -t -m input="uncertainty676edb4ec1c744bfbae8680ec5156e44" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ca8beef361e14a739e8681a9f0d1d6c5\uncertainty.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.

Checking GDAL data type and nodata value...

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

Using GDAL data type <Float32>

Exporting raster data to GTiff format...

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

r.out.gdal 完了. File <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\ca8beef361e14a739e8681a9f0d1d6c5\uncertainty.tif> created.

C:\Users\AAA\Documents>exit

Execution of <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> finished.

Cleaning up default sqlite database ...

Cleaning up temporary files...

処理は6.02秒で完了しました

結果:

{'output': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'System'}}>,

'uncertainty': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'System'}}>}

出力レイヤの読み込み

アルゴリズム 'r.fill.stats'が終了しました

以上です。

r.fill.statsはこの結果以下の2枚のレイヤが出力されました。

次にr.fill.nullsを実行した際のログを添付いたします。

QGIS version: 3.12.2-București

QGIS code revision: 8a1fb33634

Qt version: 5.11.2

GDAL version: 3.0.4

GEOS version: 3.8.1-CAPI-1.13.3

PROJ version: Rel. 6.3.1, February 10th, 2020

プロセシングアルゴリズム...

アルゴリズム 'r.fillnulls' を開始しています...

入力パラメータ:

{ 'GRASS_RASTER_FORMAT_META' : '', 'GRASS_RASTER_FORMAT_OPT' : '', 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : '38140.64,47623.1745,-138036.8719,-129837.2135 [EPSG:6676]', 'edge' : 3, 'input' : 'C:/Users/AAA/Desktop/a/nuki.tif', 'lambda' : 0.01, 'method' : 2, 'npmin' : 600, 'output' : 'TEMPORARY_OUTPUT', 'segmax' : 300, 'smooth' : 0.1, 'tension' : 40 }

g.proj -c proj4="+proj=tmerc +lat_0=36 +lon_0=138.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

r.in.gdal input="C:\Users\AAA\Desktop\a\nuki.tif" band=1 output="rast_5ed704a552b363" --overwrite -o

g.region n=-129837.2135 s=-138036.8719 e=47623.1745 w=38140.64 res=26.56172128851541

r.fillnulls input=rast_5ed704a552b363 method="rst" tension=40 smooth=0.1 edge=3 npmin=600 segmax=300 lambda=0.01 output=outputfa3d818f49be40a09f21c2b8383541a3 --overwrite

g.region raster=outputfa3d818f49be40a09f21c2b8383541a3

r.out.gdal -t -m input="outputfa3d818f49be40a09f21c2b8383541a3" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\5ff539e6a8984fa79eac57b154c19070\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

Starting GRASS GIS...

警告: マップセットの平行ロックはウィンドウズではサポートされていません。

一時ファイルを削除しています...

Executing <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> ...

C:\Users\AAA\Documents>chcp 932 1>NUL

C:\Users\AAA\Documents>g.proj -c proj4="+proj=tmerc +lat_0=36 +lon_0=138.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

警告: データム <Unknown_based_on_GRS80_ellipsoid> はGRASSに認識されていません そしてパラメータが見つかりません

Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default

Projection information updated

C:\Users\AAA\Documents>r.in.gdal input="C:\Users\AAA\Desktop\a\nuki.tif" band=1 output="rast_5ed704a552b363" --overwrite -o

警告: データム <Japanese_Geodetic_Datum_2011> はGRASSに認識されていません そしてパラメータが見つかりません

Over-riding projection check

Importing raster map <rast_5ed704a552b363>...

0..3..6..9..12..15..18..21..24..27..30..33..36..39..42..45..48..51..54..57..60..63..66..69..72..75..78..81..84..87..90..93..96..99..100

C:\Users\AAA\Documents>g.region n=-129837.2135 s=-138036.8719 e=47623.1745 w=38140.64 res=26.56172128851541

C:\Users\AAA\Documents>r.fillnulls input=rast_5ed704a552b363 method="rst" tension=40 smooth=0.1 edge=3 npmin=600 segmax=300 lambda=0.01 output=outputfa3d818f49be40a09f21c2b8383541a3 --overwrite

Using RST interpolation...

Locating and isolating NULL areas...

Growing NULL areas

Assigning IDs to NULL areas

エラー：Input map has no holes. Check region settings.

Execution of <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> finished.

Cleaning up default sqlite database ...

Cleaning up temporary files...

続行するには何かキーを押してください . . .

Starting GRASS GIS...

警告: マップセットの平行ロックはウィンドウズではサポートされていません。

一時ファイルを削除しています...

Executing <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> ...

C:\Users\AAA\Documents>chcp 932 1>NUL

C:\Users\AAA\Documents>g.region raster=outputfa3d818f49be40a09f21c2b8383541a3

エラー：ラスタマップ <outputfa3d818f49be40a09f21c2b8383541a3> が見つかりません

C:\Users\AAA\Documents>r.out.gdal -t -m input="outputfa3d818f49be40a09f21c2b8383541a3" output="C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\5ff539e6a8984fa79eac57b154c19070\output.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite

エラー：Raster map or group <outputfa3d818f49be40a09f21c2b8383541a3> not found

C:\Users\AAA\Documents>exit

Execution of <C:\Users\AAA\AppData\Local\Temp\processing_WlueKy\grassdata\grass_batch_job.cmd> finished.

Cleaning up default sqlite database ...

Cleaning up temporary files...

続行するには何かキーを押してください . . .

処理は17.32秒で完了しました

結果:

{'output': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'System'}}>}

出力レイヤの読み込み

次のレイヤは正しく生成されませんでした。<ul><li>C:/Users/AAA/AppData/Local/Temp/processing_WlueKy/5ff539e6a8984fa79eac57b154c19070/output.tif</li></ul>QGISメインウィンドウの"ログメッセージパネル"をチェックすると、アルゴリズムの実行に関する詳細情報が表示されます。

以上です。

r.fill.nullsのレイヤは結果として何も表示されませんでした。

度々迷惑をかけてしまい大変申し訳ありませんが、どうぞよろしくお願いいたします。