国土数値情報から行政区界を作成する #PostGIS

この記事とだいたい一緒です。

国土交通省のデータを使って、緯度経度から市区町村までを取り出す

ここでは、住所コードからその市区町村の形状を描いてみます。流れとしては以下のようになります。

国土数値情報ダウンロードサービスから都道府県ごとのZIPファイルをダウンロード
ZIPを展開し、Shape形式のデータをPostGISに登録
PostGISを使って形状データを市区町村ごとに結合
GeoJSONで出力して描画

まずは平成25年度のデータをダウンロードします。平成26年度のデータが公開されたらそちらを使った方が良いでしょう。クリックを繰り返すだけの単調作業なので、時間があるときに実施しておきます。47都道府県のZIPファイルを、どこか適当なディレクトリにまとめて保管します。

EC2インスタンスを起動

AWS の EC2 で作業します。インスタンスのイメージIDは Amazon Linux AMI のページに記載されています。その時の一番新しいタイプを使うと便利なことは多いですが、ライブラリのバージョン依存に悩む場合は、ちょっと古いタイプを試してみましょう。--user-data を使って自動アップデートを止める方法は下記の記事が参考になります。

Amazon Linux AMI 2014.03 Released! 早速ハマったこと (Perl・glibc編)

PostGIS を yum でインストールする場合、EPEL の gdal-libs に依存します。
しかし、2014.03 のレポジトリだと EPEL の gdal-libs が依存する poppler のバージョンが合いません。 gdal-libs は poppler-0.12 に依存しますが、最新のインスタンスには poppler-0.22.5 がインストールされています。ひとつ古い 2013.09 は poppler-0.12.4 がインストールされていますので、ソースコードからビルドするのが煩雑な場合はこちらを使うのが良いでしょう。

ということで、以下の構成で awscli のコマンドラインツールを使って EC2 を起動します。パッケージレポジトリのバージョンを固定しておきたいので、--user-data で起動時に自動アップグレードしないよう設定します。

種別	内容
AMI	2013.09 の 64bit
リージョン	アジアパシフィック東京
イメージID	ami-3561fe34
インスタンスタイプ	t1.micro
SSH鍵	aws-tokyo-sandbox
セキュリティグループ	SSHEnabledGroup

EC2インスタンスの起動

aws ec2 run-instances \
    --image-id ami-3561fe34 \
    --instance-type t1.micro \
    --count 1 \
    --key-name ec2-tokyo-sandbox \
    --security-groups SSHEnabledGroup \
    --user-data `cat << EOF | base64
#cloud-config
repo_upgrade: none
EOF`

SSH でログインしたらレポジトリの設定を変更します。

レポジトリ設定

curl -O http://yum.postgresql.org/9.3/redhat/rhel-6-x86_64/pgdg-redhat93-9.3-1.noarch.rpm
sudo rpm -ivh pgdg-redhat93-9.3-1.noarch.rpm
sudo vim /etc/yum.conf  # releasever=latest をコメントアウト
sudo vim /etc/yum.repos.d/pgdg-93-redhat.repo  # $releasever ⇒ latest

次に、PostGIS をインストールします。インストールが成功すると shp2pgsql コマンドが使えるようになっているはずです。

PostGISのインストール

sudo yum install postgresql93-server postgresql93-contrib -y
sudo yum install --enablerepo epel postgis2_93 postgis2_93_devel postgis2_93_utils -y

PostGIS がインストールできたら PostgreSQL を起動します。
作業を簡単に進めるために、アクセス権限を変更しておきましょう。SSHでログインしている場合はユーザー名を指定するだけにします。

データベース設定

sudo service postgresql-9.3 initdb
sudo service postgresql-9.3 start
sudo su - postgres
vi /var/lib/pgsql/9.3/data/pg_hba.conf    # 下記参照

/var/lib/pgsql/9.3/data/pg_hba.conf

local        all        all        trust

設定ファイルを変更したら PostgreSQL を再起動します。

sudo service postgresql-9.3 restart

データベースを作成し、PostGIS 拡張を有効にします。

データベース作成

createdb -U postgres geo
cat << EOF | psql -U postgres -d geo
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
SELECT postgis_full_version();
EOF

これでデータベースの設定ができました。

ZIPファイルの中身をPostGISに登録

国土数値情報ダウンロードサービスからダウンロードしたZIPファイルをEC2のインスタンスに転送します。転送できたら、ZIPファイルを展開して shp2pgsql でSQLファイルを作成し、PostGISに登録します。

for f in N03-*.zip
do
    unzip -d ${f%.zip} $f
done

for f in N03-*/*.shp
do
    t=`basename $f`
    sql=${t%.shp}.sql
    shp2pgsql -W "Shift_JIS" $f | psql -U postgres -d geo
done

ファイル名に応じたテーブルが生成されますので、47個のテーブルが生成されます。

データを見てみると、いくつかの都道府県のデータには注意が必要であることが分かります。 (北海道地方の支庁・地方事務所一覧 - Wikipedia)

北海道の支庁と振興局が混在している

SELECT distinct n03_002
FROM "n03-13_01_130401"
WHERE n03_002 IS NOT NULL
ORDER BY n03_002;

    n03_002
----------------
 上川支庁
 十勝支庁
 宗谷支庁
 後志支庁
 日高支庁
 根室振興局
 根室支庁
 檜山支庁
 渡島支庁
 留萌支庁
 石狩支庁
 空知支庁
 網走支庁
 胆振支庁
 釧路支庁
 釧路総合振興局
(16 rows)

香川県の郡・政令市名が香川県のレコードがある

\x

SELECT n03_001, n03_002, n03_003, n03_004, n03_007
FROM "n03-13_37_130401"
WHERE n03_001 = n03_003;

-[ RECORD 1 ]---
n03_001 | 香川県
n03_002 |
n03_003 | 香川県
n03_004 | 丸亀市
n03_007 | 37202

所属未定地の確認

SELECT n03_001, n03_002, n03_003, n03_004, n03_007 FROM "n03-13_13_130401" WHERE n03_007 IS NULL;

-- 結果は省略

たとえば、住所コードで集約したときに名称が不一致になるレコードを抽出すると、以下のレコードをリストアップできます。名寄せが必要なレコードには UPDATE 文を発行しておきましょう。

n03_001	n03_002	n03_003	n03_004	n03_007	レコード数
北海道	根室振興局		根室市	01223	76
北海道			根室市	01223	464
北海道	釧路支庁	厚岸郡	浜中町	01663	139
北海道	釧路総合振興局	厚岸郡	浜中町	01663	1
三重県		北牟婁郡	紀北町	24543	416
三重県		北牟婁郡	紀北町入会地	24543	110
香川県		香川県	丸亀市	37202	1
香川県			丸亀市	37202	14

市区町村には同一名称もありますので、集約させるキーには注意しましょう。また、テーブル定義で UNIQUE 制約を付ける場合は慎重に検討しましょう。

同一名称の市区町村一覧 - Wikipedia

PostGISを使って形状データを市区町村ごとに結合

都道府県ごとのテーブルにあるレコードを住所コードで集約します。PostGIS では複数の形状データを ST_Union で結合できます。
今回の場合は、住所コード (n03_007 フィールド) で GROUP BY し、 geom フィールドを結合させます。変換先のテーブルは以下のように定義します。

変換先のテーブル定義

CREATE TABLE address (
    code varchar(5) PRIMARY KEY CHECK (length(code) IN (2, 5)),
    geom geometry(MULTIPOLYGON, 4326) NOT NULL
);

変換クエリはSQLファイルに保存し、 psql コマンドの "-v" オプションを使って変数をループさせます。都道府県ごとに別々のテーブルになっているためです。"table" にはテーブル名、"code" には都道府県コードを指定します。

transform.sql

INSERT INTO address
SELECT n03_007, ST_Union(ST_SetSRID(geom, 4326))
FROM :"table"
WHERE n03_007 IS NOT NULL
GROUP BY n03_007;

INSERT INTO address
SELECT :'code', ST_Union(ST_SetSRID(geom, 4326))
FROM :"table"
WHERE n03_007 IS NULL
GROUP BY n03_001;

都道府県ごとにクエリを実行

for i in {01..47}
do
    psql -U postgres -d geo -f transform.sql -v table="n03-13_${i}_130401" -v code=$i
done

GeoJSONで出力して描画

PostGIS の ST_AsGeoJSON を使うと形状データを GeoJSON で出力できます。
市区町村ごとにファイルを作成してみましょう。

echo "SELECT code, ST_AsGeoJSON(geom) FROM address" | psql -U postgres -d geo -t -A --field-separator=$'\t' > address.txt
awk '{ print $2 >> $1".geojson" '} address.txt

出来上がった GeoJSON ファイルを geojson.io に貼付けてみると、市区町村の形状を確認できます。自前で描画したい場合は Leaflet や D3 を使って JavaScript を書けば良いでしょう。

住所コードに対応する名称は総務省の全国地方公共団体コードで定義されています。利用する住所データに対応する年度のコード表を使うとマッピングできます。なお、いくつか注意点があります。

年度が異なる場合は市区町村の合併などに適応できない場合があります。
総務省のデータには６桁目のチェックディジットがあります。
総務省のデータには「郡」がありません。
総務省のデータは都道府県市区町村名はカタカナですが、政令指定都市名はひらがなで併記されています。

PostgreSQL などのデータベースでデータを管理しておくと、JOIN を使った結合処理を簡単に記述できる点は便利だと言えます。表示するラベルを工夫したい場合はデータを出力するクエリを頑張って記述しましょう。