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 の測地系データをインポート(… >>
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。