2014/07/27
PostGIS >  EPSG の測地系データを R で可視化
一昨日、昨日の続き。EPSG の全世界の測地系データを、シェープファイルのポリゴンを含めてインポートしたので、確認で日本の平面直角座標の全ゾーンを R でプロットしてみる。実行環境は 2014/04/17、ただし PostGIS は最新の 2.1.3。

↓ まず epsg_coordinatereferencesystem という長いテーブル名の列 coord_ref_sys_name で検索。JGD2000 / Japan Plane Rectangular で始まり PostGIS の SRID が存在する測地系を抽出すると I〜XVIII と南鳥島の計19ゾーンあった。古い日本測地系の場合は Tokyo / Japan Plane Rectangular で検索する。
SELECT t3.srid, t1."AREA_NAME", t1."VERSION"
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 t2.coord_ref_sys_name LIKE 'JGD2000 / Japan Plane Rectangular%'
ORDER BY 1, 2 ;


↓ 次に R から呼ぶクエリを作成。rgeos パッケージの readWKT 関数で読み込むため、PostGIS の側でポリゴンを ST_AsText する（この出力が WKT 形式になる）。同時に全ポリゴンを囲む矩形も渡す。これは一行だけでいいが、クエリを一回で済ませるためウィンドウ関数を使い、全行で集約関数 ST_Extent を実行。実際使うのは先頭行の bbx 列だけ。
SELECT t3.srid
	, substring(t2.coord_ref_sys_name, '[IVX]+$') znn
	, ST_AsText(ST_Extent(geom) over()) bbx
	, ST_AsText(t1.geom) wkt
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
WHERE t2.coord_ref_sys_name LIKE 'JGD2000 / Japan Plane Rectangular%'
ORDER BY 1, 2 ;


続いて、R での作業。三つのパッケージ maptools, rgeos, RPostgreSQL を使う。バージョンは下記を参照。出力を三パッケージに絞るのが面倒だったので、そのまま。
subset(installed.packages(), select='Version')


↓ 作成したスクリプト。先のクエリを実行して結果を取得し（データフレームになる）、まず一行目の bbx を描画して全体の領域を作る。次にポリゴンを一行ずつ readWKT で読み込みプロット。例によって処理本体よりレイアウトの微調整が多い。
library(maptools)
library(RPostgreSQL)
library(rgeos)

con = dbConnect(PostgreSQL(), dbname='xxx', port=xxx, user='xxx', password='xxx')
sql = "
	SELECT t3.srid
		, substring(t2.coord_ref_sys_name, '[IVX]+$') znn
		, ST_AsText(ST_Extent(geom) over()) bbx
		, ST_AsText(t1.geom) wkt
	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
	WHERE t2.coord_ref_sys_name LIKE 'JGD2000 / Japan Plane Rectangular%'
	ORDER BY 1, 2 ;
"
dat = dbGetQuery(con, sql)
invisible(dbDisconnect(con))

windows(width=436/96, height=355/96)
par(cex=0.8, las=1, mgp=c(20, 4, 0)/10, plt=c(10, 98, 7, 98)/100, tck=0.01)
bbx = readWKT(dat$bbx[1])
plot(axes=T, border=NA, x=bbx, xaxt='n', yaxt='n')
bbr = round(sapply(bbox(bbx), unlist), -1)
xat = seq(bbr[1], bbr[3], by=10)
yat = seq(bbr[2], bbr[4], by=10)
axis(1, col=NA, col.ticks=1, at=xat, labels=sprintf('%s°E', xat))
axis(2, col=NA, col.ticks=1, at=yat, labels=sprintf('%s°N', yat))

len = nrow(dat)
col = adjustcolor(rainbow(len), alpha=0.5)
for(i in 1:len){
	ply = readWKT(dat[i,'wkt'])
	cen = gCentroid(ply)
	plot(add=T, col=col[i], x=ply)
	text(cen$x, cen$y, label=dat[i,'znn'])
}

# 凡例を両側に分けて描画
brk = 7
loc = c('topleft', 'topright')
lis = list(1:brk, (brk+1):len)
for(i in 1:2){
	j = lis[[i]]
	legend(bty='n', x=loc[i], fill=col[j]
		, legend=paste(dat$srid, dat$znn)[j]
		, title='SRID & Zone')
}

↓ 実行結果。各ゾーンのポリゴンが予想以上に細かく作られていた。二枚目は大きいサイズの出力。画像だけ開けば原寸大になる。



↓ グラフィックウィンドウの結果を EPS ファイルで保存し Illustrator で読み込み、アウトライン表示したところ。各ポリゴンの細かい境界線がきちんと接合するのが分かった。また海岸線に少しバッファを足した領域になっているので、東京湾はなくなり佐渡島はずんぐりした形。


というわけで少なくとも日本の平面直角座標系なら、昨日の最後のように経緯度を指定すれば一つのポリゴンが返る、つまりゾーンを一意に絞れることが分かった。下記のような対照表を調べる作業は不要になりそう。

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

ところで Spatial Reference で日本の平面直角座標系のゾーン IX を検索すると、南関東から下の範囲が赤い線で示される。一方、今日の出力結果では IX はもっと北まで含む。上記国土交通省告示で IX が北関東や福島県を含むので、Spatial Reference の方が不正確。

http://spatialreference.org/ref/epsg/2451/

<< 実行環境まとめ（PostGIS 更新）
EPSG の測地系データをインポート（… >>