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

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

今回のゴール

PostgreSQLで格子状の空間分析を行い、その結果をCesiumで３D表示します。

公開先 http://jsdo.it/tkama/I4qt

分析対象

　OpenStreetMap(OSM)の生活道路（Living Street）をカウントします。その理由として、OSMのwikiには下記のような記載があります。

日本の法律では、Living Streetに該当する道路区分が定められていません。日本国内ではhighway=residentialを使用してください。

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

分析準備

1.データの入手

GEOFABLIKから、エリア毎にOSMデータをダウンロードします。今回は手元のノートPCの空き容量が少ないので、関東を選択しました。

2.データベースの構築

　ダウンロードしたら、解凍しQGISで表示し、DBツールでPostGISに取り込みます。

gis.osm_roads_free_1.shpを選択してQGISで開きます

2. 読み込むと地図に表示されます

[ © OpenStreetMap contributors](http://www.openstreetmap.org/copyright) 3. DBマネージャの起動

ベクタレイヤーのインポートツールの起動
データベースに投入

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

3. インデックフィーチャーの作成

　日本には伝統的なインデックフィーチャーである標準地域メッシュがあります。今回は約10km四方である2次メッシュを作成します。
※generate_series等、PostgreSQLの方言がありますのでご注意下さい。(sqlite3では動作しないです）

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で表示すると下図のようになります。（目がチカチカします。）

東京エリアを拡大

集計の実施

基本統計

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

基本統計.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

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で表示してみます。着色は、リンク長の合計で長いほど赤くしています。

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

Cesium表現用にCZMLへ出力

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

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サイトがありますので、参考にして下さい。

<!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の強みは、大きなデータの出し入れが少ないため、データの集計に集中することも強みかと思われます。

出典

地理院タイル白地図
 © OpenStreetMap contributors