2010年11月15日月曜日

Rで座標系の変換

rgdalパッケージで読み込んだオブジェクトの座標系の情報はproj4stringというスロットに格納されている。
例えば以下のように"work.data"という名前の変数を作ってオブジェクトを格納した場合は


library(rgdal)

work.data <- readOGR("glacier_polygon.shp", layer="glacier_polygon")
work.data@proj4string 
 
 
UTM座標zone 45 northの場合は以下のようにターミナルに表示される +proj=utm +zone=45 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0   
これをWGS84地理座標系に変換したい場合は以下のように入力すると 座標系の変換された"work.data.wgs84"というオブジェクトが得られる。 new.crs <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") work.data.wgs84 <- spTransform(work.data, CRS=new.crs) 
 
  
XY座標を確認してみると以下のように表示され、 座標系がDecimal degree(十進経緯度)になっているのが確認できる。

> coordinates(work.data2)
          [,1]     [,2]
 [1,] 86.88998 27.93553
 [2,] 86.93320 27.87881
 [3,] 86.87982 27.87656
 [4,] 86.70286 28.02150
 [5,] 86.63027 28.03838
 [6,] 86.57888 28.05175
 [7,] 86.59177 27.77458
 [8,] 86.85842 27.98194

ちなみにspTransform関数で使っているCRS引数の中の文字列は、 projライブラリで使われる座標系を指定する文字列で、QGISのオプションのCRSの項目を見ると、 それぞれの座標系ごとの文字列の書き方が容易に確認できて便利。

0 件のコメント:

コメントを投稿