Automating GIS Processes 2024 写経　Lesson 3（Spatial join、前半）

https://autogis-site.readthedocs.io/en/latest/lessons/lesson-3/spatial-join.html

> Spatial joins are operations that combine data from two or more spatial data sets based on their geometric relationship. In the previous sections, we got to know two specific cases of spatial joins: Point-in-polygon queries and intersects-queries. However, there is more to using the geometric relationship between features and between entire layers.

Spatial join、空間結合＝複数の空間データセットをそれぞれのgeometryの「関連」で結合、です。
ここまで、Point-in-polygonやintersetを実践してきましたが、他にもあります。

> Spatial join operations require two input parameters: the predicament, i.e., the geometric condition that needs to be met between two geometries, and the join-type: whether only rows with matching geometries are kept, or all of one input table’s rows, or all records.

空間結合は、２つパラメータがあります。
- 述語（the predicamentと書いてあるが、the predicateと思われる）：２つのgeometryがどういう関係か
- 結合方法：一致したgeometryがある行だけを残すか、片一方のテーブルの行はすべて残すか、両方のテーブルの行をすべて残すか

> Geopandas (using shapely to implement geometric relationships) supports a standard set of geometric predicates, that is similar to most GIS analysis tools and applications:

`Geopandas`(`shapely`でgeometryの「関連」を実装)は、ほとんどのGIS分析ツールおよびアプリケーションとおなじような[標準的なgeometryの「条件」のセット](https://geopandas.org/en/stable/docs/user_guide/mergingdata.html#binary-predicate-joins) をサポートしています。

- intersects
- contains
- within
- touches
- crosses
- overlaps

> Geometric predicaments are expressed as verbs, so they have an intuitive meaning. See the shapely user manual for a detailed description of each geometric predicate.

（ここのpredicamentsもpredicatesの誤字と思われる。）
述語は動詞で直感的なものです。[shaplyの公式マニュアル](https://shapely.readthedocs.io/en/stable/manual.html#binary-predicates)に詳細があります。

:::note info
Binary geometric predicates

> Shapely supports more binary geometric predicates than those geopandas implements for spatial joins. What are they? Can they be expressed by combining the implemented ones?

`shapely`は`Geopandas`よりも多くの空間結合の「述語」サポートしてます。
では、`Geopandas`の「述語」だけで`shapely`と同じことができるでしょうか？
:::

> In terms of the join-type, geopandas implements three different options:
> - left: keep all records of the left data frame, fill with empty values if no match, keep left geometry column
> - right: keep all records of the left data frame, fill with empty values if no match, keep right geometry column
> - inner: keep only records of matching records, keep left geometry column


`Geopandas`では３種類の「結合方法」があります。

- left（左外部結合） ：左側のDataFrameはすべてに残す、マッチしない項目は空とする、左側のDataFrameのgeometryを残す
- right（右外部結合） ：右側のDataFrameはすべてに残す、マッチしない項目は空とする、右側のDataFrameのgeometryを残す
- inner（内部結合） ：マッチした行のみ残す、左側のDataFrameのgeometryを残す

:::note info
Tip

> The PyGIS book provides a great overview of spatial predicates and join-types, complete with explanatory drawings.

[PyGIS book](https://pygis.io/docs/e_spatial_joins.html)に、空間結合の「述語」と「結合条件」について、図を用いたすばらしい説明があります。
:::

# Load input data

> As a practical example, let’s find the population density at each of the addresses from earlier in this lesson, by combining the data set with data from a population grid.

> The population grid data is available from HSY, the Helsinki Region Environmental Services, for instance via their WFS endpoint or opendata.

例として、[以前のレッスン](https://qiita.com/t_t_japan/items/8e2fac435aa319bf43e7)で使った住所情報から人口密度を検出しましょう。データセットを人口グリッドのデータと組み合わせます。

人口グリッドは、 [HSY, the Helsinki Region Environmental Services](https://www.hsy.fi/en/environmental-information/open-data/)にあるものが使えます。例えばWFSエンドポイントやオープンデータポータルを通じて取得できます。 

まずは、事前準備ならびに、住所情報の読み込みです。
```
import pathlib
import geopandas


NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"

addresses = geopandas.read_file(DATA_DIRECTORY / "addresses.gpkg")

addresses.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/aeb16c8b-521c-4640-bd8c-14ac24825e77.png)


次に、WFSエンドポイントから人口グリッドを取得します。
WFSエンドポイント経由の場合は、crsの指定が必要なようです。
```
population_grid = geopandas.read_file(
    (
        "https://kartta.hsy.fi/geoserver/wfs"
        "?service=wfs"
        "&version=2.0.0"
        "&request=GetFeature"
        "&typeName=asuminen_ja_maankaytto:Vaestotietoruudukko_2020"
        "&srsName=EPSG:3879"
    ),
)
population_grid.crs = "EPSG:3879"  # for WFS data, the CRS needs to be specified manually
population_grid.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/ff6a648a-8d4e-4a5a-9b8e-ca1de69eff5e.png)

crsも見てみましょう。EPSG:3879になっていますね。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/808da448-0573-4774-a09b-f789fc7b3b40.png)

[オープンデータポータル](https://avoidatastr.blob.core.windows.net/avoindata/AvoinData/)からも取得してみましょう。
```
population_grid = geopandas.read_file(
    "https://avoidatastr.blob.core.windows.net/avoindata/AvoinData/"
    "6_Asuminen/Vaestotietoruudukko/Shp/Vaestotietoruudukko_2021_shp.zip"
)
population_grid.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/5a1cda3e-678a-4a45-b0ee-ba24d8f05750.png)

こちらの場合、WFSエンドポイント経由と全く同じというわけではないようです。こちらは列名が大文字です。また、gml_id列に相当する列も内容です。念のため、crsも見ましょう。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/724240d4-cf76-4300-bddb-cb68ff1907c0.png)

EPSG:3879ではありません。個別に定義したもののようです。まあ`name:EUREF-FIN / GK25FIN`とあるので、EPSG:3879相当なんだろうと解釈しておきます。ひとまず。


:::note info
Concatenating long strings

> In the WFS address above, we split a long string across multiple lines. Strings within parentheses are automatically concatenated (joined together) without needing any operator (e.g., +).

> For clarity, the example includes an additional set of parentheses, but the parentheses of the method call would suffice.

WFS アドレスでは、長い文字列を複数の行に分割します。括弧、”...”内の文字列は、例えば、+演算子などを使わなくても、自動的に連結（結合）されます。

今回はわかりやすくするために、この例では、"..."の外側に括弧、(...)がありますが、これは無くても成立します。言い換えると、以下でも問題なしです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/81d8ba44-7651-4104-843f-48e08c3c0a8c.png)


:::


> The population grid has many columns, and all of its column names are in Finnish. Let’s drop (delete) all of the columns except the population total, and rename the remaining to English:

人口グリッドはたくさんの列を持っています。またフィンランド語です。必要な列だけにして、列名を英語にしましょう。

**【注意】以下のコードは、オープンデータポータルから取得した場合の想定です。WFSエンドポイントから取得した場合は列抽出処理でエラーになります。前述の通り、列名が違うので。**

```
population_grid = population_grid[["ASUKKAITA", "geometry"]]
population_grid = population_grid.rename(columns={"ASUKKAITA": "population"})
population_grid.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/3a4988c3-d386-4cf8-a811-44a844329e1c.png)



> Finally, calculate the population density by dividing the number of inhabitants of each grid cell by its area in km²:

最後に人口密度（人/km²）を求めます。人数をgeometryの`area`（**面積**）で割ります。

```
population_grid["population_density"] = (
    population_grid["population"]
    / (population_grid.area / 1_000_000)
)
population_grid.head()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/4bfd7df7-27d1-40f3-81da-8075bd32d778.png)

なお、人/m²だと値が小さすぎます。0.000128等になってしまいます。人/km²としているのはそのためだと思われます。
この場合、`population_grid.area`を1,000,000（1,000m×1,000m）で割る必要があります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4076503/8649a6c0-b72f-41b9-aaf9-d7576fc94df2.png)

:::note infor
Tip: Coding style - big numbers, operators in multi-line expressions

> When using very large numbers, such as 1 million to convert m² to km², you can use underscore characters (_) as thousands separators. The Python interpreter treats numbers with underscores as regular numeric values. This syntax also works for grouping numbers by other logics, like hexadecimal or binary values.

> For expressions like mathematical formulas that spread across multiple lines, it’s good coding style to place an operator at the beginning of a new line rather than trailing it on the previous line. This improves readability, as explained in the PEP-8 styling guidelines.

かなり大きい数字（例えば100万）を㎡やkm²に変換するときは、アンダースコア（"_"）を3桁区切りとして使います。Pythonのインタプリタは、アンダースコアで数字を通常の数値として扱います。二進数や十六進数のように。

[PEP 515 – Underscores in Numeric Literals](https://peps.python.org/pep-0515/)


複数の行にまたがる数式などの場合、演算子は行の最後ではなく、次の行の先頭におくことが良いコーディングスタイルです。PEP-8のスタイリングガイドラインで説明されているように、読みやすさが向上します。

[Should a Line Break Before or After a Binary Operator?](https://peps.python.org/pep-0008/#should-a-line-break-before-or-after-a-binary-operator)

:::


後半に続きます。