インデックスフィーチャーの生成による空間統計

# インデックフィーチャー（格子状）の空間分析を全部SQLで実現

　格子状の空間に区切る事により、定量的に空間を分析することが可能です。今回は、空間を等間隔に分割する格子をデータベース内部に作成し、データ分析を行います。

## 今回のゴール
PostgreSQLで格子状の空間分析を行い、その結果をCesiumで３D表示します。
<img width="463" alt="スクリーンショット 2017-12-07 9.28.08.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/2bfd8ac7-0587-6601-fcc3-19ec2589598d.png">
公開先 http://jsdo.it/tkama/I4qt


## 分析対象 
　OpenStreetMap(OSM)の[生活道路（Living Street）](http://wiki.openstreetmap.org/wiki/JA:Tag:highway%3Dliving_street)をカウントします。その理由として、OSMのwikiには下記のような記載があります。
> 日本の法律では、Living Streetに該当する道路区分が定められていません。日本国内ではhighway=residentialを使用してください。

日本にどの程度、Living Streetが存在するのか、格子状に区切って集計してみます。（今回は関東限定ですが。。。）


# 分析準備
## 1.データの入手
 [GEOFABLIK](http://download.geofabrik.de/asia/japan.html)から、エリア毎にOSMデータをダウンロードします。今回は手元のノートPCの空き容量が少ないので、関東を選択しました。

![図1.png](https://qiita-image-store.s3.amazonaws.com/0/53918/488d085f-9a54-e808-8944-ffafba8c929b.png)

## 2.データベースの構築
　ダウンロードしたら、解凍しQGISで表示し、DBツールでPostGISに取り込みます。
1. gis.osm_roads_free_1.shpを選択してQGISで開きます
<img width="219" alt="スクリーンショット 2017-12-06 23.41.47.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/be490ca2-7738-7e3c-8089-25a171970f42.png">
2. 読み込むと地図に表示されます
<img width="614" alt="スクリーンショット 2017-12-06 23.47.19.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/be2696fa-a6b4-5f13-2a64-2c571ec52cbf.png">
[ © OpenStreetMap contributors](http://www.openstreetmap.org/copyright)
3. DBマネージャの起動
<img width="614" alt="スクリーンショット 2017-12-07 23.01.37.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/9950a0a9-795b-de2b-7f9e-791abbf83e96.png">

4. ベクタレイヤーのインポートツールの起動
![スクリーンショット 2017-12-06 23.57.55.png](https://qiita-image-store.s3.amazonaws.com/0/53918/8eeadbca-6e08-afa4-4267-bb537dff589d.png)
5. データベースに投入

 - 入力レイヤの選択：gis.osm_roads_free_1
 - スキーマの選択： スキーマの意味がわからない人はpublic
 - テーブル名: kanto_osm_roads

## 3. インデックフィーチャーの作成
　日本には伝統的なインデックフィーチャーである[標準地域メッシュ](https://ja.wikipedia.org/wiki/%E5%9C%B0%E5%9F%9F%E3%83%A1%E3%83%83%E3%82%B7%E3%83%A5)があります。今回は約10km四方である2次メッシュを作成します。
※generate_series等、PostgreSQLの方言がありますのでご注意下さい。(sqlite3では動作しないです）

```sql:CREATE_TABLE_JMESH_L2.sql
DROP TABLE IF EXISTS public.jmesh_l2;
CREATE TABLE public.jmesh_l2 AS
SELECT 
  TRUNC( y*1.5 ) *10000  + 
  TRUNC(x - 100)*100 + 
  TRUNC(TRUNC( (((y*60)::int)%40 )/5)*10) +
  TRUNC (( x -TRUNC(x))/(45.0/360.0) )
 AS meshcode 
 ,  ST_MakeEnvelope( x , y , x+(45.0/360.0) ,  y+(5.0/60.0) , 4326 )  AS the_geom
FROM 
 generate_series( 20.000000, 46.0 - (5.0/60.0) , 5.0/60.0) AS y
LEFT OUTER JOIN 
 generate_series( 122.000000, 154 -(45.0/360.0), (45.0/360.0) ) AS x  
ON y > 0
 ;

--空間インデックスの作成
CREATE INDEX jmesh_l2_gidx
  ON public.jmesh_l2
  USING gist
  (the_geom);

```
ささっと作ったので、色々問題があるかもしれません。
バグを見つけたら、優しく指摘頂けると助かります。

　上記のSQLを実行すると、79872行のメッシュがテーブル内に格納されます。QGISで表示すると下図のようになります。（目がチカチカします。）

<img width="609" alt="スクリーンショット 2017-12-06 23.21.35.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/6a55739c-a9bb-ca4d-4700-a00128193720.png">

東京エリアを拡大
<img width="747" alt="スクリーンショット 2017-12-06 23.34.43.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/7884ea9a-bfdb-573c-64ee-b6938346a1e7.png">

# 集計の実施

## 基本統計
道路種別毎の基本的な統計を取得します。

```sql:基本統計.sql
SELECT 
 fclass 
 , COUNT( * ) AS num
 , SUM( ST_Length( geom::geography ) )::int AS sum
 , TRUNC( MIN( ST_Length( geom::geography ) ) :: numeric, 2) AS min
 , MAX( ST_Length( geom::geography ) )::int AS max
 , AVG( ST_Length( geom::geography ) )::int AS avg
FROM public.kanto_osm_roads
GROUP BY 1
ORDER BY 2 DESC 
LIMIT 15 ;

```
| fclass | num | sum | min | max | avg |
|:------:|:---:|:---:|:---:|:--:|:----:|
|residential|544747|83669427|0.00|12365|154|
|unclassified|475896|84770621|0.00|14215|178|
|service|104568|12795888|0.13|12629|122|
|tertiary|72046|26301809|0.19|14751|365|
|footway|71160|7149357|0.21|5185|100|
|path|62031|18269730|0.08|16812|295|
|track|59694|15389094|0.78|9963|258|
|steps|19175|406398|0.21|502|21|
|primary|14917|9768399|0.71|23138|655|
|secondary|13443|9128186|0.77|17720|679|
|living_street|13205|943751|2.31|4378|71|
|trunk|11649|7676985|0.04|16794|659|
|pedestrian|7498|1005841|0.30|4307|134|
|motorway|6026|3726151|9.43|15534|618|
|track_grade1|5277|1725858|1.83|7966|327|

今回注目している、living_street は１１位。
関東には、13205本あることがわかります。

## 空間統計
ポイント
 `ST_Length( ST_Intersection( r.geom , m.the_geom )::geography )`

ST_Intersection( 空間形状１ , 空間形状２ )関数は２つの引数で指定した空間形状同士の重なり合う部分だけを抽出しています。

　今回は、OSMの道路の線形状と2次メッシュの四角形との重なり合った部分を切り出します。結果は、線形状が返却されます。
また、WHERE句でliving_streetだけを抽出しています。

```sql:空間統計.sql
SELECT
   meshcode 
 , TRUNC( SUM(  ST_Length( ST_Intersection( r.geom , m.the_geom )::geography ) ) )
 , COUNT( * )
 , ST_AsText( m.the_geom ) AS wkt
FROM 
      public.kanto_osm_roads AS r
    , public. jmesh_l2 AS m
WHERE 
 ST_Intersects( r.geom , m.the_geom )
AND
 fclass = 'living_street'
GROUP BY 1 , the_geom;
```
手元の環境、３年目のmac(MacBook Pro (Retina, 13-inch, Mid 2014))とPostgreSQL9.6 で実行時間は２秒弱で返却されます。

 ST_AsTextで格子状の形状（２次メッシュ）を出力したので、QGISで表示してみます。着色は、リンク長の合計で長いほど赤くしています。

![3.png](https://qiita-image-store.s3.amazonaws.com/0/53918/42ddd7cc-d6cf-536d-df29-c7be0a8a4ef2.png)

画像を見れば、どこにliving_streetのリンクが多いのか一目瞭然ですね。

## Cesium表現用にCZMLへ出力

Webブラウザで３D表現が可能な地図ライブラリである、Cesiumで３D棒グラフとして描画を試みます。
Cesium用のJSONフォーマットであるCZMLに変換してみます。

```sql:CZML変換.sql
WITH ans AS (
		SELECT
		   meshcode 
		 , TRUNC( SUM(  ST_Length( ST_Intersection( r.geom , m.the_geom )::geography ) ) ) AS length
		 , COUNT( * ) AS c
		 , the_geom 
		FROM 
		      public.kanto_osm_roads AS r
		    , public.mesh2 AS m
		WHERE 
		 ST_Intersects( r.geom , m.the_geom )
		AND
		 fclass = 'living_street'
		GROUP BY 1 , the_geom 
		ORDER BY 2 DESC 
	)
	
	SELECT
	array_to_string(
		array_agg(
	 '{
	    "id" : "' || meshcode || '",
	    "name" : "' || meshcode || '",
	    "rectangle" : {
		"coordinates" : {
		    "wsenDegrees" : [' 
			     || ST_XMIN(the_geom ) 
		    || ', ' || ST_YMIN( the_geom ) 
		    || ', ' || ST_XMAX( the_geom ) 
		    || ', ' || ST_YMAX( the_geom ) ||']
		},
		"height" : ' || length || ',
		"extrudedHeight" : 0,
		"fill" : true,
		"material" : {
		    "solidColor" : {
			"color": {
			    "rgba" : ['||   TRUNC( length /clevel.l256 ) || ', 0 , '||  255 - TRUNC( length /clebel.l256 ) || ', 210 ]
			}
		    }
		},
		"outline" : true,
		"outlineColor" : {
		    "rgba" : ['||   TRUNC( length /clevel.l256 ) || ', 0, '|| 255 - TRUNC( length /clebel.l256 ) || ' , 255 ]
		},
		"outlineWidth" : 1
	    }
	}'
		), ',')
	 
	 AS packet
	FROM ans , 
	(SELECT MAX(length / 256)+1 AS l256   FROM ans )  AS clevel 
```

## HTMLへCZMLの埋め込み

下のHTMLのコードの中の`/* ココにSQLの実行結果を貼り付ける */`にSQLの出力結果を貼り付けます。
実際に動作している[Webサイト](http://jsdo.it/tkama/I4qt)がありますので、参考にして下さい。

```html
<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8" />
<title>Cesium CZML rectangle 3D Japan mesh level2</title>
<meta name="viewport" content="width=device-width, user-scalable=no, initial-scale=1, maximum-scale=1">
<style type="text/css">@import url("https://cesiumjs.org/Cesium/Build/Cesium/Widgets/widgets.css");
body {
    padding: 0;
    margin: 0;
    overflow: hidden;
}
#cesiumContainer {
    position: absolute;
    top: 0;
    left: 0;
    width: 100%;
    height: 100%;
    margin: 0;
    padding: 0;
    overflow: hidden;
}
</style></head>
<body>
<script src="https://cesiumjs.org/Cesium/Build/Cesium/Cesium.js"></script>
<div id="cesiumContainer"></div>
<div id="loadingOverlay"><h1>Loading...</h1></div>
<div id="toolbar"></div>
<script type="text/javascript">
var t1 = [{
    "id" : "document",
    "name" : "",
    "version" : "1.0"
},
 /*
 ココにSQLの実行結果を貼り付ける 
 */
];

//デフォルトカメラ位置の登録
var extent = Cesium.Rectangle.fromDegrees(122, 20, 154, 46);
Cesium.Camera.DEFAULT_VIEW_RECTANGLE = extent;
Cesium.Camera.DEFAULT_VIEW_FACTOR = 0;

//Cesium地図ビューワの作成
var viewer = new Cesium.Viewer('cesiumContainer');
viewer.dataSources.add(Cesium.CzmlDataSource.load(t1));
</script>
</body>
</html>
```

# まとめ
 データベースの中だけで、インデックフィーチャーを作成し空間統計を行いました。
空間統計では、１０００万行のOSMのテーブルと8万行弱の２次メッシュのテーブルとExcelでは中々ハンドリングが難しいデータ量を扱いましたが、最終的なSQLの実行時間は２秒程度とスピーディーな集計が実現できました。
　in-Databaseの強みは、大きなデータの出し入れが少ないため、データの集計に集中することも強みかと思われます。

# 出典
[地理院タイル 白地図](https://maps.gsi.go.jp/development/ichiran.html)
[© OpenStreetMap contributors](http://www.openstreetmap.org/copyright)