EPSG の測地系データをインポート（2） : 研究に使うポスグレ

2014/07/26
PostGIS > EPSG の測地系データをインポート（2）

昨日の続き。実行環境は 2014/04/17 を参照、ただし PostGIS は最新の 2.1.3。昨日やり残したシェープファイルのインポートを行い、任意の経緯度から使える測地系を検索できるか試す。↓ の EPSG supporting files、Area polygons in Shapefile format をダウンロードする。

↓ ZIP に圧縮されていて、解凍すると約 50 MB になる。

↓ 2013/11/30（PL/Python で shp2pgsql を使う）と同じ方法でインポートする。事前準備として PL/Python のイントールが必要。解凍したシェープファイルを適当な場所に置き、そのパスを下記の無名コードブロック中に設定する。ZIP にある全てのファイルは必要なく、dbf、shp、shx があればいい。測地系はインポートが成功したら別途行う。

-- PL/Python インストール
CREATE EXTENSION plpython3u ;

-- PL/Python から shp2pgsql を起動してインポート
DO LANGUAGE plpython3u $D$
shp = 'R:/TMP/EPSG_Polygons'	# ファイルパス
enc = 'SJIS'			# 文字コード
tbn = 'public.polygons'		# インポート先テーブル。スキーマ名含む。引用符なし
tmp = 'R:/TMP/epsg.sql'		# shp2pgsql.exe が出力する SQL ファイル

import subprocess
import codecs
import os

cmd = ['cmd', '/c', 'C:/Program Files/PostgreSQL/9.3/bin/shp2pgsql.exe'
, '-i', '-k', '-W', enc, shp, tbn, '>', tmp]
prc = subprocess.Popen(cmd)
prc.wait()
sql = codecs.open(tmp, 'r', 'utf-8').read()
sql = sql.replace('BEGIN;', '')	# shp2pgsql の出力にある BEGIN; COMMIT; を省く
sql = sql.replace('COMMIT;', '')
plpy.execute(sql)
os.remove(tmp) # SQL ファイル残す場合はコメントアウト
$D$ ;

上を実行して、8.5 秒かかったが無事完了した。なぜか文字コードを LATIN1 や UTF-8 にすると失敗した。シェープファイルの中身ではなく OS と関係があるのかもしれない。ともかく ↓ が作成されたテーブル。行数は 3000 余でそれほど多くない。元のファイルサイズが約 50 MB なのでデータ密度が濃いのだろう。

次に geom 列の測地系を設定する。シェープファイルの prj をエディタで開くと ↓ 設定内容から SRID=4326 と分かる。

↓ 普通に SQL で geom 列を UPDATE し、念のためデータ型を限定する。MULTIPOLYGON があるので POLYGON 型ではエラーになる。データ量が多いので約 7.7 秒かかった。

UPDATE polygons SET geom = ST_SetSrid(geom, 4326) ;
ALTER TABLE polygons ALTER geom TYPE geometry(MULTIPOLYGON, 4326) ;

適当な地点を設定して、geom 列が含む（交差する）エリアを検索するテスト。地点は、一昨日の例に使ったさいたま市桜区の地図からとった。その時の測地系が 4612 だったので、一応 4326 に変換した（実際はほとんど同じ）。

WITH a AS (
	-- 対象地点を設定
	SELECT ST_Transform(
		ST_SetSrid(ST_Point(139.6, 35.86), 4612)
	, 4326) poi
)
SELECT * FROM polygons, a
WHERE ST_Intersects(geom, poi)
ORDER BY 1 ;

↑ 24 行もあった。確かに日本付近のエリアがあるが世界全体のもあり、これだけでは経緯度から測地系をきちんと設定できるか分からない。そこで昨日の最後のクエリと同様に、EPSG 本体のデータおよび PostGIS のテーブル spatial_ref_sys と結合して、下のように検索する。

SELECT t3.srid, t1."AREA_NAME", t2.coord_ref_sys_name
FROM polygons t1
JOIN epsg_coordinatereferencesystem t2
	ON t1."AREA_CODE" = t2.area_of_use_code
JOIN spatial_ref_sys t3
	ON t2.coord_ref_sys_code = t3.srid
			-- PostGIS の SRID と結合
WHERE ST_Intersects(t1.geom
    , ST_Transform(ST_SetSrid(ST_Point(139.6, 35.86), 4612), 4326))
ORDER BY 1, 2 ;

↑ 先ほどの 24 エリアすべてについて、対応する SRID および EPSG データがあった。先頭の 2451 は Japan Plane Rectangular で平面直角座標のゾーン IX。念のため下記と対照すると、確かに埼玉県はゾーン IX に当たる。

■ 平面直角座標系（平成十四年国土交通省告示第九号）
http://www.gsi.go.jp/LAW/heimencho.html

平面直角座標系のゾーン区分は行政界が使われているなど複雑で、今日のシェープファイルで完全にカバーできるか不明だが、上の結果からは概ね使えそう。次回、日本付近の測地系エリアを R で可視化してチェックする予定。

<< EPSG の測地系データを R で可視化

EPSG の測地系データをインポート（… >>