Go to Qiita Advent Calendar 2024 Top
LoginSignup
7
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

FOSS4G Advent Calendar 2017Day 7
@takahi

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

Last updated at Posted at 2017-12-07

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

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

今回のゴール

PostgreSQLで格子状の空間分析を行い、その結果をCesiumで３D表示します。
スクリーンショット 2017-12-07 9.28.08.png
公開先 http://jsdo.it/tkama/I4qt

分析対象

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

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

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

分析準備

1.データの入手

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

図1.png

2.データベースの構築

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

  1. gis.osm_roads_free_1.shpを選択してQGISで開きます
スクリーンショット 2017-12-06 23.41.47.png 2. 読み込むと地図に表示されます スクリーンショット 2017-12-06 23.47.19.png [ © OpenStreetMap contributors](http://www.openstreetmap.org/copyright) 3. DBマネージャの起動 スクリーンショット 2017-12-07 23.01.37.png
  1. ベクタレイヤーのインポートツールの起動
    スクリーンショット 2017-12-06 23.57.55.png
  2. データベースに投入
  • 入力レイヤの選択：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で表示すると下図のようになります。（目がチカチカします。）

スクリーンショット 2017-12-06 23.21.35.png

東京エリアを拡大
スクリーンショット 2017-12-06 23.34.43.png

集計の実施

基本統計

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

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

3.png

画像を見れば、どこに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

7
8
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?