lldbEC  <a href="http://xxrpafythekc.com/">xxrpafythekc</a>, [url=http://vexboyxdqjdt.com/]vexboyxdqjdt[/url], [link=http://qhsvdbvbxlah.com/]qhsvdbvbxlah[/link], http://topdvhhwjcas.com/
** 第2章プログラム [#qf7c2471]
*** 全体を通じて使用するパッケージ [#zee0a3c0]
 library(spdep)
 # library(sp)

*** 2.1 空間データの基本構造 [#effb48ad]
 # 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 単一レイヤでの操作 [#u52cd25d]
 # 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 複数レイヤでの操作 [#f791f7cd]
 # 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 [#bcba6a4a]
 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 [#aaf66e4e]
 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)

*** 備忘録 [#k03f256a]
 # データ読み込み
 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