PostGISで気象庁の細分区域&市町村等に対応する重心を求める

震源地と各地の震度情報をOpenLayers3で表示するために、気象庁が独自に分けた区域の重心をPostgreSQLとPostGISで求めようというお話。
今回は、国土数値情報(行政区域データ)をPostgreSQLに入れてQGISで閲覧した前回の続きです。

1. 気象庁の資料からテーブル作成

1. 気象庁防災情報XMLフォーマット 技術資料から個別コード表(jmaxml_20150430_Code.zip)をダウンロード
2. 地震火山関連コード表.xlsのシート24の内容をコピー、新しくファイルを作る
3. AreaInformationCityのコードを隣の列にコピー、下2桁を削って6桁にし、5桁(先頭0埋め)にする
4. CSVで書き出し

テーブルの作成

CREATE TABLE jma_areacode
(
  areaforecastlocale_code integer,
  areaforecastlocale_name character varying(10),
  areaforecastlocale_kana character varying(50),
  areainformationcity_code integer,
  areainformationcity_code_s character varying(5),
  areainformationcity_name character varying(20),
  areainformationcity_kana character varying(50),
  pointseismicintensity_code integer NOT NULL,
  pointseismicintensity_name character varying(30),
  pointseismicintensity_kana character varying(100),
  CONSTRAINT point_pkey PRIMARY KEY (pointseismicintensity_code);

データベースへCSVを読み込む

COPY jma_areacode FROM '/Users/keiichiro/jmaxml_code_s.csv' WITH ENCODING 'sjis' HEADER CSV;

試しに市町村等区域でグループ化してみる

SELECT
  min(AreaForecastLocalE_code) AS area_code, 
  min(AreaForecastLocalE_name) AS area_name, 
  AreaInformationCity_code_s AS city_code, 
  min(AreaInformationCity_name) AS city_name 
FROM
  jma_areacode
GROUP BY
  AreaInformationCity_code_s
ORDER BY
  city_code;

2. 市町村等区域で重心を求める

行政区域(国土数値情報)を気象庁の市町村等区域で結合

CREATE TABLE jma_city_shapes AS
SELECT 
  min(areainformationcity_code) AS city_code,
  n03_007 AS city_code_s,
  min(areainformationcity_name) AS city_name,
  ST_Union(geom) AS geom
FROM
  jma_areacode,
  shapes
WHERE
  jma_areacode.areainformationcity_code_s=shapes.n03_007
GROUP BY
  n03_007
ORDER BY
  n03_007 ASC;
ALTER TABLE jma_city_shapes ADD PRIMARY KEY (city_code);

Core i7-3820QM 2.7GHz, 16GB RAMなMacBook Proで537121ミリ秒、だいたい9分掛かりました。

ここで問題点が2つ。

その1. 五島と甑島の区分
気象庁区分では、
・長崎県佐世保市と長崎県佐世保市宇久島
・鹿児島県薩摩川内市と鹿児島県薩摩川内市甑島
は別域になっているので、それぞれ分割してarea_code,nameとcity_code,nameを修正、追加

730, 長崎県北部, 4220201, 佐世保市
737, 長崎県五島, 4220202, 佐世保市宇久島
770, 鹿児島県薩摩, 4621501, 薩摩川内市
775, 鹿児島県甑島, 4621502, 薩摩川内市甑島

その2. 気象庁と国交省のデータ更新日の差
栃木県下都賀郡岩舟町が栃木市に編入合併された
→栃木市にST_Unionする

気象庁の市町村等区域で重心を求める

ALTER TABLE jma_city_shapes
ADD centroid geometry;
UPDATE jma_city_shapes
  SET centroid = t.c
  FROM (SELECT city_code, ST_Centroid(geom) AS c FROM jma_city_shapes) AS t
WHERE jma_city_shapes.city_code = t.city_code;

3. 細分区域で重心を求める

市町村等を気象庁の細分区域(188区域)で結合

CREATE TABLE jma_area_shapes AS
SELECT 
  areaforecastlocale_code AS area_code,
  max(areaforecastlocale_name) AS area_name,
  ST_Union(geom) AS geom
FROM
  jma_areacode,
  jma_city_shapes
WHERE
  jma_areacode.areainformationcity_code=jma_city_shapes.city_code
GROUP BY
  areaforecastlocale_code
ORDER BY
  areaforecastlocale_code ASC;
ALTER TABLE jma_area_shapes ADD PRIMARY KEY (area_code);

これは819513ミリ秒(だいたい14分)で終了。

気象庁の細分区域(188区域)で重心を求める

ALTER TABLE jma_area_shapes
ADD centroid geometry;
UPDATE jma_area_shapes
  SET centroid = t.c
  FROM (SELECT area_code, ST_Centroid(geom) AS c FROM jma_area_shapes) AS t
WHERE jma_area_shapes.area_code = t.area_code;

テーブルを出力

COPY (SELECT area_code, area_name, ST_Y(centroid) AS lat, ST_X(centroid) AS lon FROM jma_area_shapes ORDER BY area_code ASC) TO '/temp/jma_area_centroid.csv' DELIMITER ',';
COPY (SELECT city_code, city_name, ST_Y(centroid) AS lat, ST_X(centroid) AS lon FROM jma_city_shapes ORDER BY city_code ASC) TO '/temp/jma_city_centroid.csv' DELIMITER ',';

生成したデータはGitHubに置いてます。
https://github.com/9SQ/jma-eqarea-centroid

次回は、OpenLayers3で実際に震源地・震度の表示をします。

NEXT→ OpenLayers 3で気象庁発表の震度をマッピングする

PostgreSQLとPostGISで国土数値情報(行政区域)を扱ってみる

先日(5/30)、小笠原諸島西方沖を震源とする地震がありましたが、そういや気象庁防災情報XMLで震源地や各地の震度情報なども配信していたなぁと思いだし...

震源地と各地の震度情報をOpenLayers3で表示するのを作ろう！

と思い立ったのが4日前。
早速作ろうと色々データを集め始めたのですが、表示の際に一番必要な各地の座標(緯度,経度)は気象庁から公開されていません。
地震が発生したときに気象庁から送信されるXML電文は、震度速報、震源に関する情報、震源・震度に関する情報などがあり、例えば震源・震度に関する情報では...

(前略)
<Pref><Name>東京都</Name><Code>13</Code><MaxInt>5+</MaxInt>
 <Area><Name>小笠原</Name><Code>359</Code><MaxInt>5+</MaxInt>
  <City><Name>小笠原村</Name><Code>1342100</Code><MaxInt>5+</MaxInt>
   <IntensityStation><Name>小笠原村母島</Name><Code>1342103</Code><Int>5+</Int></IntensityStation>
   <IntensityStation><Name>小笠原村父島西町</Name><Code>1342100</Code><Int>4</Int></IntensityStation>
   <IntensityStation><Name>小笠原村父島三日月山</Name><Code>1342101</Code><Int>4</Int></IntensityStation>
  </City>
  </Area>
 <Area><Name>東京都23区</Name><Code>350</Code><MaxInt>4</MaxInt>
  <City><Name>東京千代田区</Name><Code>1310100</Code><MaxInt>4</MaxInt>
(後略)

このようになっています。
ここで、Areaは気象庁が独自に分けた188区域を、Cityは総務省の市区町村コードをベースにした1898区域が指定されます。

参考：気象庁 | 震度情報や緊急地震速報で用いる区域の名称

このXMLには区域コードと区域名だけで、気象庁が公開している技術資料などにも各区域のコードと名称のリストはありますが、座標は無い...

というわけで、無いなら作りましょう。

今回はPostgreSQLと、その拡張でGISオブジェクを格納することができるPostGIS、GISの閲覧、編集、分析機能を持ったQGISを使って、国交省が公開している国土数値情報(行政区域データ)を扱ってみます。

1. PostgreSQL, PostGIS, QGISなどをインストール (Mac/homebrew)

まずはQGISをインストール

brew install pyqt
brew install caskroom/cask/brew-cask
brew cask install qgis

次に、PostgreSQLとPostGIS、pgAdmin3のインストール

brew install postgresql
brew install postgis
brew cask install pgadmin3

2. データベースの作成、データの読み込み

データベースを作ってログインします

postgres -D /usr/local/var/postgres
createdb -E utf8 gis_dataset
psql -U keiichiro -d gis_dataset

エクステンションを設定してPostGISを使えるようにします

CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;

PostGISが読み込めているか確認してみる

SELECT PostGIS_Version();

読めてたらpsqlを抜ける(\q)

国土数値情報 行政区域データ(全国)をダウンロードする
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03.html

N03-140401_GML.zipを展開して、SQLに変換、UTF-8に変換する

unzip N03-140401_GML.zip
cd N03-140401_GML
shp2pgsql -W SHIFT-JIS -c N03-14_140401.shp shapes > N03-14_140401.sql
nkf -W N03-14_140401.sql > N03-14_140401_UTF8.sql

先ほど作ったデータベースに読み込む

psql gis_dataset < N03-14_140401_UTF8.sql

3. QGISで見てみる

QGISを開き、レイヤ->レイヤの追加->PostGISレイヤの追加

[新規]をクリックして、ホスト(localhost)、先ほど作ったデータベース名、ユーザ名を入力して[OK]

先ほど取り込んだデータは、publicスキーマのshapesテーブルに入っているので、これを選択して[追加]

測地系を選択(WGS 84がデフォルトで選ばれているはず)して[OK]

これで、取り込んだ国土数値情報 行政区域データがレイヤとして追加され、画面に表示されるはず。

元となる行政区域データの準備が整ったので、次回は気象庁のデータを取り込んで、重心を求める計算をします。

NEXT→ PostGISで気象庁の細分区域&市町村等に対応する重心を求める

参考
http://qiita.com/nakamods/items/7f25aaaba950c8b8a458
http://morphocode.com/how-to-install-postgis-on-mac-os-x/
http://qiita.com/yellow_73/items/0845451b792f4bc33e90