Db2地理空間分析: 空間処理変数を使ってみる

`Db2地理空間分析`シリーズ


当記事は「7. 空間処理変数を使ってみる」です。

1. [Db2の地理空間分析機能](https://qiita.com/nishikyon/items/2f8bc4858e2d54bf9b1a)
1. [地理空間分析機能を有効にする](https://qiita.com/nishikyon/items/4cd0430e4ebb7d541a82)
1. [JGD2011座標系の設定](https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
1. [SHAPEファイルのインポート](https://qiita.com/nishikyon/items/fd6d166fdbf6523930d4)
1. [緯度・経度情報をインポート](https://qiita.com/nishikyon/items/cf90796a13fa7d2cf769)
1. [DBeaverでDb2の地理空間情報を表示してみる](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5)
1. [**空間処理変数を使ってみる**](https://qiita.com/nishikyon/items/7eeb7ee541e7d4f76a7f)

# 1. 空間処理関数
[空間処理関数](https://www.ibm.com/docs/ja/db2/11.5?topic=analytics-spatial-functions)を使用して、形状間の距離やある形状が別の形状に含まれるかどうかの判別などが可能です。 これをSQLのWHERE句として使用することができます。

こちらページに一覧があります:
https://www.ibm.com/docs/ja/db2/11.5?topic=analytics-spatial-functions

例として2つの関数を紹介します。

## ST_Distance関数　　
使い方詳細： https://www.ibm.com/docs/ja/db2/11.5?topic=sf-st-distance

ST_Distance 関数は、2 つの形状およびオプションとして単位を入力パラメーターとして取り、1番目の 形状内の任意のポイントと 2番目の形状内の任意のポイントとの間の最短距離を戻します。

つまり、ある地点からxkm以内の場所のリストを出す、というようなSQLに使用できます。

JR東京駅(緯度35.6809591, 経度139.7673068)から距離1km以内のLOCATION(ST_Point属性の列)
```sql
-- 2011は自分で設定したJGD2011のSRS ID(参考　https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
WHERE ST_Distance(
ST_Point (139.7673068, 35.6809591, 2011), LOCATION, 'KILOMETER') <= 1
```

## ST_Within関数　　
使い方詳細： https://www.ibm.com/docs/ja/db2/11.5?topic=sf-st-within

ST_Within 関数は、ある形状が完全に別の形状内にあるかどうかを判別するのに使います。

つまり、市区町村境界などの囲まれた形状の中にある、場所のリストを出す、というようなSQLに使用できます。


中央区内のLOCATION(LOCATIONはST_Point属性の列、中央区の形状はSHAPEにあるとする)
```sql
WHERE ST_Within (LOCATION, SHAPE)
```
↑ちょっとイメージ湧かないかもしれないので、次のセクションで今までロードしたはずのデータでのSQL実行例を紹介します。

# 2. SQLと視覚化のサンプル
### 前提: データの準備
以下でDb2のテーブルに入れたものを使用します:
- [4. SHAPEファイルのインポート](https://qiita.com/nishikyon/items/fd6d166fdbf6523930d4)
- [5. 緯度・経度情報をインポート](https://qiita.com/nishikyon/items/cf90796a13fa7d2cf769)

## 2.1 JR東京駅(緯度35.6809591, 経度139.7673068)から距離1km以内ののシェアサイクルポートの情報
[ST_Distance関数](https://www.ibm.com/docs/ja/db2/11.5?topic=sf-st-distance)

```sql
-- 2011は自分で設定したJGD2011のSRS ID(参考　https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
ST_Distance(ST_Point (139.7673068, 35.6809591, 2011), ST.GEO_POINT, 'KILOMETER') <=1
```

で1km以内の条件となる

SQL:
```sql
-- 2011は自分で設定したJGD2011のSRS ID(参考　https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
SELECT ST.NAME, ST.CAPACITY, 
       ST_X(ST.GEO_POINT) LON, ST_Y(ST.GEO_POINT) LAT
FROM SHAREBIKE_STATION ST 
WHERE ST_Distance(ST_Point (139.7673068, 35.6809591,  2011), ST.GEO_POINT, 'KILOMETER') <=1
```

 結果： 21件
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/6771beaf-2c49-47f3-8faa-d801efed39a0.png)


ではこの結果をJupyter notebook + python + [folium]() + [Db2 Jupyter Notebook Extensions](https://github.com/IBM/db2-jupyter)を使って視覚化してみます。
 
### Jupyter notebookのコード
```python
# Db2 Jupyter Notebook Extensions のロードと実行
import os
path = "./db2.ipynb"
if not os.path.exists(path):
    !wget https://raw.githubusercontent.com/IBM/db2-jupyter/master/db2.ipynb
else:
    print("db2.ipynb is existing.")
%run db2.ipynb
```
```python
# Db2接続情報のセット
# 下記にDb2接続情報をセットしてください
user = "xxxx" # db2ユーザーid
password = "xxxx" #db2パスワードを入れる
databese = "xxxx" # db2 Db名
host =  "xxxx" # db2 ホスト名
port =  9999 #Db2ポート番号
```


```python
# Db2 Jupyter Notebook Extensionsを使用してDb2に接続（SSL接続）
%sql CONNECT TO {databese} USER {user} USING {password} HOST {host} port {port} SSL TRUE
```

```python
# JR東京駅(緯度35.6809591, 経度139.7673068)から距離1km以内ののシェアサイクルポート
# 2011は自分で設定したJGD2011のSRS ID(参考　https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
df =%sql \
SELECT ST.NAME, ST.CAPACITY, \
       ST_X(ST.GEO_POINT) LON, ST_Y(ST.GEO_POINT) LAT \
FROM SHAREBIKE_STATION ST  \
WHERE ST_Distance(ST_Point (139.7673068, 35.6809591,  2011), ST.GEO_POINT, 'KILOMETER') <=1
df 
```
```python
import folium
# 地図生成
folium_map = folium.Map(location=[35.6809591, 139.7673068], zoom_start=14, height = 500)

# マーカープロット　
for i, row in df.iterrows():
    folium.CircleMarker(
        location=[row['LAT'], row['LON']],
        tooltip=row['NAME'],
        radius=5,
        color='green',
        fill_color='green'
    ).add_to(folium_map)

#東京駅から1km 円と東京駅地点プロット
folium.Circle(
        [35.6809591, 139.7673068], #東京駅
        radius=1000,
    ).add_to(folium_map)
folium.Marker(
        location=[35.6809591, 139.7673068],
        tooltip='東京駅'
    ).add_to(folium_map)

folium_map
```
こんな表示ができます！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/f54b91e8-3040-3f37-6b4b-956002c9b8fe.png)



### DBeaverで表示

[6. DBeaverでDb2の地理空間情報を表示してみる](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5)で解説したDBeaverで表示するSQLは以下になります:

```sql
-- 2011は自分で設定したJGD2011のSRS ID(参考　https://qiita.com/nishikyon/items/fc97391f02b4e9e4bb9c)
SELECT ST.NAME, ST.CAPACITY, 
        CAST( ST_AsText(ST.GEO_POINT)AS VARCHAR(50) ) AS LOC
FROM SHAREBIKE_STATION ST 
WHERE ST_Distance(ST_Point (139.7673068, 35.6809591, 2011), ST.GEO_POINT, 'KILOMETER') <=1
```

結果が表示されたら[4.2 結果パネルに表示された空間情報列のformatを変更](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5#42-%E7%B5%90%E6%9E%9C%E3%83%91%E3%83%8D%E3%83%AB%E3%81%AB%E8%A1%A8%E7%A4%BA%E3%81%95%E3%82%8C%E3%81%9F%E7%A9%BA%E9%96%93%E6%83%85%E5%A0%B1%E5%88%97%E3%81%AEformat%E3%82%92%E5%A4%89%E6%9B%B4)を実施後、 [5. Spatialビューの表示](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5#5-spatial%E3%83%93%E3%83%A5%E3%83%BC%E3%81%AE%E8%A1%A8%E7%A4%BA)を行なってください。

以下のように表示されます（1km円は表示されません、点のみ）：
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/8b6ee512-3bd9-70a5-fd6d-97eb5a18fe3f.png)


## 2.2 千代田区と中央区と港区のシェアサイクルポートの情報
[ST_Within関数](https://www.ibm.com/docs/ja/db2/11.5?topic=sf-st-within)

よくばり（？）にも、3つのエリア情報(Polygon)が必要なので、まずは3つのエリア千代田区と中央区と港区をSELECTするSQLは以下です：
```sql
SELECT GEO
FROM CITY_SHAPE 
WHERE N03_004 IN ('中央区', '千代田区', '港区') # N03_004には市区町村名が入っている
```
この結果を[ST_Within関数](https://www.ibm.com/docs/ja/db2/11.5?topic=sf-st-within)で使います。最終的なSQL全体は以下のようになります:
```sql
with AREA AS (
    SELECT GEO
    FROM CITY_SHAPE 
    WHERE N03_004 IN ('中央区', '千代田区', '港区')
)
SELECT ST.NAME, ST.CAPACITY, ST_X(ST.GEO_POINT) LON, ST_Y(ST.GEO_POINT) LAT
FROM SHAREBIKE_STATION ST, AREA
WHERE ST_Within(ST.GEO_POINT, AREA.GEO) 
--- ↑ここで'中央区', '千代田区', '港区'のPolygonに内包されたシェアサイクルポートのPoint(経度、緯度)という条件にしている

```
結果：　331件
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/5ccab7ff-a5de-2f8b-436a-2d8b9f2c0961.png)

ではこの結果をJupyter notebook + python + [folium]() + [Db2 Jupyter Notebook Extensions](https://github.com/IBM/db2-jupyter)を使って視覚化してみます。

### Jupyter notebookのコード
```python
# Db2 Jupyter Notebook Extensions のロードと実行
import os
path = "./db2.ipynb"
if not os.path.exists(path):
    !wget https://raw.githubusercontent.com/IBM/db2-jupyter/master/db2.ipynb
else:
    print("db2.ipynb is existing.")
%run db2.ipynb
```
```python
# Db2接続情報のセット
# 下記にDb2接続情報をセットしてください
user = "xxxx" # db2ユーザーid
password = "xxxx" #db2パスワードを入れる
databese = "xxxx" # db2 Db名
host =  "xxxx" # db2 ホスト名
port =  9999 #Db2ポート番号
```


```python
# Db2 Jupyter Notebook Extensionsを使用してDb2に接続（SSL接続）
%sql CONNECT TO {databese} USER {user} USING {password} HOST {host} port {port} SSL TRUE
```
```python
# 千代田区と中央区と港区のシェアサイクルポートの情報
df = %sql \
    with AREA AS ( \
    SELECT GEO \
    FROM CITY_SHAPE  \
    WHERE N03_004 IN ('中央区', '千代田区', '港区')\
    )\
    SELECT ST.NAME, ST.CAPACITY, ST_X(ST.GEO_POINT) LON, ST_Y(ST.GEO_POINT) LAT \
    FROM SHAREBIKE_STATION ST, AREA \
    WHERE ST_Within(ST.GEO_POINT, AREA.GEO)
display(df)
```

```python
# '中央区', '千代田区', '港区'の境界線情報(Polygon)取得
df_area =%sql \
    SELECT OBJECTID ,N03_001 ,N03_002 ,N03_003, N03_004, N03_007,\
    CAST(ST_AsText(GEO) AS VARCHAR(30000)) LOC  \
    FROM CITY_SHAPE WHERE N03_004 IN ('中央区', '千代田区', '港区')

#取得情報の表示
display(df_area)
```
```python
import folium
import geopandas as gpd
import 
# 地図生成
folium_map = folium.Map(location=[35.6636391,139.7606268], zoom_start=12, height = 800)

#境界線の描画 geopandas　GeoDataFrameを使用
gdf = gpd.GeoDataFrame(df_area, geometry = gpd.GeoSeries.from_wkt(df_area.LOC.astype('str')))　
geo_j = gdf['geometry'].to_json()
folium.GeoJson(geo_j, style_function=lambda x: {'fillColor': 'orange'}).add_to(folium_map)
```

```python
#シェアサイクルポートマーカープロット＋地図表示

for i, row in df.iterrows():
    folium.CircleMarker(
        location=[row['LAT'], row['LON']],
        tooltip=row['NAME'],
        radius=5,
        color='green',
        fill_color='green'
    ).add_to(folium_map)

folium_map
```
こんな表示ができます！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/72f90d42-fd8a-c926-21ee-e94be0f1ea64.png)

### DBeaverで表示

[6. DBeaverでDb2の地理空間情報を表示してみる](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5)で解説したDBeaverで表示するSQLは以下になります:
```sql
with AREA AS (
    SELECT N03_004, GEO
    FROM  CITY_SHAPE 
    WHERE N03_004 IN ('中央区', '千代田区', '港区')
)
SELECT ST.NAME, ST.CAPACITY, CAST( ST_AsText(ST.GEO_POINT)AS VARCHAR(50) ) AS LOC, 
AREA.N03_004, CAST( ST_AsText(AREA.GEO)AS VARCHAR(30000) )AS AREA
FROM  SHAREBIKE_STATION ST, AREA
WHERE ST_Within(ST.GEO_POINT, AREA.GEO) ;
```

結果が表示されたら[4.2 結果パネルに表示された空間情報列のformatを変更](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5#42-%E7%B5%90%E6%9E%9C%E3%83%91%E3%83%8D%E3%83%AB%E3%81%AB%E8%A1%A8%E7%A4%BA%E3%81%95%E3%82%8C%E3%81%9F%E7%A9%BA%E9%96%93%E6%83%85%E5%A0%B1%E5%88%97%E3%81%AEformat%E3%82%92%E5%A4%89%E6%9B%B4)を実施後、 [5. Spatialビューの表示](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5#5-spatial%E3%83%93%E3%83%A5%E3%83%BC%E3%81%AE%E8%A1%A8%E7%A4%BA)を行なってください。

以下のように表示されます（AREAは'中央区', '千代田区', '港区'3レコード分のみ選択）：
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294516/358ac8ef-f877-780f-319d-1171d1106882.png)

# 3. まとめ
いかがでしたか？　シリーズ7本目にしてやっとDb2地理空間分析としてSQLの結果を地図にマッピングして眺められた感があると思います(!)

Db2での地理空間分析の最大の利点は、Db2にデータを持っている場合、いつも使っているテーブルのデータと位置情報を組み合わせた分析がDb2にデータを置いたままできることです。

ぜひみなさんのDb2データに空間情報を付加して分析してみてください。

なお、`Db2地理空間分析`シリーズ

- 　4. [SHAPEファイルのインポート (ファイル作成部分のみ)](https://qiita.com/nishikyon/items/fd6d166fdbf6523930d4)
- 　5. [緯度・経度情報をインポート](https://qiita.com/nishikyon/items/cf90796a13fa7d2cf769)
- 　6. [DBeaverでDb2の地理空間情報を表示してみる](https://qiita.com/nishikyon/items/60e23324b1ec1f65deb5)
- 　7. [空間処理変数を使ってみる](https://qiita.com/nishikyon/items/7eeb7ee541e7d4f76a7f)

の内容はIBM Database Dojoとして 
[Db2でシェアサイクルポートの地理情報分析をやってみよう!](https://ibm-developer.connpass.com/event/259583/)
というハンズオンを実施しています。
ハンズオンで使用したJupyter Notebookを公開していますので、ぜひ参考にしてみてください。
- セッション資料: https://speakerdeck.com/kyokonishito/20220928-db2-spatial
- セッション動画: https://video.ibm.com/recorded/132165139
- セッションで使用したJupyter Notebook: https://github.com/kyokonishito/Db2_Spatial_Analysis

以上です。