gdalを使ってラスタを作りたい

[i401...@gmail.com]

初めまして。

以下の投稿

https://groups.google.com/forum/#!topic/qgisshitumon01/bsbdidijRVI

を読んで

CSVファイルに、

x座標,y座標,B1,B2,B3　（BはBandのB)

のようなファイルを私も持っており、そのファイルからgeotiffをコマンドライン上から作りたいと思っております。

当該ページにはgdal_glidコマンドでgeotiffが作れること前提の会話になっており、細かいパラメータの設定方法などが残念ながら記載だれていませんでした。

(マルチバンドTiffを作るの前提の話だったみたいですね。)

そこで、gdal_glidを用いた場合の座標系や解像度の指定をどのように制御してやればいいのか、ご指導いただけたら幸いです。

QGISは2.18.20を使っているのですが、QGISで点群データとして読み込ますには2000万行程度ありハングアップしてしまいました。

そのため、GDALという機能でできればうれしいなという希望を持っています。

持っている点群の座標系は、青森県のデータで、1mグリッド間隔の解像度tiffを作りたいと思っております。1mグリッドの中に入っている点群は1点とは限らず

複数入っている場合にはその中の最小値、平均値、最大値をとった画像を作りたいです。

このような都合の良い方法があれば非常にありがたいのですが、良い方法がありますでしょうか？

[adachi]

コマンドではなくQGISを使っての操作ですが…

ファイルを読み込ませずともできるかと思います。

「ラスタ」メニュー　→　解析　→　グリッド（補間）

これが「gdal_grid」に相当します。

入力ファイルの右の「選択」をクリックすることで、読み込んでいるレイヤではなくパソコン内にあるファイルを指定できます。

いちおうアルゴリズムから「データメトリクス」を選択することで「最大」「最小」等選べるようになりますが

これらが実際にどういうアルゴリズムで値を設定するのか、恥ずかしながらよくわかっておりません。。

中途半端な回答ですみません、どなたか補足いただければ…

[ありた]

こんにちは

手元でテスト等できておらず、誤りが含まれているかもしれませんが。

まず gdal_grid のオプションについては下記のとおりです。オプション等はこちらに載っています。

https://www.gdal.org/gdal_grid.html

データとして X,Y,B1,B2,B3, ... という CSV データをお持ちとのことですが、

（簡単にキーワードを拾いつつ）読む感じだと B1, B2, B3, ... の平均等の処理の実行ではなく、

近傍（ある半径内の点群）の平均等といった補間アルゴリズムとしてのオプションのようです。

もし、  X,Y,B1,B2,B3, ... のように Z にあたるデータが複数ある場合や、 Z と判別できないような

名称（たとえば elevation ）のヘッダ行がある場合は、 -zfield オプションで明示します。

もしくは上記ドキュメントを参考に vrt ファイルを作る必要があるかもしれません。

B1, B2, B3, ... というデータがあるとして、これらのデータに対し平均等の処理を行いたい場合、

-sql オプションで演算させることで、不可能ではないかと思われますが、 Python や R などで

前処理を行ったほうがよいかもしれません。

座標系ですが、ドキュメントには出力ファイルの座標系を上書きする -a_srs のみの

記載しかないため、（ GIS としての）座標系は考慮されず、座標あるがままを認識し、解析し、

その素の座標値として出力ファイルが生成される気がします。

そのうえで、出力ファイルに（ GIS としての）座標系を割り当てるため -a_srs を使います。

つまり gdal_grid のなかでは（たとえば）経度緯度の入力データから平面直角座標系で 1m グリッドの

出力といった変換（再投影）を行うことはできなさそうです。

そのため、あらかじめ ogr2ogr などを利用して、入力データの座標変換を行う必要があると思います。

解像度ですが、 -outsize で出力のグリッド数（ピクセル数）を指定します。

またそれに対応する範囲指定ですが、ドキュメント末尾にある例では -txe と -tye が使われています。

そのほか、 -clipsrc や -spat は入力データを制限するオプションもありますが、境界における

補間の対象外になるかもしれません。（詳しい違いは把握してません）

不規則点群に対する補間ですが、 gdal_grid では近傍（半径1、半径2、傾きによる楕円内）でみつかる

点群に対し、 IDW （ invdist / 逆距離加重）や平均といった補間法や、最大値といった点群数といった

手法を選択することができます。ドキュメント末尾のコマンド例も参考にしてください。

https://www.gdal.org/gdal_grid.html#gdal_grid_algorithms

https://www.gdal.org/grid_tutorial.html

QGIS のラスタメニューあるいはプロセッシングで実行できるグリッドの解析は、下部の

コンソールをみると、直接 gdal_grid を実行しているようなので、ダイアログで指定できる

オプションの範囲内ですと、 QGIS からでも問題ありません。

ただし、 QGIS からではなく、直接 GDAL のコマンドを実行する際のメリットとして、

CSV データをレイヤとして読み込む負荷がない点、また上記細かなオプションを指定できるほか

メモリのキャッシュ上限や処理スレッド数の設定を変更することができます。

（ QGIS のダイアログ上で設定できる任意オプションは、出力ファイル生成オプション -co の

　key / value ペアのようで、その他のオプションを自由に追加することは難しそうです。）

gdal_grid では GDAL_NUM_THREADS が使えるとあります。

https://www.gdal.org/gdal_grid.html

また、その他の設定は

https://trac.osgeo.org/gdal/wiki/ConfigOptions

これらの設定は、環境変数として以下のように設定するか

set GDAL_CACHEMAX=4096

以下のように --config オプションを使用します。

gdal_translate --config GDAL_CACHEMAX 4096 in.tif out.tif

[i401...@gmail.com]

adachi様、ありた様

お答えありがとうございます。QGISからでもファイル指定ができたのですね。

今回のファイル、書き方もちょっと悪かった感じではあります。

実際には投稿文末尾のようなデータがあり、(実際にはこんなにきれいにｘｙは昇順降順に並んでいなく、順番すらもバラバラです。)

X,Y,RとX,T,GとX,Y,Bという3種類かつ、それぞれの画像のグリッドごとに含まれる点群の最大値、平均値、最小値の3種類、合計9種類が作りたいという意味合いでした。

イメージはhttps://qiita.com/mski_iksm/items/48c7d9d30233ad8c4f4f

のような感じで、右の緑色の点がグリッド内の最大、平均、最小になればいいなという感じです。(正確には、ポリゴンでいうと、ポリゴンの重心に左側の青い点の最大値、平均値、最小値を与えたいという質問内容です。)

そして、その時にこの値を計算してくれるグリッドサイズを自由に決めれるのかなと。

私自身、GISの知識はArcDesktopでベクトルを扱った知識ぐらいしかなく、点群やラスタ処理は素人同然です。GISが使えるならできるだろという業務命令で今回の処理にあたっております。そんな時にここの掲示板を見つけて投稿させていただきました。

OSはwindows10の64を使っております。構文を教えていただければ幸いです。

sampling.csv

X,Y,R,G,B
1,0,20.25060085,0.057214561,83.81735352
2,0,75.17559236,28.70552018,15.11220429
3,0,49.00856867,78.56766848,40.0182276
4,0,23.77558023,25.7974754,90.92075407
5,0,45.88342739,79.15951527,36.67598384
6,0,73.80046135,57.97756118,33.99453479
7,0,84.92431308,95.0682696,98.62427508
8,0,29.95630814,82.47577436,53.21702999
9,0,66.27211805,13.23983441,88.66759417
10,0,81.61791642,92.26151543,5.808710984
1,1,81.61904899,48.14488151,45.60155497
2,1,55.23254914,20.22754835,90.37982769
3,1,71.92841276,10.98504685,38.35280686
4,1,54.93763376,73.4277076,99.76393017
5,1,84.62196746,38.16199577,24.991753
6,1,63.14945192,6.096032362,35.29786325
7,1,39.07929293,75.17300995,60.45357998
8,1,38.43979177,57.21124517,6.593622189
9,1,62.70174371,74.1934077,76.46360947
10,1,59.14064693,50.6482346,5.404922172
1,2,67.11150272,12.97674121,60.53968467
2,2,82.07384268,92.30429027,90.03386836
3,2,51.04890846,54.93637755,80.11969542
4,2,19.47189368,72.57486547,65.9406836
5,2,78.41441889,97.59106583,25.27193303
6,2,95.597638,24.87065752,50.15175585
7,2,54.17365251,28.54976447,90.89241339
8,2,58.82471567,6.841408837,68.2988348
9,2,28.56922999,70.82493156,76.53219371
10,2,37.79007583,45.21032461,38.5938781

[ありた]

ご質問されている項目（対象とするデータ（ここでは R, G, B のいずれか）の指定、補間法（最大、

平均、最小）、グリッドサイズの指定）は昼の投稿に記載しておりますのでご参考ください。

また、具体的なコマンドの例はドキュメントの下部に載っております。

https://www.gdal.org/gdal_grid.html

QGIS のラスタ／プロセッシングで実行できるグリッドでも、ダイアログの下部に実行するコマンドが

表示されますので、参考になるかと思います。

そのうえで、先程の投稿を何点か補足いたします。

X,Y,Z というヘッダをつけた CSV ファイル（ VRT ファイル作成なし）でテストしてみたところ、

ジオメトリの認識失敗のメッセージがでたため、座標値を認識させるために CoordX,CoordY,Real と記載した

CSVT ファイルを作成する必要がありました。（ VRT ファイルを作成させる場合には不要）

（ ogr2ogr の CSV 入力で使用できる -oo X_POSSIBLE_NAMES 等は使用できないようです。）

同様にデータが Z という名前でも -zfield Z は必要なようです。

また、使用している GDAL のバージョンによってはエラーが発生することもありえます。

gdalinfo --version などで、バージョンの確認を行ってみてください。（現在の最新版は 2.4.0 ）

なお、2000万行と膨大なデータですので、あらかじめ Python 等で一定距離四方＋境界の補間のための余裕分

（例： -2m - 12m / 8m - 22m / 18m - 32m ... ）でファイル分割を行っておき、そのデータからそれぞれ

ラスタ化を行い、出力ラスタをマージするという手法のほうが処理時間等の点でも、必要となるスペックの

点でもよろしいかと思います。

[i401...@gmail.com]

ありた様

ご説明ありがとうございます。

csvファイルを短くして試しているのですが、どうにも画像ができてくれる感じがしておりません。

csvファイルがカンマ区切りであるということをgdalに認識させれていないのじゃないかと思っております。

カンマ区切りのcsvファイルであると認識させる方法はあるのでしょうか？

QGISの方で実行させる場合ですが、その時に「フィールド」からのZ値に列を入れることができませんでした。ファイルがカンマ区切りだと何かと都合が悪いということでしょうか？添付図は平均値がほしいという仮定のもと「平均距離」で計算は10m×10mの間にあるものにするため半径5mで探索させるとしました。

しかし、この場合の半径5mというのはどの座標を起点とした半径5mになっているのでしょうか？

VRTファイルというのも調べたのですが、CSVTファイルというもののほうが扱いやすかったのでそちらを作ることにしました。

コンピュータ自体は相当スペックに余裕があるものを使うため大丈夫かなとは思っていたのですが(XeonのCPUでメモリ64GB、Ubuntuシステム)分割したほうがいいのですか。

こちらの勉強からやることになるとは、相当な業務なんだと今日初めて実感しました。

[ありた]

こんにちは

「どうにも画像ができてくれる感じがしておりません」とありましても困ってしまいますので、

下端の GDAL のコマンドと、実行したあとのログタブに記載されている内容を載せていただくと

解決につながるかもしれません。

もしくは、レイヤは作成されたが、（たとえば）真っ白で意図する結果でないという意味なのでしょうか。

Z 値に選択できるものがない場合、数値としてではなくテキストの形式として認識された可能性もありえます。

無効値等で数値以外が入り込んだデータではないでしょうか。

読み込んだレイヤのプロパティを開き、フィールドタブで確認できます。

https://docs.qgis.org/2.18/ja/docs/user_manual/working_with_vector/vector_properties.html#fields-properties

指定する半径ですが、リサンプリングされる格子点からの半径となり、指定半径の楕円内に

含まれる点群が、補間アルゴリズムの対象となります。（半径 0.0 の場合は全点群）

https://www.gdal.org/grid_tutorial.html#grid_tutorial_ellipse

なお、添付画像では「平均距離」を選択しておりますが、点群までの距離の平均ですので、

おそらく目的の結果は得られないかと思います。 Z 値の平均は「移動平均」にあたります。

ラスタメニューからデータメトリクスのグリッドではなく、移動平均のものを選択してみてください。

データの分割の件ですが、もちろん十分なスペックの計算機を用意できるのであれば

問題ないかと思います。

[i401...@gmail.com]

ありた 様

CSVTファイルというもので、

CoordX,CoordY,Real

とした後に点群をQGISに読み込ませたら実験用のファイルからは画像が出力されるようになりました。

もう少し教えてください。

>指定する半径ですが、リサンプリングされる格子点からの半径～～

この部分ですが、リサンプリングは、格子の大きさだととらえているのですが、それでよいでしょうか？

また、リサンプリングされる格子の大きさを設定することってできるのでしょうか？今回は1m×1mの格子に含まれる点群の最大、平均、最小にしたいのですが。

補間アルゴリズムはわかりました。今回は格子に含まれるすべての点を対象にするため、ゼロ、ゼロがよさそうです。

平均値を使いたい場合は移動平均を使えばよいのですか。最大、最小の場合は「データメトリクス」の方で選べばよいでしょうか？

いろいろと質問だらけで申し訳ありません。

[ありた]

たとえば、 -txe 0 1000 -tye 0 1000 -outsize 100 100 で実行すると、 (0, 0) - (1000, 1000) の

領域を、サイズ 100x100 でリサンプリングせよという指示になるので、リサンプリング時の各点の座標は

(5, 5), (15, 5), (25, 5), ..., (5, 15), (15, 15), (25, 15), ... となります。この各点から指定した半径内に存在する

点群を補間時の計算対象として扱われます。

なお、この結果、作成された TIFF 画像は (0, 0) - (10, 10) でひとつのピクセルとして表示されるはずです。

ラスタは、メッシュ統計や土地利用など、一定の領域における代表値という扱いもしますが、

標高などでは、規則正しく並んだ点群データを扱いやすくするため画像として保存しているパターンも

あります。この場合、画像のピクセルは隣の点まで隙間を埋めて表示しているくらいのイメージです。

私も「格子点」と、各ピクセルの四隅の点のような書き方をして混乱させて申し訳ありません。

ですので、格子が解析の対象範囲とかそういった意味合いはここではないはずです。

なお、ドキュメントを読む限り、対象範囲の指定方法は楕円のみで、矩形で計算範囲を指定することは

（少なくとも） gdal_grid では対応していないようです。

また、その際、半径を 0, 0 として指定すると全点群を対象になりますので、この場合、2000万データすべてを

対象とするという意味になるかと思われます。（すみませんテストはできていません。）

それでも逆距離加重であれば、遠いほどマイナスの重みが付与されますので、近傍データの影響を強く受けますが、

平均の場合ですと全点において同じ値となります。

最大、最小はデータメトリクスの中にあるものですね。

# 今日、 QGIS 2.18 の「グリッド（補間）」ダイアログを確認しましたが、解像度（ -txe -tye -outsize ）が

# 指定できますね。というより、 QGIS 3.x ではダイアログからできないのはなんでなんだろう……

# ついでにコマンド自体の編集もできるので、 --config の指定もできるのでは

[i401...@gmail.com]

ありた様

いろいろとありがとうございます。

だんだんと形になってきております。

もう一つ教えてください。

今回ですが、

-txe 0 1000 -tye 0 1000 -outsize 100 100

と書いていただいておりますが、点群から点群が入っている領域の4隅を取得するのに良い方法はありませんでしょうか？

点数が多すぎて表計算上でソートするのも一苦労になっておりまして。その3DCGなどではバウンダリボックスというのですが、点群領域が入るであろう

4隅を知りたいなということです。

このようなときに良い方法があれば教えてください。

[ありた]

こんにちは

（仕事の関係で）遅くなりました。

少なくとも gdal_grid を使う場合、矩形で対象を指定できないため、バウンダリボックスは（計算しても、

補間対象は楕円内であることに変わりはなく）直接関係ないため、意図するところがちょっとわかりませんが、

TIFF 画像でいうところのピクセルの四隅の座標ということでしょうか。

-txe xmin xmax -tye ymin ymax -outsize width height としたとき、各ピクセルの

横幅は w = (xmax - xmin) / width 、縦幅は h = (ymax - ymin) / height なので、

横 i 番目、縦 j 番目のリサンプリング点（ピクセル）の四隅の座標は

w * (i-0.5), h * (j-0.5), w * (i+0.5), h * (j+0.5) といったことでしょうか？

[i401...@gmail.com]

ありた様

前回お聞きしたことは、そうです。作成されるtiff画像の4隅の座標を知った方が良いのかなと思っての質問でした。あまり意味がないもの？なのですね。

色々とご指導ありがとうございました。

点群データの図化処理、うまく実行することができました。

お返事が遅くなり大変申し訳ありませんでした。

QGISおよび付属しているコマンドライン群、すごいことができるものなのですね。

そう考えるとすごく高いE社の商用GIS、いらないんじゃないかと思えてきました。なんだかんだで使ってきたソフトにこういうのなんですけど

バイナリデータばっかりで書式も公開していないもの、今後のバージョンアップでいきなり互換が切られかねなくて怖くて使っていられない。

ほとんどの機能がQGISで代替あるいはもっと上位の事が実施できるしこうやってコミュニティだってある。私はこちらに乗り換えようと思います。

しいて言うなら、ラスタに日本語使うとエラーが起きてしまう事ですけど、コピーしてリネームすれば対応できるので些細な問題ですね。

今回は教えてもらうだけという事になって申し訳ないですが、年度明けからいろいろと勉強して、回答できる知識を得ようと思いました。

これからもいろいろとご指導お願いします。