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