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 のクエリツールを電卓代わ… >>
×

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