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