PostGISで空間クラスタリング関数 ST_ClusterKMeansの試用

# PostGISに空間クラスタリング関数の登場
　PostGISのバージョン2.3より[ST_ClusterKMeans](https://www.finds.jp/docs/pgisman/2.3.0/ST_ClusterKMeans.html)と[ST_ClusterDBSCAN](https://www.finds.jp/docs/pgisman/2.3.0/ST_ClusterDBSCAN.html)の2つの空間クラスタリング関数が追加されました。

　日本語のドキュメントにつきましては、[農研機構](http://www.finds.jp/)様のサイトにはバージョン[2.3](http://www.finds.jp/docs/pgisman/2.3.0/postgis-ja.html)はもちろん、[2.5](http://www.finds.jp/docs/pgisman/2.5.0/postgis-ja.html)のドキュメントが公開されております。この場を借りて厚く御礼申し上げます。

## データ準備
国土数値情報ダウンロードサービスから、[国土数値情報　避難施設データ](http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P20.html)を用いてみました。
災害対策基本法で定められた避難所の場所の地点を記録したデータで、埼玉県には3691個あります。

![base.png](https://qiita-image-store.s3.amazonaws.com/0/53918/e4ff9e59-3335-5bdc-5d4f-4e49032232e8.png)

# クラスタリング
 k-measn法のアルゴリズムついては「k個のクラスタに分割する」というルールを予め指定する必要があります。

## ST_ClusterKMeans関数
 この関数では実行結果として0〜(k-1)までのクラスタ番号が返却されます。またWindow関数であるためOVER句を付ける必要があります。

> integer ST_ClusterKMeans(geometry winset geom, integer number_of_clusters);

戻り値： k-means法のクラスタ番号（INTEGER型）
第1引数：地物（geometry型 ※geography型非対応）
第2引数：クラスタ数（INTEGER型）


### 3分割
k=3として埼玉県を3つのクラスタに分割してみます。

```sql:3-means
SELECT 
	  ST_ClusterKMeans(the_geom,3 )OVER() AS cluster_id 
	, 名称
	, the_geom
FROM shelter_saitama;

```

表示してみましょう。
![3means.png](https://qiita-image-store.s3.amazonaws.com/0/53918/7c214a7f-af2f-c580-a198-a2b3ef95b878.png)

ザッと埼玉県が3分割されました。

### 10分割
k=10として埼玉県を10個のクラスタに分割してみます。

```sql:10-means
SELECT 
	  ST_ClusterKMeans(the_geom,10 )OVER() AS cluster_id 
	, 名称
	, the_geom
FROM shelter_saitama;
```
![10.png](https://qiita-image-store.s3.amazonaws.com/0/53918/7a8c2e15-f1b3-f515-97de-87d4977332df.png)

ボロノイ図に見えてきました。。。

### 行政区域毎に2分割する
Window関数の練習として、行政区域毎にそれぞれ2つのクラスタに分割します。

```sql:city-2-means
SELECT 
	  ST_ClusterKMeans(the_geom,2 )over(PARTITION BY 行政区域) AS cluster_id 
	, 行政区域 
	, 名称
	, ST_AsText(the_geom)
FROM shelter_saitama
```
　OVERE句には、カッコの中にPARTITION BY句 とORDER BY句が記述できます。
　PARTITION BYでは集計対象のグループを定義する値が入る列を列挙することができ、GROUP BYと非常によく似ています。WINDOW関数は、SELECTの結果自体を集約しないという違いがあります。

![city_2.png](https://qiita-image-store.s3.amazonaws.com/0/53918/8d3ed823-d8d5-ed29-5bad-2cfbf10f5533.png)

## クラスタ結果の空間的表現
クラスタをより直感的に地図表現してみましょう。

### 各クラスタの重心にクラスタの所属地点数を取得

```sql:10-means-centroid-count
WITH t1 AS(
	SELECT 
		  ST_ClusterKMeans(the_geom,10 )OVER() AS cluster_id 
		, the_geom
	FROM shelter_saitama
)
SELECT 
 cluster_id
 ,ST_Centroid(ST_Union(the_geom)) AS 重心
 ,COUNT(cluster_id)
FROM t1
GROUP BY 1
ORDER BY 1
```
![centroid_count.png](https://qiita-image-store.s3.amazonaws.com/0/53918/9dce776f-c1a4-0bd2-002e-17f3af76d30f.png)

### 重心と各地点を線分で結ぶ

```sql:star
WITH 
--クラスタリング
t1 AS( 
	SELECT 
		  ST_ClusterKMeans(the_geom,10 )OVER() AS cluster_id 
		, 名称
		, the_geom
	FROM shelter_saitama
)
, 
--重心計算
t2 AS( 
SELECT 
 cluster_id
 , ST_Centroid(ST_Union(the_geom))
 ,COUNT(cluster_id)
FROM t1
GROUP BY 1
ORDER BY 1
)
--クラスタリング済みの点と重心点の線分作成
SELECT 
 t2.cluster_id 
 , count
 , ST_Union( ST_Makeline(ST_Centroid,the_geom))			 
FROM t1,t2
WHERE 
 t1.cluster_id = t2.cluster_id
GROUP BY 1,2
```
　WITH句を利用して、t1でクラスタ分類リストを作成し、そこから重心の取得結果をt2として保持し、更にt1とt2をJOINすることで重心と各点をグルーピングしています。

![star.png](https://qiita-image-store.s3.amazonaws.com/0/53918/8e409b5b-b1f2-87d7-4987-8d474127ddc1.png)



# 付録

## 環境構築
　PostGISはPostgreSQLの拡張機能です。インストールしてバージョンをチェックしてみましょう。今回の環境はAmazonのクラウドサービス、AWSの[RDS](https://aws.amazon.com/jp/rds/?sc_channel=PS&sc_campaign=acquisition_JP&sc_publisher=google&sc_medium=rds_rlsa_b&sc_content=rds_e&sc_detail=rds&sc_category=rds&sc_segment=278520267021&sc_matchtype=e&sc_country=JP&sc_brand=brand&ef_id=Cj0KCQiA_4jgBRDhARIsADezXci1LbMZgiRzgCjsND-jyKmVLW0uTedFDoPj2qpsKJfhgTq5SWIl7a4aAkzhEALw_wcB:G:s)から postgreSQL 11（評価版）を利用しています。

```sql:インストールとバージョンチェック
/**
 * PostGISのインストール
 **/
CREATE EXTENSION postgis;
-- クエリが 1 secs 399 msec で成功しました

/**
 * PostGISのバージョンチェック
 **/
SELECT PostGIS_Full_Version();
/* POSTGIS="2.5.0 r16836" 
[EXTENSION] 
PGSQL="110" 
GEOS="3.7.0-CAPI-1.11.0 673b9939" 
PROJ="Rel. 5.2.0, September 15th, 2018" 
GDAL="GDAL 2.3.1, released 2018/06/22" 
LIBXML="2.9.1" 
LIBJSON="0.12" 
RASTER
*/
```


　PostGISのバージョンが2.3以上であること、GEOSがインストールされていることを確認します。

## テーブル定義

```sql:createtable
CREATE TABLE public.shelter_saitama
(
    "行政区域" text,
    "名称" text,
    "住所" text,
    "施設の種類" text ,
    "収容人数" integer,
    "施設規模" integer,
    "災害分類_地震災害" boolean,
    "災害分類_津波災害" boolean,
    "災害分類_水害" boolean,
    "災害分類_火山災害" boolean,
    "災害分類_その他" boolean,
    "災害分類_指定なし" boolean,
    no integer ,
    the_geom geometry
)
WITH (
    OIDS = FALSE
)
;

```