分析用データ(作成中)

このページの位置づけ

空間分析用パッケージのインストール

install.packages("ctv")
library(ctv)
install.views("Spatial")

ライブラリ呼び出し

library(maptools)
library(splancs)
library(rgdal)
library(gpclib)

全国市区町村境界データ

  • ESRIジャパンのHPより全国市区町村境界データをダウンロードし解凍する
  • Windowsでは、C:/Program Files/R/R-2.7.2/japan_ver62/ にファイルを置く
  • Macでは、/Library/Frameworks/R.framework/japan_ver62/ にファイルを置く

shapeファイルの読み込みとプロット

# Windows
jpn <- readShapeSpatial("C:/Program Files/R/R-2.7.2/japan_ver62/japan_ver62.shp",
IDvar="JCODE",proj4string=CRS("+proj=tmerc +ellps=GRS80 +units=m +datum=JGD2000"))
# jpn <- readShapeSpatial("jpn_pref.shp",IDvar="PREF_CODE",
# proj4string=CRS("+proj=tmerc +ellps=GRS80 +units=m +datum=JGD2000"))
# jpn <- readShapeSpatial("jpn_pref.shp",IDvar="PREF_CODE",
# + proj4string=CRS("+proj=longlat +ellps=GRS80 +units=m +datum=JGD_2000"))
lps <- coordinates(jpn)
plot(jpn, xlim=c(min(lps[,1]),max(lps[,1])), 
ylim=c(min(lps[,2]),max(lps[,2])))
# jpn <- readShapeSpatial("C:/Program Files/R/R-2.7.2/japan_ver62/japan_ver62.shp",
# IDvar="JCODE",proj4string=CRS("+proj=longlat +ellps=GRS80 +units=m +datum=WGS84"))
# jpn <- readShapeSpatial("C:/Program Files/R/R-2.7.2/japan_ver62/japan_ver62.shp",
# IDvar="JCODE",proj4string=CRS("+proj=eqdc +ellps=GRS80 +units=m"))
# jpn <- read.shape("C:/Program Files/R/R-2.7.2/japan_ver62/japan_ver62.shp")
# jpn <- readShapeSpatial("C:/Program Files/R/R-2.7.2/japan_ver62/japan_ver62.shp")
# Mac
# jpn <- read.shape("/Library/Frameworks/R.framework/japan_ver62/japan_ver62.shp")
# jpn <- readShapeSpatial("/Library/Frameworks/R.framework/japan_ver62/japan_ver62.shp",
# IDvar="JCODE",proj4string=CRS("+proj=loglat +ellps=GRS80 +units=m"))
# プロット
# plot(jpn)

空間データの処理1

# 座標系などを定義
jpn.prj = spTransform(jpn,CRS("+proj=longlat +ellps=GRS80 +units=m +datum=WGS84"))
jpn.prj.ps = SpatialPolygons2PolySet(jpn.prj)
plotPolys(jpn.prj.ps, proj =  TRUE, col="wheat1",
xlab="longitude",ylab="latitude",
xlim=c(min(lps[,1]),max(lps[,1])),  ylim=c(min(lps[,2]),max(lps[,2])))
attr(jpn.prj.ps, "projection") <- "LL"
# attr(pref.ps, "zone") <- 53
# 面積を計算(UTM zone 53でkm^2)
jpn.prj.area = calcArea(jpn.prj.ps,rollup=1)
jpn.prj.area
# 市町村境界データを都道府県境界データにディゾルブする
# 都道府県名でディゾルブする方法
pref <- unionSpatialPolygons(jpn, jpn$PREF)
sapply(slot(pref, "polygons"), function(i) slot(i, "ID"))
# 都道府県コードでディゾルブする方法
# jpn1 <- jpn
# row.names(jpn1) <- trunc(jpn1$JCODE / 1000)
# pref1 <- unionSpatialPolygons(jpn1, round(jpn1$JCODE / 1000))
# 属性テーブルの結合
# 社会経済データの読み込み
pref_soc <- read.table("pref_soc.csv", sep=",",header=T)
str(pref_soc)
# 属性テーブルのマッチング
id <- match(pref_soc$PREF1, sapply(slot(pref, "polygons"),
function(i) slot(i, "ID")))
pref_soc1 <- pref_soc[id,]
row.names(pref_soc1) <- pref_soc$PREF1
# row.names(pref_soc1) <- pref_soc1$CODE1
# ディゾルブしたマップに属性データを定義
pref_soc1.df <- as(pref_soc1, "data.frame") 
pref.spdf <- SpatialPolygonsDataFrame(as(pref, "SpatialPolygons"), data=pref_soc1.df) 
# 都道府県庁緯度経度ポイントデータ
pref.att <- read.table("pref_gov.csv",header=T,sep=",")
pref.mat <- cbind(pref.att$COORD.X, pref.att$COORD.Y)
row.names(pref.mat) <- 1:nrow(pref.mat)
str(pref.mat)
pref.sp <- SpatialPoints(pref.mat,
proj4string = CRS("+proj=longlat +ellps=GRS80"))
points(pref.sp)

空間データの処理2

# 横浜市だけを切り取る
yoko <- jpn[jpn$JCODE > 14100,]
yoko <- yoko[yoko$JCODE < 14120,]
writeSpatialShape(yoko, "C://Program Files/R/R-2.7.2/yoko")
# 属性テーブルを結合する
# 横浜市社会経済属性データを読み込む
yoko_soc <- read.table("yoko_soc.csv", sep=",", header=T)
row.names(yoko_soc) <- yoko$JCODE
# IDをマッチングさせる
ID.match <- match(yoko$JCODE, yoko_soc$JCODE)
yoko_soc1 <- yoko_soc[ID.match, ]
row.names(yoko_soc1) <- yoko$JCODE
# 属性テーブルを結合させる
yoko1 <- spCbind (yoko, yoko_soc1)
names(yoko1)
summary(yoko1)
# セントロイド(重心)を作成する
# 既存シェープファイルの読み込み
yoko.shp <- read.shape("yoko.shp")
yoko.ply <- Map2poly(yoko.shp)
plot(yoko.ply)
# セントロイドの作成
yoko.pc <- get.Pcent(yoko.shp)
CentCoords <- matrix(c(yoko.pc[,1],yoko.pc[,2]),
nrow=length(yoko.pc[,1]),ncol=2,byrow=FALSE) 
ProjectedCentroids <- spTransform(SpatialPoints(coordinates(CentCoords)))
yoko.cnt <- SpatialPoints(coordinates(CentCoords))
plot(yoko.cnt,add=T)

都道府県

# plot(pref, xlim=c(min(lps[,1]),max(lps[,1])), 
# ylim=c(min(lps[,2]),max(lps[,2])))
pref.ps <- SpatialPolygons2PolySet(pref)
dev.set(1) 
plotPolys(pref.ps, proj = TRUE, 
col="salmon1",xlab="longitude",ylab="latitude")
# 面積を計算
attr(pref.ps, "projection") <- "LL"
# attr(pref.ps, "zone") <- 53
pref.area = calcArea(pref.ps ,rollup=1)
# writeSpatialShape(pref.ps, "C://Program Files/R/R-2.7.2/pref.ps")
pref1.fid <- spChFIDs(pref.spdf, as.character(pref.spdf[[1]]))
pref2 <- spCbind(pref1.fid, id)
pref2 <- merge(pref.spdf, pref_soc1.df)
pref2 <- spCbind(pref1.fid, pref_soc1.df)

神奈川県

# library(maptools)
# 属性から切り取る場合
jpn.ply<-Map2poly(jpn)
jpn.att<-jpn$att.data
jpn14.ply<-subset(jpn.ply, jpn.att$JCODE>14100)
jpn14.att<-subset(jpn.att, jpn.att$JCODE>14100)
pref14.ply<-subset(jpn14.ply,jpn14.att$JCODE<14999)
pref14.att<-subset(jpn14.att, jpn14.att$JCODE<14999)
plot(pref14.ply)
# plot(pref14.ply, proj="+proj=eqdc +ellps=GRS80 +units=m")
writeSpatialShape(yoko.ply,yoko.att,"C:/Program Files/R/R-2.7.2/yoko0")
# writePolyShape(yoko.ply,yoko.att,"C:/Program Files/R/R-2.7.2/yoko0")
kana.shp <- read.shape("kana.shp")
# kana.shp <- readShapeSpatial("kanagawa.shp", IDvar="JCODE",
# proj4string=CRS("+proj=longlat +ellps=GRS80"))
kana.ply <- Map2poly(kana.shp)
# library(maptools)
kana.Pcent <- get.Pcent(kana.shp)
CentCoords <- matrix(c(kana.Pcent[,1],kana.Pcent[,2]),
nrow=length(kana.Pcent[,1]),ncol=2,byrow=FALSE) 
ProjectedCentroids <- spTransform(SpatialPoints(coordinates(CentCoords))
kana.cnt <- SpatialPoints(coordinates(CentCoords))
plot(kana.shp)
plot(kana.cnt,add=T)
# library(gpclib)
# np<-length(kana.ply)
# kana.union <- as(kana.ply[[1]],"gpc.poly")
# for(i in 2:np){
#   poly <- as(kana.ply[[i]],"gpc.poly")
#   kana.union<-union(kana.union, poly)
# }
# 面積計算

横浜市

# library(maptools)
jpn.ply<-Map2poly(jpn)
jpn.att<-jpn$att.data
jpn14.ply<-subset(jpn.ply, jpn.att$JCODE>14100)
jpn14.att<-subset(jpn.att, jpn.att$JCODE>14100)
yoko.ply<-subset(jpn14.ply,jpn14.att$JCODE<14120)
yoko.att<-subset(jpn14.att, jpn14.att$JCODE<14120)
plot(yoko.ply)
writeSpatialShape(yoko.ply,yoko.att,"C:/Program Files/R/R-2.7.2/yoko")
# writePolyShape(yoko.ply,yoko.att,"C:/Program Files/R/R-2.7.2/yoko")

都道府県データ

# 都道府県庁緯度経度ポイントデータ
pref.att <- read.table("pref_gov.csv",header=T,sep=",")
pref.mat <- cbind(pref.att$COORD.X, pref.att$COORD.Y)
row.names(pref.mat) <- 1:nrow(pref.mat)
str(pref.mat)
pref.sp <- SpatialPoints(pref.mat,
proj4string = CRS("+proj=longlat +ellps=GRS80"))
# 既存シェープファイルの読み込み
# 都道府県境界データ
jpn_pref.shp <- read.shape("jpn_pref.shp")
jpn_pref.ply <- Map2poly(jpn_pref.shp)
plot(jpn_pref.ply)
# 都道府県庁緯度経度ポイントデータ
pref_gov.shp <- read.shape("pref_gov.shp")
pref_gov.pnt <- Map2points(pref_gov.shp)
points(pref_gov.pnt)

参考文献・サイト

商標・免責

  • RはR Foundationの商標です。R Development Core Team (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org を参照してください。Rを開発されたコアチームの方々、パッケージの開発をされた方々、普及に尽力されている方々に感謝の意を表します。
  • その他、本書で言及している製品名、商標、および登録商標は、権利所有者が権利を有します。各ソフトウェアのライセンスサイトなどを参照してください。
  • 本書で記載しているソフトウェアなどの利用に関して、万一障害が生じても、出版社および筆者は一切責任を負いません。
  • 日本の市町村境界データの利用に関しては、ESRIジャパン株式会社の全国市区町村界データを使用しました。データの著作権はESRIジャパン株式会社に帰属します。

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-05-25 (火) 00:51:49 (3345d)