QGISのフィールド計算機などにおける式についてまとめてみた

:::note info
これは [MIERUNE AdventCalendar 2022](https://qiita.com/advent-calendar/2022/mierune)　20日目の記事です。
昨日は @northprint さんによる [Cesium for Unity クイックスタート](https://qiita.com/northprint/items/ba9c0760ad359cadf66d) でした。
:::

## **はじめに**

QGISでは、SQLのような式を使ってデータ処理できる機能が用意されています。
式を使うことで、レイヤのスタイルやラベルといったシンボロジに関わる部分を制御したり、属性テーブルの値やジオメトリを処理したりするといったさまざまな処理が可能です。

意外に便利なのがこの式ですが、日本語の解説記事はあまり多くないように思うため、この記事では式の使い方についてまとめてみたいと思います。

## 式を使用できる場所
以下では代表的な場所をいくつか挙げてみます。なお、この他にもラベルだったりシンボロジなどさまざまな場所で使用可能です。

### **フィールド計算機**

まず一番利用することが多い、フィールド計算機です。

フィールド計算機では、フィールドのデータ型を文字列から数値に変換したり、複数のフィールドの文字列を結合したりといったスプレッドシートの関数と同じようなことが可能です。それ以外にも、ラインの長さやポリゴンの面積を算出したり、ジオメトリの重心にポイントを生成といったようなGIS特有の関数も用意されています。

**●使用方法**
レイヤの属性テーブル＞フィールド計算機を開く
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/cdac9667-fd32-5a45-d45e-cdadc01ebdff.png)

### 式によるジオメトリ
式によるジオメトリは、式を使って新しいジオメトリを作成することが可能です。

たとえば、ラインの始点や終点にポイントを作成するなど、プロセシングツールでやろうとすると以外に面倒な処理が簡単に実現可能だったりします。

**●使用方法**
プロセシングツールボックス＞ベクタジオメトリ＞式によるジオメトリ
![Untitled 1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/bdffa632-e562-2770-10c0-147432a781ff.png)


### **ジオメトリジェネレータ**

ジオメトリジェネレータは、式によるジオメトリと同様に式を使ってジオメトリを作成できます。

式によるジオメトリと異なるのは、ジオメトリジェネレータで作成されるジオメトリは、既存レイヤのシンボルとしてジオメトリが生成される点です。

なので、ジオメトリジェネレータでラインの始点にシンボルを表示させていた場合、ラインの始点の位置を修正すればシンボルの位置も動的に変更されることになります。

**●使用方法**
レイヤのシンボロジ＞シンボルレイヤタイプ＞ジオメトリジェネレータ
![Untitled 2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/97fda825-5f6b-dd7a-0a4f-9a20f49c36ba.png)

### **QGIS式による選択 or 抽出**
フィールドに特定の値が入っている地物のみ選択や抽出したい場合に利用します。例えば、面積フィールドが5,000㎡よりも大きい（”area” > 5000）地物のみを抽出したいといった場合に有用です。

**●使用方法**
プロセシングツール＞ベクタ選択＞QGIS式による選択or抽出
![Untitled 3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/d4070d51-18df-ada8-039f-caa93eab223f.png)

## 式文字列ビルダー
上記の場所に配置されている`ε`のボタンをクリックすると「式文字列ビルダー」という式を記述するためのウィンドウが表示されます。
中央の枠内にはQGISに組み込まれている関数が分類されており、ここから目的の関数を探すことができます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/aa80603a-baf6-3d34-576e-5de880c4f07e.png)

## **フィールド計算機**
それでは、ここからはフィールド計算機の使い方を紹介します。

### **フィールドの設定**

フィールド計算機を使用する場合、まずはレイヤのどこのフィールドに処理結果を書き出すのかを設定します。

新しいフィールドに結果を出力したい場合は、「新規フィールドを作成」に✓をいれます。
この場合、フィールド型（文字列、整数、少数など）を求める結果に応じて適切に入力する必要があります。出力結果は文字列なのに、フィールド型が数値となっているなど型が一致していない場合は、結果がNULLとなるため注意が必要です。

既存のフィールドに結果を出力したい場合は、「既存のフィールドを更新する」に✓をいれて、出力先のフィールドを選択してください。
![Untitled 4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/0ef918ab-01ca-b369-0b0c-7daa69efaaec.png)

なお、地物を選択した状態でフィールド計算機を開くと、左上の「選択されているn個の地物のみ更新する」がデフォルトで✓が入る仕様になっています。

これに✓をいれた状態でフィールド計算機を実行すると選択地物のみしか処理結果が反映されません。そのため、想定外の結果となるのを避けるためにも地物の選択を全解除したうえで実行するのが安全でしょう。

### **基本的事項**
フィールド計算機を使用するにあたっての基本的事項です。

```sql
--何も囲わないと数値
123

--単一引用符で囲むと文字列
'string' 

--二重引用符で囲むとフィールド名
"field"
```

フィールドの値を使用した計算例です。
```sql
--平方メートルの数値を平方キロメートルに変換
"area" /　1000000

--人口フィールドと面積フィールドから人口密度を計算
"population" / "area" 
```

式にはコメントの記述も可能です。次の2通りの方法があります。

```sql
-- コメント記載
/* コメントを記載 */
```

### **データ型の変換**
特定のカラムの値を文字列から数値に変換したり、数値を文字列に変換するなどが可能です。

```sql
--文字列のフィールドを数値に変換
to_int(string) --整数型
to_real(string) --実数型

--数値のフィールドを文字列に変換
to_string(number) 

--数値のフィールドを文字列に変換し3桁で0埋め
lpad(to_string(number),3,'0')
```

### **文字列操作**
文字列を操作する関数です。

```sql
--文字列を大文字に変換
upper(string) 
--文字列を小文字に変換
lower(string)

--文字列の左端のn文字分を返す
left(string, n)  
--文字列の右端のn文字分を返す
right(string, n) 

--文字列を指定した別の文字列で置換した文字列を返す
replace(string, before, after): 

--文字列の文字数を返す
length(string)

--複数の文字列を1つに結合したものを返す
concat(string1, string2)
```

### **条件分岐による値の入力**
```sql
--条件をテストし、それに応じて異なる結果を返す
if(チェックする条件,正の場合の結果,負の場合の結果)

--cdをテキストに置換する（国土数値情報の都市地域データでの例）
case
	when "layer_no" = 1 then '市街化区域'
	when "layer_no" = 2 then '市街化調整区域'
	when "layer_no" = 3 then 'その他用途地域'
	when "layer_no" = 4 then '用途未設定'
end

--価格に応じてランク付け（国土数値情報の地価公示データでの例）
case
	when "L01_006" < 100000 then 'low'
	when "L01_006" < 300000 then 'mid'
	else  'high'
end

```

- 都市地域データ
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A09.html](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A09.html)
- 地価公示データ
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_1.html#prefecture01](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_1.html#prefecture01)

### **ジオメトリ関係**
それぞれの地物に応じた結果が返されます。

```sql
--現在の地物のジオメトリを返す（単体では使わず、他の関数と組み合わせて使用する）
$geometry

--ポリゴンの面積を返す
$area

--ラインの長さを返す
$length

--ポリゴンの周長を返す
$perimeter

--ポイントのXY座標を返す
$x
$y

--ジオメトリのwkt表現を返す
geom_to_wkt($geometry)
```

次のように複数の関数を組み合わせることも可能です。

```sql
--WGS84から平面直角座標系に変換し面積を算出
area(transform($geometry,'EPSG:4326','EPSG:6675'))

--ラインやポリゴンの重心点のXY座標の取得
x(centroid($geometry))
y(centroid($geometry))

--平面直角座標系のポイントレイヤをWGS84に変換して、緯度経度をカンマで区切りで取得（||は文字列の結合）
y(transform($geometry,'EPSG:6675','EPSG:4326'))|| ','||x(transform($geometry,'EPSG:6675','EPSG:4326'))

--二点のポイントから角度を計算（ここではラインの始点と終点を使用）
degrees(azimuth(start_point($geometry),end_point($geometry)))
```

### **空間結合**
式を使って空間結合を行うことも可能です。
なお、大量の地物を処理しようとすると結構時間がかかるので、その場合はプロセシングツールの空間結合を使いましょう。

ここで使用する`overlay関数`は一対多の結合となる場合が想定されるためデータが配列で返ってくる仕様となっています。

配列は直接フィールドに出力できないため、一対一の結合の場合は[0]を式の後ろにつけて配列の先頭の要素を指定をします。一対多の場合は`array_to_string`で配列を文字列に変換させてから書き出しましょう。
```sql
--空港から最も近い駅名を取得
overlay_nearest('駅',"N02_005")[0]

--空港を含む行政区域から都道府県を取得
overlay_within('行政区域',"N03_001")[0]

--都道府県の行政区域に含まれる空港を抽出(一対多の結合)
array_to_string(overlay_contains('空港',"C28_005"))
```

- 行政区域データ（都道府県名でディゾルブしてから使用）
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html#prefecture00](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html#prefecture00)
- 駅データ（鉄道データ内）
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v3_0.html](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v3_0.html)
- 空港データ（ポイントに変換して使用）
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C28-v3_0.html](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C28-v3_0.html)

### **集計関数**
プロセシングツールの「空間結合（集計付き）」と概ね同様の機能です。

```sql
--市区町村の行政区域に含まれる地価公示のポイント数を集計
aggregate(
	layer:='地価公示',
	aggregate:='count',
	expression:="L01_001",
	filter:=contains(geometry(@parent), $geometry)
)

--市区町村の行政区域に含まれる地価公示のうち一番高い価格を取得
aggregate(
	layer:='地価公示',
	aggregate:='max',
	expression:="L01_006",
	filter:=contains(geometry(@parent), $geometry)
)
```

- 行政区域データ
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html#prefecture00](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html#prefecture00)
- 地価公示データ
    - [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_1.html#prefecture01](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_1.html#prefecture01)

### その他

最後は、施設を示すポイントデータと道路を示すラインデータから、その施設が接面する道路の方位を取得してみます。
`with_variable`は第3引数の式で使われる変数を設定することができる関数です。第1引数で設定した変数名の前に`@`をつけて使用します。
```sql
--施設（①）と施設から最も近い道路への最短となるラインの終点（②）を基に角度を計算
with_variable(
	'angle',
	degrees( --ラジアンから度に変換
		azimuth( -- ①から見た②の方位をラジアンで返す
			$geometry, --①
			end_point( --②
				shortest_line(
					$geometry, 
					overlay_nearest('line',$geometry)[0])
				)
			)
		),
--case文で角度に応じて方位を入力する
case
	when @angle < 22.5  then '北'
 	when @angle < 67.5  then '北東'
	when @angle < 110.5 then '東'
	when @angle < 157.5 then '南東'
	when @angle < 202.5 then '南'
	when @angle < 247.5 then '南西'
	when @angle < 292.5 then '西'
	when @angle < 337.5 then '北西'
	else  '北'
end
)
```
![Untitled 5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/3eebde9b-82c6-1415-0858-1e36f5aa9b2b.png)
※施設が角地にある場合などは、適切に取得できない可能性があります。

## ラベル
ラベルでの式の利用例です。値の欄に直接式を入力することで、ラベル表示用に新規フィールドを作成することなく目的のラベルを表示することができます。
```sql
--人口フィールドを文字列化して単位をつけて表示
to_string("JINKO") || '人'
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/862f927f-2f56-a63d-8af9-7f5022271a95.png)

- 境界データ | 13112 世田谷区（e-stat）
https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2020&serveyId=A002005212020&prefCode=13&coordsys=1&format=shape&datum=2000


## ジオメトリ作成関係
以下の式の結果はすべてジオメトリを返すので、「式によるジオメトリ」と「ジオメトリジェネレーター」どちらでも動作します。下記はジオメトリジェネレーターで表示している例です。

```sql
--100mのバッファを作成
buffer($geometry,100)
```
![Untitled 6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/d9efb8f2-ab1f-5985-8e2d-2566cbf687b2.png)


```sql
--ラインの始点にシンボル作成（青色）
start_point($geometry)

--ラインの終点にシンボル作成（赤色）
end_point($geometry)
```
![Untitled 7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/f34ba342-13d7-02d1-dc35-cdc82514184f.png)


```sql
--ラインの重心点にシンボル作成（青色）
centroid($geometry)

--ラインの中点にシンボル作成（赤色）
line_interpolate_point($geometry,$length/2)
```

![Untitled 8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/38a859b2-f7b0-7e6d-c805-bdacc7f47d14.png)

```sql
--ラインの指定した頂点にポイント作成
make_point($x_at(1),$y_at(1))
```
![Untitled 9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/8dc6d5e9-9044-3238-9cfd-1ac379295662.png)


```sql
--ターゲットレイヤの最も近い地物に対して最短となるラインを返す（黒の破線）
shortest_line($geometry, overlay_nearest('line', $geometry)[0])
```
![Untitled 10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/4513fd1a-2afd-8d8d-74b4-5a86e11acc6e.png)


```sql
--元のラインから50メートル間隔でラインを作成（黒の破線）
collect_geometries (
    array_foreach (
        generate_series(50,200,50), --50～200まで50間隔で配列を作る
        offset_curve($geometry, @element)
    )
)
```
![Untitled 11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/0ceb4245-5044-2a12-4161-0b24aeb1acbb.png)

```sql
--ポイントから16方位を描く（generate_seriesの引数を”0,315,45”にすれば8方位）
collect_geometries(
	array_foreach(
		generate_series(0,337.5,22.5),
		make_line(
			$geometry,
			project($geometry, 100, radians(@element))
    		)
    	)
    )
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/f7bc5a45-1ae3-1c5d-e966-421b53cf66fb.png)


最後に、[こちらの記事](https://static.geotribu.fr/articles/2022/2022-07-11_qgis_joy_division/)で紹介されていた、Joy Divisionの「Unknown Pleasures」のジャケットのような地図を作る、というのを試してみます。

このような図は「Joy Plot」と呼ばれているようで、30DayMapChallengeでも[作成している方](https://twitter.com/anhsirkwijaya/status/1587877272802902016?t=C8c6Er8uSBcmYx1u6EmxFQ&s=19)がいたり、[PythonやRでも作成できたり](https://qiita.com/kazetof/items/a9645951e9579991cadb)するようです。

羊蹄山付近で試してみると以下のような結果になりました。

![Untitled 12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/670401/ae7ef9d6-d949-5166-8aa5-68be001a5b67.png)


## おわりに
以上、式を使ってできることをいろいろと紹介させていただきました。式では思った以上に多くのことに使えるので覚えておくといずれ役に立つ場面があるかと思います。

この記事を書いている時に知りましたが、QGISの公式ドキュメントに式や関数について詳しくまとまっていました。こちらも合わせて見ていただくと、より理解が深まると思います。

- QGISのドキュメント | 式について
https://docs.qgis.org/3.22/ja/docs/user_manual/expressions/expression.html

- QGISのドキュメント | 関数一覧
https://docs.qgis.org/3.22/ja/docs/user_manual/expressions/functions_list.html#functions-list 

また、こちらの記事もおすすめです。
- QGISのフィールド演算機でできること
https://medium.com/@mitsuhamiyake/qgisのフィールド演算機でできること-530e7c86ad12

- QGISのフィールド演算機を使いこなす
https://medium.com/@mitsuhamiyake/qgisのフィールド演算機を使いこなす-66459176c822

- QGISのフィールド演算機でデータ型変換と文字列操作
https://medium.com/@mitsuhamiyake/qgisのフィールド演算機を使いこなす-2-31eb889a4502

- QGISのフィールド演算機でジオメトリの属性を調べる、ジオメトリジェネレータでシンボルを表現する
https://medium.com/@mitsuhamiyake/qgisのフィールド演算機でジオメトリの属性を調べる-ジオメトリジェネレータでシンボルを表現する-c076dbd31450

:::note info
明日は@asahina820さんです！お楽しみにー
:::