polarsで地理データを扱う


Polars には plugin という、拡張機能を簡単に開発できるものがあります[^plugin]。

plugin のひとつとして、[polars_st](https://oreilles.github.io/polars-st/) という地理空間情報について扱えるものがあります。本記事では polars_st の簡単な使い方を紹介します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/484ae6a8-466f-0485-efca-b597acc87990.png)

https://oreilles.github.io/polars-st/


まだalpha版なのでちょっと使い勝手は良くないですが、早速使ってみました。

[^plugin]: pluginについてはこちら。  https://docs.pola.rs/api/python/stable/reference/plugins.html 少し前ですがいろんな便利ツールやプラグインについて紹介した記事があるので、そちらも一緒にどうぞ！  https://qiita.com/_jinta/items/e135c48676e08c993c92



# インストール

```
pip install polars-st
```

```python
import polars as pl
import polars_st as st
```


# データ読み込み
- GeoJSONを読みます
    - たぶん他の地理データファイル形式でも同様に行けるはず
- CRSの設定は適当にしてます

### 若干の修正
インストールしたままだとエラーを吐くので、一部修正してしまいます。

※なお、google colabでやってます

<details><summary>GeoJSON読むときにおこるエラーへの一次的な対処</summary>

```
%%writefile /usr/local/lib/python3.10/dist-packages/polars_st/utils/srid.py

from __future__ import annotations

from logging import warning

from polars_st._lib import get_crs_auth_code


def get_crs_srid_or_warn(crs: str) -> int | None:
    try:
        _auth, code = get_crs_auth_code(crs)
        if code.isdigit():
            return int(code, base=10)
        # warning.warn(
        #     f"Found an authority for {crs} but couldn't"
        #     f'convert code "{code}" to an integer srid. ',
        #     "The geometries SRID will be set to 0.",
        # )
    except ValueError:
        # warning.warn(
        #     f'Couldn\'t find a matching crs for "{crs}". The geometries SRID will be set to 0.'
        # )
        pass
    return None
```
</details>


### １．兵庫県のWiFiデータ[^hyogo]

```python
url_wifi = "https://www.geospatial.jp/ckan/dataset/06f57e8f-4abf-4d06-a499-62230a7968b5/resource/f4728608-efb7-4f72-a3d1-c90c74e0b431/download/lan.geojson"

# データ読み込み
df = st.read_file(url_wifi)

# 座標参照系の設定
df = df.with_columns(st.set_srid(srid=4326))
```

[^hyogo]: 兵庫県作成の公衆無線LAN（HYOGO Free Wi-Fi）の設置施設一覧を使用しています: https://www.geospatial.jp/ckan/dataset/lan/resource/f4728608-efb7-4f72-a3d1-c90c74e0b431

### ２．日本の都道府県の地図[^jmap]

```
!curl -L -o japan.geojson https://raw.githubusercontent.com/dataofjapan/land/master/japan.geojson
```

```python
df_ja = st.read_file("/content/japan.geojson")
df_ja = df_ja.with_columns(st.set_srid(srid=4326))
```

[^jmap]: 地球地図日本( http://www.gsi.go.jp/kankyochiri/gm_jpn.html ) を変換したデータを使用させていただいています: https://github.com/dataofjapan/land

# データ確認

`geometry` 列に地理情報が入る。binary で持つのでやや見にくい。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/db1104f2-77e3-bff1-4fe1-d16233c26ca4.png)

plot 機能もありそうでしたが、使い方分からなかったので geopandas に変換して描画しています。


<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/ec1e0273-a6a6-97f7-b0ee-f2c95b1182ea.png" width=50%>

重ねるとこんな感じ。

```python
fig, ax = plt.subplots(figsize=(6, 6))

df_ja.filter(pl.col("nam_ja")=="兵庫県").st.to_geopandas().plot(ax=ax, color="grey")
df.st.to_geopandas().plot(ax=ax);
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/f7ef046d-94ed-2983-8679-8e7b1c898d40.png" width=50%>





# geometryをいじる

- Point を中心にした半径の円を作成する：`st.buffer`

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/41aa0c83-8819-6d5e-fc97-9efff636f0ba.png" width=50%>

- 他のデータとのJOIN: `st.sjoin`
    - JOIN の条件はいろいろ指定できるけど、デフォルトは intersection（重なってるものをジョイン）
    - 結合後は geometry に相当するカラムが複数できるので、そこでミスが起きそう

```python
df_joined = df.st.sjoin(df_ja)
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/baf91a51-4139-846c-f3e0-cea9e53f7051.png" width=50%>


- 色々合わせてやってみる
    - 結合したジオデータフレームについて、
    - 右側（兵庫県の地図）に対して以下とのアンドのエリアを取得
        - WiFiの位置を中心に半径nで作った円
    - それをグラフにする

```python
(
    df_joined.select(
        pl.col("geometry_right").st.intersection(
            pl.col("geometry").st.buffer(0.1)
        )
    )
    .rename({'geometry_right': "geometry"})
    .st.to_geopandas().plot()
)
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/665792/cc7d44d9-8e6b-168a-e0c6-7ad61feec3ce.png" width=50%>


# その他

### チートシート
- sf（Rの同じ役割のライブラリ）のチートシートが分かりやすいです
- https://rstudio.github.io/cheatsheets/sf.pdf
- だいたい同じ処理が同じ名前で出来る気がします

### geopolarsとの関係
- [geopolars](https://geopolars.org/latest/) はもっと先鋭的に速くなるらしいです
- でもまだ開発は先になる模様
- https://github.com/oreilles/polars-st?tab=readme-ov-file#about-geopolars

# おわりに
まだドキュメントなど未整備な部分多そうですが、調べながらであれば使えそうです。引き続き追っていきたいと思います。ではでは！