2014/07/24
PostGIS >  メートル距離を度単位に変換
経緯度座標で作った地図に、メートル単位で何らかの距離や面積を入れる際の自分の方法。典型的にはスケールバーや 1km 四方など。例として ↓ のデータを R で作業している場合。元データは平成22年国勢調査のさいたま市桜区の境界データで、詳細は 2013/11/30（PL/Python で shp2pgsql を使う）にある。また今日の実行環境は 2014/04/17 を参照。


本題ではないが、↓ が R の maptools パッケージでプロットするまでのコード。色塗りはランダムで意味なし。余白や軸目盛りの微調整が結構手間。
library(maptools)
shp = readShapeSpatial('R:/h22ka11106.shp'
, proj4string=CRS("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"))

windows(width=436/96, height=355/96)
par(cex=0.8, mgp=c(30, 5, 0)/10, plt=c(14, 93, 9, 95)/100)
plot(axes=T, shp, col=rgb(0, 0, 1, runif(nrow(shp))), xaxt='n', yaxt='n')
xat = seq(139.58, 139.64, by=0.02)
yat = seq(35.84, 35.88, by=0.02)
axis(1, at=xat, col=NA, col.ticks=1, labels=sprintf('%.2f°E', xat), tck=0.02)
axis(2, at=yat, col=NA, col.ticks=1, labels=sprintf('%.2f°N', yat), las=1, tck=0.02)

地図の左下が空いているので、（139.58, 35.84）を起点に 1000 m の距離を図示したいとする。PostGIS がインストールされたデータベースからクエリツールを起動して ↓ のクエリを実行すると、概ね 1000 m 相当の緯度・経度距離が出る。
WITH a AS (
	SELECT 139.58 x		-- 経度
		, 35.84 y	-- 緯度
		, 4612 srid	-- 測地投影系
		, 1000 met	-- メートル距離
), b AS (
	SELECT *, ST_Transform(
			ST_SetSrid(ST_Point(x, y), srid)
		  , 4326) :: geography geog
	FROM a
), c AS (
	SELECT *, ST_Transform(
			ST_Envelope(ST_Buffer(geog, met) :: geometry)
		, srid) env
	FROM b
)
SELECT abs(ST_XMax(env) - x) deg_x
	, abs(ST_YMax(env) - y) deg_y
FROM c ;


経度 deg_x が約 0.011°、緯度 deg_y が約 0.009°。この値を R のコードで適宜使う。↓ がクエリの流れ。基本はメートルを単位とする geography 型を使うこと。


↓ 日本語マニュアルでの geography 部分。現行バージョン 2.1 より前だが基本的に同じ。

■ PostGIS 2.0.0マニュアル日本語訳 : ジオグラフィ型
http://www.finds.jp/docs/pgisman/2.0.0/using_postgis_dbmanagement.html#PostGIS_Geography

↓ 先の算出結果を使い、起点（139.58, 35.84）から 1000 m 四方を範囲を追加してみた。文字位置の微調整が結構手間。
xy = c(139.58, 35.84)
xy = c(xy, xy + c(0.011, 0.009)) # <- PostGIS で算出した 1000m 相当の経緯度差
rect(xy[1], xy[2], xy[3], xy[4], lty=2)
xx = xy[c(1,3)]
yy = xy[c(2,4)]
text(adj=c(0.5, 1), x=mean(xx), y=yy[1]-diff(yy)*0.1, label='1000m')
text(adj=c(1, 0.5), x=xx[1]-diff(xx)*0.1, y=mean(yy), label='1000m')


↓ 地図の測地投影系が OpenStreetMap などと同じ世界測地系（SRID=4326）なら、直接 geography 型に変換できるのでクエリが少し短くなる。
-- 測地投影系が 4326 なら、少し短くなる
WITH a AS (
	SELECT 139.58 x		-- 経度
		, 35.84 y	-- 緯度
		, 1000 met	-- メートル距離
), b AS (
	SELECT *, ST_Envelope(ST_Buffer(
			ST_SetSrid(ST_Point(x, y), 4326) :: geography
		, met) :: geometry ) env
	FROM a
)
SELECT abs(ST_XMax(env) - x) deg_x
	, abs(ST_YMax(env) - y) deg_y
FROM b ;

同じことを R 単独で行う方法もあるかもしれないが、特に調べていない。
<< EPSG の測地系データをインポート（…
pgAdmin のクエリツールを電卓代わ… >>