QGIS:属性テーブルの文字列を置換する

 GISの属性テーブルはデータベースと理解していたので、中身の文字列を置換できるとは思っていませんでした。ので、チマチマ個別に入力したり、CSVで出力→Excelで編集→CSVを読み込み、みたいな面倒くさいことをしていました。今回、数千件の文字列の書式を変換する作業が発生したので、いろいろ調べてみたら、普通にできることが分かったので共有します。

参考1:https://qiita.com/kntoshiya/items/217cab81e4e6b37d8918

参考2:https://groups.google.com/g/qgisshitumon01/c/oWX9suovJM0?pli=1

 サンプルデータは国土数値情報の公共施設@富山県。このポイントデータの属性テーブルのフィールドは、行政区域コードや分類、名称、所在地などからなっています。

公共施設データ。行政区域コードで色分けした。
属性テーブルのフィールド。
属性テーブルの一部。P02 _009は管理者。

 ↑に示したように、P02_005やP02_006には分類や名称が記述されています。ここではまず、P02_005の全角カッコを半角にしてみます。

 まず、このファイルを編集可能にします。属性テーブルの鉛筆マークをクリックするか、レイヤを右クリックして編集モードを切り替えます。次に、フィールド計算機を開きます。

左上の鉛筆ボタンをクリックして編集可能にする。
そろばんぽいボタンをクリックしてフィールド計算機を開く。

 フィールド計算機を開いたら、「式」のスペースに以下のように入力します。

replace( "P02_005" ,'(','(')

 replaceの書式はreplace(string, before, after)となっています。stringは文字列で、今回の場合はP02_005を指定しています。次に、beforeとして全角左カッコを、afterとして半角左カッコを指定しています。基本はこれだけで、あとは新たにフィールドを作るか、既存のフィールドを更新(上書き)するかを選びます。OKをクリックすると、全角左カッコが半角に変換されます。

 なお、私の環境では式のフィールドで全角入力ができないので、テキストエディタなどに入力したものをコピペしています。全角での入力ができない場合は、コピペを試してみてください。

フィールド計算機のダイアログ。「既存のフィールドを更新」にチェックを入れている。
置換結果。左カッコが半角になっている。

 右カッコを全角にするには、もう一度replaceをかければ良いのですが、複数の置換を一度に行うことも可能です。例えば左右カッコを一度に半角にするには、次の式を入力します。

replace( replace("P02_005" ,'(','('),')',')')

 つまり、replaceを重ねがけすれば良いということです。結果は↓の通りで、両カッコが半角になっています。また、隣りのP02_006には影響を与えていません。もちろん、replaceはいくつでも重ねられます。

置換結果。

 もちろん、ある文字を削除することもできます。例えばカッコを削除したいと思ったら、

replace( replace("P02_005" ,'(',''),')','')

 とすれば、↓のようにカッコが削除されます。

 アンダーバーで繋いでみましょう。以下の式を入力します。

replace( replace("P02_005" ,'(','_'),')','')

 結果は↓の通りです。

 これで、文字列の細かい置換ができるようになりました。その他、役に立ちそうな式も紹介しておきます。

"P02_005" || '_' || "P02_009"

 P02_005とP02_009をアンダーバーで繋ぎました。

 定型文をくっつけることもできます。

concat( "P02_005" ,'_富山県')

 とすれば、「富山県」という文字列が各セルの末尾につきます。ちなみに、これと一つ前の式でやっていることは同じなので、例えば

concat( "P02_005" , '_', "P02_009")

 あるいは

"P02_005" || '_富山県'

 と書くこともできます。定型文をくっつけるなんて、普通はやらないですよね。が、例えば、全てのセルの末尾に単位をつけるといった用途に使えます。フィールド計算機は奥が深いので、また何か見つけたら共有します。

 最後に、文字列、とくに記号は色々なところに潜んでいるので、一括置換すると思わぬ結果を生む場合があります。置換対象のフィールドをよく見て式を適用してください。

 

QGISのマップを3Dで表示する

 QGISには3Dマップビューという機能があり,そちらも使いこなせば便利なんだと思うんですが,私はちょっと使いにくいので,別の方法を紹介します.

 今回紹介する方法は,2つのプラグインを使います.まずはメニュー/プラグイン/プラグインの管理とインストールをクリックし,Qgis2threejsElevationTile4JPというプラグインをインストールします.前者は3Dビューの,後者はDEM(数値標高モデル)をダウンロードするためのプラグインです.作者の皆さまに感謝いたします.

プラグインのインストール画面.検索窓に単語を入力すると見つけやすい.
プラグインのインストール画面.目的のプラグインを選択し,右下のインストールをクリックする.

 プラグインをインストールしたら,まずDEMをダウンロードします.何らかの地図タイルを表示して,3D表示したい範囲にズームしましょう.

地理院タイル「淡色地図」で富士山周辺を表示してみた.

 この範囲のDEMをダウンロードします.メニュー/プラグイン/ElevationTile4JPをクリックします.操作はかなりシンプルで,ズームレベルと座標系を設定し,出力先を決めればOKです.ズームレベルの解説は国土地理院のページをご覧ください.基本的には,レベルを上げると解像度が上がりますが,範囲が広すぎるとエラーが出ます.その場合はズームレベルを下げましょう.ちなみに,上の範囲(縮尺5万分の1)だとレベル12でした.座標系は(地形解析などに使うのでなければ)何でも良いと思います.最後にOKをクリックすると,ダウンロードが始まります.ダウンロードが終わると,DEMが表示されます.が,DEMは何もしないとグレースケールで表示されるので,3D映えしません.レイヤパネルで先ほどの地図レイヤをDEMの上に持ってきましょう.

ElevationTile4JPのダイアログ.出力ファイルの…をクリックして出力先とファイル名を入力する.
ダウンロードされたDEM.何もしないとグレースケールで表示される.

 いよいよ3D表示に移ります.メニュー/Web/Qgis2threejsを起動します.起動したては,なんか空みたいな画面で肩透かしをくらうと思います.左側のパネルの上の方に,DEMレイヤが表示されているので,先ほど保存したDEMファイルをチェックします.これで,マップが3D表示されるはずです.

3D表示.

 画面の操作は,左クリック&ドラッグで回転,右クリックで移動,スクロールでズームになります.鉛直倍率が足りない場合は,QGIS2threejsのメニュー/Scene/Scene Settingsで,Z exaggerationを大きくしましょう.デフォルトは1倍になっているので,例えば2,3といった具合に数値を増やします.上の図の例では,2倍にしています.他の設定もScene Settingsでできます.

 表示画面の保存は,QGIS2threejsのメニュー/File/Save Scene As/Image (png)でできます.同じメニューでgltfを選ぶと,3Dモデルをまるごと保存できます.このモデルはこちらのページにドラッグすると閲覧できます.

 ちなみに,3D表示はマップに表示されているものを反映します.マップ側を変更すると,3D表示も変わります.もちろん,ラスタやタイルだけでなく,ベクタも表示できます.なお今回は簡単のため,プラグインを使ってDEMをダウンロードしましたが,手持ちのDEMでも同じことができます.

シームレス地質図+地理院タイルの例.
活断層図+傾斜量図の例.断層変位地形がよく分かる.

QGISで段丘面を抽出する(試行)

 段丘面区分は,地形判読の大事な要素なんですが,実際やるとなると手作業の手間がかかって大変です.こんなこと言うと怒られそうですが,自動でベクタ化(ポリゴン化)と区分ができると良いいなと思っていました.すごいマニアックなテーマですが,ある程度うまくいったのでその方法を公開します.

 対象地域は日本を代表する河岸段丘が見られる群馬県の沼田付近です.国土地理院による日本の典型地形でも紹介されている場所です(https://www.gsi.go.jp/kikaku/tenkei_kasen.html#河岸段丘及び段丘崖).QGISで地理院タイルを表示して,沼田付近に移動してみます.

沼田付近.地理院タイル(傾斜量図)を使用.利根川水系沿いに多段化した段丘面が広がる.右下は赤城山.

 この場所のDEMを用意します.今回はMIERUNEさんが公開されている,ElevationTile4JPというプラグインを使って国土地理院のDEMをダウンロードしました.このプラグインの概要やインストール方法はこちらで紹介されています.すっごい便利です.感謝です.

ElevationTile4JPのダイアログ.
ダウンロードしたDEM(ラスタ).たぶん10mメッシュDEM.

 QGISのメニュー/ラスタ/解析/傾斜を使って傾斜量図(ラスタ)を作ります.今回はDEMを平面直角座標系で出力したのですが,地理座標系の場合はGRASS GISのr.slope.aspectを使う必要があります.傾斜量図を作ったら,レイヤプロパティ/シンボルで色分けしてみましょう.段丘面の傾斜は概ね数度なので,上限を10°とかにして色分けすると,段丘面の範囲が強調されます.ポイントは,地形面の境界(崖)がはっきり押さえられる傾斜を見極めることです.この範囲の場合は,6°で段丘面境界が押さえられると判断しました.

 ちなみに,下の図を見ると,下位の面ほど明るく,平均的な傾斜が緩いことが分かります.後ほど,これを使って段丘面区分にもチャレンジします.

傾斜量図.色が濃いほど斜度が高い.

 この傾斜をしきい値として,2値化したラスタを作ります.メニュー/ラスタ/ラスタ計算機で,傾斜6°以下のピクセルを1,それ以外を0とするラスタを生成します.計算式は以下の通りです.

"傾斜量ラスタ" <= 6
ラスタ計算機のダイアログ.
結果.傾斜6°以下の場所が1となり,白く表示される.

 この2値化したラスタをポリゴンにしていきます.まず,ラスタの各ピクセルをポリゴン化します.私はメニュー/プロセシング/ツールボックスから呼び出せる「ラスタをベクタ化(pixels to polygons)」を使いました.実行すると,全ピクセルが四角形のポリゴンになり,各ピクセルにはラスタの属性値が与えられます.今回の場合,0か1の属性値が与えられることになります.

pixels to polygonsのダイアログ.対象は2値化したラスタ.
実行後.ラスタの全ピクセルがポリゴン化される.
上のポリゴンを0と1で色分けしたもの.

 このポリゴンに,メニュー/ベクタ/空間演算ツール/融合(dissolve)をかけます.これにより,同じ属性値を持つセルが結合され,段丘面ごとのポリゴンが生成されます.結構きれいに抽出できて嬉しいです.

「融合」ダイアログ.
実行結果.面ごとにポリゴンが分かれる.

 これで段丘面の抽出ができました.次は,区分にチャレンジしてみます.先ほど指摘しましたが,下位の面ほど傾斜が緩いことを利用すれば区分できるのでは,と考えました.そのために,面のポリゴンごとに傾斜の統計量をとります.ツールボックスからGRASS GISのv.rast.statsを呼び出し,実行します(SAGA GISのRaster statistics for polygonsでも代用できます).これにより,傾斜量の統計量(平均や分散など)が付与された新しいポリゴンファイルが生成されます.

v.rast.statsダイアログ.
v.rast.statsの実行結果.

 ここまで来たらあと一歩です.まずクエリで段丘面以外の部分を非表示にします.ポリゴンファイルのレイヤプロパティを開いて,設定タブからクエリビルダを呼び出します.クエリの式は以下の通り.

"VALUE" = 1
レイヤプロパティ/設定タブ.
クエリビルダのダイアログ.

 最後に,レイヤプロパティのシンボルで傾斜の統計量で色分けします.下の図は,中央値を連続値で色分けしたものです.下位の面ほど緩い傾向が出せたと思います.

今のところの結果.

 以上,手間はかかりますが,段丘面の抽出と区分がなんとかできそう,な感じです.ただ,いくつか課題はあります.

・上の図で顕著だが,意図しない境界ができている(右側).融合がうまくできていない?DEMのデータ配列のせい?

・本当はもっと細かく分かれているはず.より解像度の高いDEMを使えば抽出できる?あるいは,傾斜よりうまく分かれる地形量がある?

 まだ完成にはほど遠いですが,今後改善に向けて取り組んでいきます.2値化やラスタのポリゴン化,ポリゴン統計など,他の目的でも活用できそうな手法もあるのでご紹介しました.