# 地理空間データの種類

# (1)ベクタモデル・・・地形を構成する「ポイント・ライン・ポリゴン」のデータからなる

# ポイントデータ:(x,y)座標のデータフレーム,RではSpatialPointsDataFrameというクラス

# ラインデータ:連結する順序に並んだポイントデータ,RではSpatialLinesDataFrameというクラス

# ポリゴンデータ:始点と終点が同一座標のラインデータ,RではSpatialPolygonsDataFrameというクラス

# (2)ラスタモデル・・・リモートセンシングのような,グリッドを用いて連続面を表現するデータ

# Rでは多くのパッケージでspクラスが利用されている

#

# シェープファイル

# ESRI社ArcViewのデジタル地図データ形式が広く普及している

# ファイル形式は公開されており,属性データにdBASE形式を採用している

# 幾何データのメインファイル(.shp),インデックスファイル(.shx),属性データ(.dbf)からなる

# Rでは,maptools,rgdal,shapefilesなどのパッケージを使って読み込む

# maptools

# デモファイルを多く含んでおり,system.file()でフルパスを得ることができる

# getinfo.shape()でシェープファイルの情報を得る

# シェープファイルの読込は種類に応じてreadChapePoints(),readShapeLines(),readShapePoly()を使う

# readShapeSpatial()は種類を問わず利用可能

# sidsデータの例(シェープファイルの読込)

library(maptools)

f <- system.file("shapes/sids.shp",package="maptools") # デモファイルのフルパスを得る

getinfo.shape(f) # シェープファイルの情報を得る

map <- readShapeSpatial(f) # シェープファイルを読み込む

plot(map) # 地図を作図

# 準拠楕円体と測地系

# 位置として使われる(緯度,経度)は準拠楕円体と測地系で決まる

# 準拠楕円体と測地系が異なれば,同じ値の(緯度,経度)でも別の場所を示す

# 日本ではこれまで「準拠楕円体=ベッセル楕円体,測地系=日本測地系(Tokyo)」が使われてきた

# 2002/4/1から「準拠楕円体=GRS80,測地系=世界測地系(日本測地系2000)」に移行

# 日本の地図データにはこれら2種類の測地系が混在するため注意が必要(長崎辺りで420mずれる)

# 地図投影法

# 地球は3次元立体であるため,2次元平面定の地図に表現するには投影が必要

# 投影により面積・角度・距離・形に歪みが生じる

# 目的(どの情報の歪みを少なくするか)に応じて,様々な投影法が提案されている

# 座標系

# 空間内の点の位置を数値の組で表すために定められるパラメータセット

# Rのspクラスの座標系はCRSクラスで定義されている

# UTM:横メルカルトル図法に基づく座標系

# 平面直角座標系:狭い範囲なら平面と扱った測量が用いられる(日本では,1万分の1より大きな地図で採用)

# 座標系の変換

# 投影法・座標系を相互に変換することを幾何変換という

# sidsデータの例(地理座標系の地図を作成)

library(maptools)

f <- system.file("shapes/sids.shp",package="maptools")

pj <- CRS("+proj=longlat +ellps=clrk66") # 投影法=地理座標系,準拠楕円体=Clarke1866

map <- readShapeSpatial(f,proj4string=pj) # 座標系情報をもつSpatialPolygonsDataFrameの作成

plot(map,axes=T)

title(main="米国ノースカロライナ州(地球座標系)",xlab="Longitude",ylab="Latitude")

# sidsデータの例(rgdalのspTransform()を用いて地理座標系からUTM第17帯に幾何変換)

library(rgdal)

pj2 <- CRS("+proj=utm +zone=17 +ellps=WGS84 +units=km") # 投影法=UTM第17帯,準拠楕円体=WGS84,単位距離=km

map2 <- spTransform(map,pj2)

plot(map2,axes=T)

title(xlab="Easting(km)",ylab="Northing(km)",main="米国ノースカロライナ州(UTM zone 17)")

# sidsデータの例(投影法をランベルト正角円錐図法に変更)

pj3 <- CRS("+proj=lcc +lat_1=40 +lon_0=-95 +datum=NAD27 +units=km")

map3 <- spTransform(map2,pj3) # 投影法=LCC,測地系=NAD27,地図単位=km,原点=(緯度32,経度-95)

plot(map3,asex=T)