第2章プログラム

全体を通じて使用するパッケージ

library(spdep)
# library(sp)

2.1 空間データの基本構造

# jpn_pref <- readShapePoly("jpn_pref.shp")
# pref_pnt <- readShapePoints("pref_pnt.shp")
# pref_gov <- read.table("pref_gov.txt",",",header=TRUE, row.names=2)
# coords <- as.matrix(cbind(pref_gov$X,pref_gov$Y))
# pref.tri.nb <- tri2nb(coords, row.names=rownames(pref_gov))
# par(mfrow=c(1,3))
# plot(pref_pnt)
# plot(pref.tri.nb, coords)
# plot(jpn_pref, col="grey")
# ライブラリ
library(gpclib)
# 空間オブジェクトの作成
# 点データ
px <- runif(30,0,10)
py <- runif(30,0,10)
pz <- as.data.frame(px+py)
colnames(pz) <- c("pz")
pnt_xy <- cbind(px, py)
pnt_sp <- SpatialPoints(data.frame(px, py))
pnt_spdf <- SpatialPointsDataFrame(pnt_sp, pz)
# 図2.1
plot(pnt_spdf, cex=2, pch=19, lwd=3)
# 線データ
line1 <- cbind(runif(4,0,10), runif(4,0,10))
line2 <- cbind(runif(4,0,10), runif(4,0,10))
line3 <- cbind(runif(4,0,10), runif(4,0,10))
line1_ln <- Lines(list(Line(line1)), "line1")
line2_ln <- Lines(list(Line(line2)), "line2")
line3_ln <- Lines(list(Line(line3)), "line3")
line_sp <- SpatialLines(list(line1_ln, line2_ln, line3_ln))
line_spdf <- SpatialLinesDataFrame(line_sp, data.frame(c(1:3), 
row.names=c("line1", "line2", "line3")))
# 図2.2
plot(line_spdf, lty=1, lwd=3)
# 面データ
yoko <- readShapePoly("yoko.shp")
yoko_coord <- coordinates(yoko)
# 図2.3
plot(yoko, col="grey", border="black", lwd=3)
# 参考:自分で作る場合
poly1 <- cbind(c(0, 5, 5, 0, 0), c(0, 0, 5, 5, 0))
poly2 <- cbind(c(5, 10, 10, 5, 5), c(0, 0, 5, 5, 0))
poly3 <- cbind(c(0, 5, 5, 0, 0), c(5, 5, 10, 10, 5))
poly4 <- cbind(c(5, 10, 10, 5, 5), c(5, 5, 10, 10, 5))
poly1_pl <- Polygons(list(Polygon(poly1)), "poly1")
poly2_pl <- Polygons(list(Polygon(poly2)), "poly2")
poly3_pl <- Polygons(list(Polygon(poly3)), "poly3")
poly4_pl <- Polygons(list(Polygon(poly4)), "poly4")
poly_sp <- SpatialPolygons(list(poly1_pl, poly2_pl, poly3_pl, poly4_pl))
poly_spdf <- SpatialPolygonsDataFrame(poly_sp, data.frame(c(1:4), 
row.names=c("poly1", "poly2", "poly3", "poly4")))
plot(poly_spdf, col="grey") 
# par(mfrow=c(3,1))
# plot(pnt_spdf)
# plot(line_spdf)
# plot(poly_spdf, col="grey") 
# グリッドデータ
grid_topo <- GridTopology(cellcentre.offset=c(1, 1), 
cellsize=c(2, 2), cells.dim=c(5, 5)) 
grid_sgdf <- SpatialGridDataFrame(grid_topo, 
data=as.data.frame(runif(25, min=1, max=50)))
# image(grid_sgdf, col=grey((0:50)/51)) 
# グリッドデータの面オブジェクトへの変換
sp_grd <- as.SpatialPolygons.GridTopology(grid_topo) 
# 図2.4
plot(sp_grd, lwd=3)
text(getSpPPolygonsLabptSlots(sp_grd), 
getSpPPolygonsIDSlots(sp_grd)) 
sp_spdf<- SpatialPolygonsDataFrame(sp_grd, 
data=data.frame(c(1:25),
row.names=sapply(slot(sp_grd, "polygons"), 
function(i) slot(i, "ID")))) 
spplot(sp_spdf)
# ASCII grid
# library(maptools)
# x <- readAsciiGrid("fig2-6.txt")
# spplot(x)

2.2 単一レイヤでの操作

# 2.2.1 面積の計算
# SpatialPolygonsDataFrameを作成した際に、面積やセントロイドは計算済みである
poly3_gpc <- as(poly3, "gpc.poly")
area.poly(poly3_gpc)
# パッケージ geosphere を使ってでも計算できる。
# 2.2.2 セントロイド
library(geosphere)
plot(yoko, col="white", border="grey", lwd=3)
points(centroid(yoko), lwd=3, pch=19)
# plot(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)),
# type="l", main="", xlab="", ylab="")
# points(centroid(pol)) <- pol要確認
# 2.2.3 距離計算
# ユークリッド距離
# dist(pnt_xy)
# ポイントデータの座標
x <- c(5, 2, 6, 8,10)
y <- c(6, 4,10, 0, 7)
xy <- t(rbind(x,y))
# 作図
# plot(xy, xlim=c(0,10), ylim=c(0,10), cex=2)
# abline(h=0:10,v=0:10,lty=2, col="grey")
plot(xy,xlim=c(0,10),ylim=c(0,10),cex=2, pch=19)
abline(h=0:10,v=0:10,lty=2)
# ユークリッド距離
dist(xy, method="euclidean")
# マンハッタン距離
dist(xy, method="manhattan")
# 2.2.4 属性テーブルの結合
library(maptools)
jpn_pref <- readShapePoly("jpn_pref.shp",IDvar="PREF_CODE")
jpn_COD <- read.table("COD.csv",sep=",",header=TRUE)
ID.match <- match(jpn_pref$PREF_CODE, jpn_COD$PREF_CODE)
ID.match
jpn_COD1 <- jpn_COD[ID.match,]
summary(jpn_COD1)
jpn_pref_COD <- spCbind(jpn_pref,jpn_COD1)
names(jpn_pref_COD)
summary(jpn_pref_COD) 
# 2.2.5 ディゾルブ
library(maptools)
jpn <- readShapeSpatial("japan_ver62.shp", IDvar="JCODE",
proj4string=CRS("+proj=tmerc +ellps=GRS80 +units=m +datum=JGD2000"))
library(gpclib)
gpclibPermit()
jpn1 <- jpn
pref1 <- unionSpatialPolygons(jpn1, round(jpn1$JCODE / 1000))
# 2.2.6 切り取り
yoko <- jpn[jpn$JCODE > 14100,]
yoko <- yoko[yoko$JCODE < 14120,]
# writeSpatialShape(yoko, "yoko")
# 図2.14
plot(jpn, xlim=c(139.4, 139.8), ylim=c(35.02, 35.8), col="grey95", border="grey50",lwd=2)
plot(yoko, col="grey", border="black", add=TRUE, lwd=3)
# 2.2.7 バッファリング
# パッケージの読み込み
library(geosphere)
library(sp)
library(maptools) 
library(spatstat)
# 例1:任意の点
# 参照→http://r-sig-geo.2731867.n2.nabble.com/compute-buffer-from-point-shapefile-to-have-shapefile-td4574666.html
x <- c(5, 2, 6, 8,10)
y <- c(6, 4,10, 0, 7)
temp <- as.data.frame(cbind(x,y))
colnames(temp) <- c("x", "y")
polys<-list()
for(i in 1:nrow(temp)) {
 discbuff<-disc(radius=1, centre=c(temp$x[i], temp$y[i]))
    discpoly<-Polygon(rbind(cbind(discbuff$bdry[[1]]$x,
y=discbuff$bdry[[1]]$y), c(discbuff$bdry[[1]]$x[1],
y=discbuff$bdry[[1]]$y[1])))
    polys<-c(polys, discpoly)
}
spolys<-list()
for(i in 1:length(polys)) {
 spolybuff<-Polygons(list(polys[[i]]), ID=row.names(temp)[i])
 spolys<-c(spolys, spolybuff)
}
polybuffs<-SpatialPolygons(spolys)
# 図2.15
plot(polybuffs, col="grey", lwd=2)
points(temp, col="black", pch=20, cex=3) 
# 例2:任意の点(geoshpere)
x <- c(139.8836, 139.0608, 139.6489, 140.1233, 139.6917, 139.6425)
y <- c(36.56583, 36.39111, 35.85694, 35.60472, 35.68944, 35.44778)
coords <- cbind(x, y)
dist = 10000
angles = seq(1,360, by=5)
crds = list()
for (i in 1:nrow(coords)) {
       d = destPoint(coords[i,], angles, dist)
       crds[[i]] = rbind(d, d[1,])
}
p = lapply(crds, Polygon)
pp = list()
for (i in 1:length(p)) pp[i] = Polygons(p[i], i)
spdf = SpatialPolygonsDataFrame(SpatialPolygons(pp), data.frame(id=1:length(p)))
plot(coords, xlim=c(138.5,140.5), ylim=c(35.0,37.0))
plot(spdf, col=2, add=T)
points(coords, pch=20)
# 例3:都道府県庁ポイントデータ
# ポイントデータの座標
pref_gov <- read.table("pref_gov.txt",",",header=TRUE, row.names=2)
coords <- cbind(pref_gov$X,pref_gov$Y)
dist = 40000
angles = seq(1,360, by=5)
crds = list()
for (i in 1:nrow(coords)) {
       d = destPoint(coords[i,], angles, dist)
       crds[[i]] = rbind(d, d[1,])
}
p = lapply(crds, Polygon)
pp = list()
for (i in 1:length(p)) pp[i] = Polygons(p[i], i)
spdf = SpatialPolygonsDataFrame(SpatialPolygons(pp), data.frame(id=1:length(p)))
plot(spdf)
points(coords, pch=20)
# 2.2.8 ボロノイ分割
x <- c(5, 2, 6, 8,10)
y <- c(6, 4,10, 0, 7)
plot(deldir(x,y), axes=FALSE, cex=3)

2.3 複数レイヤでの操作

# 2.3.1 ポリゴンデータ内のポイントデータの属性の集計
overlay(pnt_spdf, sp_spdf)
overlay(pnt_spdf, sp_spdf, fn=colSums)
overlay(pnt_spdf, sp_spdf, fn=mean)
overlay(pnt_spdf, sp_spdf, fn=sd)
# overlay(pnt_spdf, poly_spdf)
# overlay(pnt_spdf, poly_spdf, fn=colSums)
# overlay(pnt_spdf, poly_spdf, fn=mean)
# overlay(pnt_spdf, poly_spdf, fn=sd)
# 2.3.2 面オブジェクトに対するオーバーレイ操作
library(gpclib)
# オーバーレイする面オブジェクト
p1 <- cbind(c(0, 5, 5, 0, 0), c(0, 0, 5, 5, 0))
p2 <- cbind(c(5, 7.5, 5, 2.5, 5), c(2.5, 5, 7.5, 5, 2.5))
# gpc.polyに変換
p1_gpc <- as(p1, "gpc.poly")
p2_gpc <- as(p2, "gpc.poly")
# 図2.19
plot(p1_gpc, poly.args = list(lwd=8), axes=FALSE, xlab="", ylab="")
plot(p2_gpc, poly.args = list(lwd=8), axes=FALSE, xlab="", ylab="")
#
plot(append.poly(p1_gpc, p2_gpc), poly.args = list(lwd=8), axes=FALSE,
xlab="", ylab="")
# erase -> setdiff()
poly_diff <- setdiff(p1_gpc, p2_gpc)
plot(append.poly(p1_gpc, p2_gpc), poly.args = list(lwd=8), axes=FALSE,
xlab="", ylab="")
plot(poly_diff, poly.args=list(col="grey", lwd=8), add=TRUE)
# union -> union()
poly_union <- union(p1_gpc, p2_gpc)
plot(append.poly(p1_gpc, p2_gpc), poly.args = list(lwd=8), axes=FALSE,
xlab="", ylab="")
plot(poly_union, poly.args=list(col="grey", lwd=8), add=TRUE)
# intersect -> intersect()
poly_intersect <- intersect(p1_gpc, p2_gpc)
plot(append.poly(p1_gpc, p2_gpc), poly.args = list(lwd=8), axes=FALSE,
xlab="", ylab="")
plot(poly_intersect, poly.args=list(col="grey", lwd=8), add=TRUE)
# identity -> append.poly()
poly_append <- append.poly(p1_gpc, p2_gpc)
plot(poly_append, poly.args=list(col="grey", lwd=8), axes=FALSE,
xlab="", ylab="")
# plot(append.poly(p1_gpc, p2_gpc), poly.args=list(col="grey"))
# update
poly_update <- append.poly(setdiff(p1_gpc, p2_gpc), intersect(p1_gpc, p2_gpc))
plot(append.poly(p1_gpc, p2_gpc), poly.args=list(lwd=8), axes=FALSE,
xlab="", ylab="")
plot(poly_update, poly.args=list(col="grey", lwd=8), add=TRUE)

演習問題1

yoko <- readShapePoly("yoko.shp")
yoko_coord <- coordinates(yoko)
plot(yoko, col="grey", border="white")
text(yoko_coord[,1], yoko_coord[,2], yoko$JCODE, cex=0.7, col="black")

演習問題2

px <- runif(100,0,10)
py <- runif(100,0,10) 
pz <- as.data.frame(px+py)
colnames(pz) <- c("pz")
pnt_xy <- cbind(px, py)
pnt_sp <- SpatialPoints(data.frame(px, py))
pnt_spdf <- SpatialPointsDataFrame(pnt_sp, pz)
plot(pnt_spdf, pch=20, cex=2, axes=TRUE, cex.axis=1.8)

備忘録

# データ読み込み
pnt <- read.table("chap02_1.csv", sep=",", header=TRUE)
# xy座標
pnt_xy <- cbind(pnt$x, pnt$y)
# ポイントデータ
# Spatial形式
pnt_sp <- SpatialPoints(data.frame(pnt$x, pnt$y))
bbox(pnt_sp)
# Spatial Data Frame 形式
pnt_spdf <- SpatialPointsDataFrame(pnt_xy, pnt)
head(pnt_spdf)
plot(pnt_spdf)
# ラインデータ
line1 <- cbind(c(64500, 64500, 67500, 67500, 64500), c(21000, 24000, 24000, 21000, 21000))
line2 <- cbind(c(67500, 67500, 70000, 70000, 67500), c(21000, 24000, 24500, 21500, 21000))
line3 <- cbind(c(70000, 70000, 72000, 73500, 73500, 70000), c(21500, 24500, 24500, 23000, 21000, 21500))
line1_ln <- Lines(list(Line(line1)), "line1")
line2_ln <- Lines(list(Line(line2)), "line2")
line3_ln <- Lines(list(Line(line3)), "line3")
line_sp <- SpatialLines(list(line1_ln, line2_ln, line3_ln))
line_spdf <- SpatialLinesDataFrame(line_sp, data.frame(c(1:3),  row.names=c("line1", "line2", "line3")))
plot(line_spdf)
# ポリゴンデータ
poly1 <- cbind(c(64500, 64500, 67500, 67500, 64500), c(21000, 24000, 24000, 21000, 21000))
poly2 <- cbind(c(67500, 67500, 70000, 70000, 67500), c(21000, 24000, 24500, 21500, 21000))
poly3 <- cbind(c(70000, 70000, 72000, 73500, 73500, 70000), c(21500, 24500, 24500, 23000, 21000, 21500))
poly1_pl <- Polygons(list(Polygon(poly1)), "poly1")
poly2_pl <- Polygons(list(Polygon(poly2)), "poly2")
poly3_pl <- Polygons(list(Polygon(poly3)), "poly3")
poly_sp <- SpatialPolygons(list(poly1_pl, poly2_pl, poly3_pl))
poly_spdf <- SpatialPolygonsDataFrame(poly_sp, data.frame(c(1:3), row.names=c("poly1", "poly2", "poly3")))
plot(poly_spdf, col="grey")
par(mfrow=c(3,1))
plot(pnt_spdf)
plot(line_spdf)
plot(poly_spdf, col="grey")
# plot(poly_spdf, col="grey")
# plot(line_spdf, add=T)
# グリッドデータ
grid_topo <- GridTopology(cellcentre.offset=c(64500, 21000), cellsize=c(500, 500), cells.dim=c(18, 6))
grid_topo
summary(grid_topo)
# grid_sgdf <- SpatialGridDataFrame(grid_topo, data=as.data.frame(c(1:108)))
grid_sgdf <- SpatialGridDataFrame(grid_topo, data=as.data.frame(runif(108, min=1, max=10)))
image(grid_sgdf, col=grey((0:30)/31))
# coordinates(grid_topo)
# grid_sp <- SpatialGrid(grid_topo)
# gridded(grid_sp) =TRUE
# grid_sgdf <- as(grid_sp, "SpatialGridDataFrame")
# plot(grid_sgdf)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-11-08 (土) 18:24:42 (1077d)