# CRSクラス

# PROJ.4 : 100以上の座標系をサポートする地図投影変換ライブラリ

# RのspクラスはPRO.J4を利用して地図座標系を操作する

# RではCRSクラスが用意されていて,proj4stringオプションで指定する

# SpatialPoints(),SpatialPlygons()でspクラスを作成

# spTransform()でspクラスの座標系変換を行う

# PROJ.4の表現式はRとは異なり,パラメータは「+パラメータ=値」で表される

# ここでいうパラメータとは「地図投影法,準拠楕円体,測地系,原点,縮尺係数」など

# 値のないパラメータは「+パラメータ」とされる

# http://trac.osgeo.org/proj/

# CRSクラスオブジェクトの作成

# spパラメータのCRS()を使う

# 地理座標系

library(sp)

pj <- CRS("+proj=longlat +datum=WGS84")

class(pj)

# rgdalパッケージのCRSargs()で詳細が確認できる

library(rgdal)

CRSargs(pj)

# UTM

pj <- CRS("+proj=utm +zone=52 +datum=WGS84")

CRSargs(pj)

# 平面直角座標系

# 日本を19の地域で分割する(系番号:I~XIX)

# CRSには系番号は登録されていないが,対応するEPSGが利用可能

# 広島の平面直角座標系はⅢ系(EPSG=30163)←旧日本測地系(Tokyo測地系)

pj <- CRS("+init=epsg:30163")

CRSargs(pj)

# 地図単位をkmに変更

pj.km <- CRS("+init=epsg:30163 +units=km")

CRSargs(pj.km)

# EPSGのリスト

library(rgdal)

epsg <- make_EPSG()

# 日本の旧測地系(Tokyo測地系)と新測地系(日本測地系2000)に対応するEPSG

epsg[grep("Japan Plane Rectangular CS III",epsg$note),]

# 広島市地図の例

pj <- CRS("+init=epsg:30163 +units=m") # 座標系は「旧日本測地系のⅢ系」

map <- readShapeSpatial("river.shp",proj4string=pj) # river.shpの出典は佐藤裕哉先生

pj.km <- CRS("+init=epsg:30163 +units=km +x_0=-26721 +y_0=+178395") # 座標系の定義

map2 <- spTransform(map,pj.km) # 座標系の変換

plot(map2,axes=T)

points(0,0,pch=19,col="red")

xy <- unique(cbind(data$xx,data$yy))/10 # 被爆時所在地の座標(km)

points(xy,pch=20,col="gray80")

title(main="原点が爆心の平面直角座標系(km)")