PostGIS の ST_AsGeoJSON を使うと、空間情報を簡単に GeoJSON 形式で出力できます。しかし、空間情報に属性情報を付与したい場合もあります。人口などの量的情報はもちろんですが、Mapbox の Simplestyle に沿った描画情報 - An open platform > Styling features | Mapbox - が挙げられます。
PostgreSQL には row_to_json などの関数があり、SQL の結果を JSON で出力できます。これを利用すると、PostGIS に保存してある空間情報を GeoJSON で出力するときに属性を付与できます。バッチ処理でデータ形式を変換するときには便利なのでメモしておきます。
PostgreSQL: Documentation: 9.3: JSON Functions and Operators
データとしては、国土数値情報の将来推計人口メッシュを使います。
2010年の国勢調査等に基づき、2050年までの1km四方のメッシュ別の将来人口の試算を行ったものである.
データの準備
国土数値情報ダウンロードサービスの「データ形式」で「統一フォーマット(SHP・GML)」を選択し、「5. 各種統計」セクションに「将来人口推計メッシュ」があります。このチェックボックスを ON にして、ページ最下部の「選択」でダウンロードページに進みます。利用用途などの簡単なアンケートに答えて suikei140704-48.zip というファイルをダウンロードします。
ZIPファイルを展開すると Shape ファイルがあります。QGIS などのソフトウェアで簡単に表示を確認しておきましょう。各属性の意味は冒頭のリンク先に記載されています。文字エンコードは Shift-JIS です。
PostGISへの登録
Shape ファイルのデータを PostGIS に登録します。まずは "geotest" というデータベースを作成し、PostGIS 拡張を有効にします。PostgreSQL サーバーのホスト名、ポート番号、接続ユーザー情報などは適宜設定してください。
$ createdb geotest
$ cat << EOF | psql -d geotest
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
SELECT postgis_full_version();
EOF
次に、shp2pgsql で SQL の INSERT 文を作成してデータを登録します。先ほどダウンロードしたファイルは、ホームディレクトリ配下の ~/Downloads/suikei140704-48 に展開してあるものとしています (MacOSXのデフォルト設定)。
注意点として、Shapeファイルで値の無いフィールドは "-" になっていますので、これを NULL にしておきます。
$ shp2pgsql -W "Shift_JIS" ~/Downloads/suikei140704-48/suikei140704_00.shp |
sed "s/'-'/NULL/g" | psql geotest
"suikei140704_00" テーブルが作成されていますので、カラム名や件数などを確認しておきます。テーブル名は Shape ファイルのレイヤー名と一緒です。
$ psql geotest
> \d suikei140704_00
> SELECT count(*) FROM suikei140704_00;
> SELECT city_code / 1000 AS pref, count(*) FROM suikei140704_00
GROUP BY city_code / 1000 ORDER BY pref;
"city_code" は５桁の市区町村コードなので、先頭２桁が都道府県コードになります。テーブル定義を特に指定しない場合は整数型になりますので1000で除算して先頭２桁を取り出しています。文字列として扱うには、型をキャストするか別のテーブルにデータを流し直しましょう。
また、人が住んでいない（人口未登録）メッシュには市区町村コードが割り振られていないようです。NULL のカラムがあるレコードの扱いは事前に決めておきましょう。
GeoJSONで出力
こちらの記事を参考にして GeoJSON を組み立てていきます。
まずは "geom" 列を GeoJSON 形式にして geometry と名付けます。なお、データ確認用に出力は10行に制限しておきます。
SELECT ST_AsGeoJSON(geom)::json AS geometry
FROM "suikei140704_00" AS t
LIMIT 10;
次に、"pop2010" 列を population と名付けて JSON 形式に変換した列を properties とします。ここで row_to_json 関数を使います。
SELECT
ST_AsGeoJSON(geom)::json AS geometry,
row_to_json((
SELECT p FROM (
SELECT pop2010 AS population
) AS p
)) AS properties
FROM "suikei140704_00" AS t
LIMIT 10;
これらを更に row_to_json 関数で集約してひとつの列にします。"feature" という別名を付けておき、その別名に対して関数を適用させます。また、形状種別として "Feature" を固定で割り当てます。
SELECT row_to_json(feature)
FROM (
SELECT
'Feature' AS type,
ST_AsGeoJSON(geom)::json AS geometry,
row_to_json((
SELECT p FROM (
SELECT pop2010 AS population
) AS p
)) AS properties
FROM "suikei140704_00" AS t
) AS feature
LIMIT 10;
これで出力は GeoJSON っぽくなってきたと思います。
あとは列同士の関係を計算して属性 (properties) を追加します。たとえば、市区町村コードの付与、推計人口の差分計算などが挙げられます。ここでは推計人口の差分に応じて表示色を変更してみましょう。人口が増える見込みの場合は赤、減る場合は青、変わらない場合は緑、不明の場合は灰色にしておきます。
SELECT row_to_json(feature)
FROM (
SELECT
'Feature' AS type,
ST_AsGeoJSON(geom)::json AS geometry,
row_to_json((
SELECT p FROM (
SELECT
pop2010 AS population,
lpad(city_code::varchar, 5, '0') AS city_code,
pop2050 - pop2010 AS population_diff,
CASE
WHEN pop2050 - pop2010 > 0 THEN '#ff0000'
WHEN pop2050 - pop2010 < 0 THEN '#0000ff'
WHEN pop2050 - pop2010 = 0 THEN '#00ff00'
ELSE '#cccccc'
END AS "fill",
0.8 AS "fill-opacity",
CASE
WHEN pop2050 - pop2010 > 0 THEN '#ff0000'
WHEN pop2050 - pop2010 < 0 THEN '#0000ff'
WHEN pop2050 - pop2010 = 0 THEN '#00ff00'
ELSE '#cccccc'
END AS "stroke",
0.9 AS "stroke-opacity"
) AS p
)) AS properties
FROM "suikei140704_00" AS t
) AS feature
LIMIT 10;
最後に、都道府県ごとにひとつの GeoJSON にまとめます。複数の Feature を持つ GeoJSON になりますので、"FeatureCollection" でラッピングします。複数行のデータをひとつに集約するには array_to_json と array_agg を使います。都道府県コードは WHERE 句の部分で指定します。
最終的には以下のような SQL になります。なお psql の機能を使って都道府県コードを変数化しています。
SELECT row_to_json(featurecollection)
FROM (
SELECT
'FeatureCollection' AS type,
array_to_json(array_agg(feature)) AS features
FROM (
SELECT
'Feature' AS type,
ST_AsGeoJSON(geom)::json AS geometry,
row_to_json((
SELECT p FROM (
SELECT
pop2010 AS population,
lpad(city_code::varchar, 5, '0') AS city_code,
pop2050 - pop2010 AS population_diff,
CASE
WHEN pop2050 - pop2010 > 0 THEN '#ff3333'
WHEN pop2050 - pop2010 < 0 THEN '#3333ff'
WHEN pop2050 - pop2010 = 0 THEN '#33ff33'
ELSE '#cccccc'
END AS "fill",
0.6 AS "fill-opacity",
CASE
WHEN pop2050 - pop2010 > 0 THEN '#ff0000'
WHEN pop2050 - pop2010 < 0 THEN '#0000ff'
WHEN pop2050 - pop2010 = 0 THEN '#00ff00'
ELSE '#cccccc'
END AS "stroke",
0.8 AS "stroke-opacity"
) AS p
)) AS properties
FROM "suikei140704_00"
WHERE substr(lpad(city_code::varchar, 5, '0'), 1, 2) = :'prefcode'
) AS feature
) AS featurecollection;
都道府県コード "01" (北海道) のデータを出力する場合は以下のように実行できます。
$ psql -t -f export-as-geojson.sql -v prefcode="01" geotest > population-01.geojson
まとめて実行するには以下のようになります。SQLは複雑になりますが、別途プログラムを記述しなくても GeoJSON を生成できることが分かりました。
$ for i in `seq -w 47`; do
psql -t -f export-as-geojson.sql -v prefcode=$i geotest > population-$i.geojson
done
あとは geojson.io などで表示を確認しましょう。GitHub の地図も Mapbox なので、Gist を作成しても確認できます。自前で描画したい場合は Leaflet や D3 を使って JavaScript を書けば良いでしょう。
試しに東京のデータを Gist に貼付けると、きちんと地図上に表示されます。 - tokyo-population.geojson
東側と北側では人口が増える予測、それ以外の大部分では人口が減少する見込みであることが分かります。
あくまでも推計値ではありますが、東京以外の道府県でも表示を確認したり、変化量を分析してみるのも興味深いかもしれません。
参考：