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で気象庁発表の震度をマッピングする