演習の概要

  • 空間データの基本構造を理解する
  • 空間レイヤーの操作（単一レイヤー）

注意事項

  • 慶應義塾大学SFCで開講している「環境ガバナンスのデータサイエンス」「空間モデリング特論」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.5.3を使っています。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

Rの利用方法

  • フリー統計ソフトRを利用するには、Rをダウンロードする、Webブラウザ上で利用する、などの方法があります。
  • RはCRANからダウンロードし、インストール可能です。
  • RStudioには、より操作が容易なGUIが用意されています。
  • RStudio Cloudを使うと、RをダウンロードすることなくWebブラウザ上でRを操作できます。
  • RStudioは日本語が使えない場合があります。起動時にLC_CTYPE failedなどのエラーメッセージが出る場合は、以下の一行を実行しRStudioを再起動すると、日本語が使えるようになることがあります。
system("defaults write org.R-project.R force.LANG en_US.UTF-8")

データ分析の準備

データのダウンロード

演習で用いるデータはここからダウンロードしてください。

この演習では~/直下に作業用のgisdataフォルダを作成し、そこに演習用データがある想定としています。

この場合、setwd("~/gisdata")とすることで、ディレクトリを指定できます。

パッケージのインストール

install.packages()でパッケージをインストールし、library()でパッケージを呼び出すことができます。

この演習全体を通じて、以下のようなパッケージを用います。データ分析を行うにあたっては、適宜パッケージのインストールと読み込みを行ってください。

install.packages("spdep")
install.packages("sp")
install.packages("maptools")
install.packages("rgdal")
install.packages("sf")
install.packages("classInt")
install.packages("RColorBrewer")
install.packages("readr")
install.packages("magrittr")
install.packages("dplyr")
install.packages("kokudosuuchi")
install.packages("estatapi")
install.packages("purrr")
install.packages("needs")
install.packages("DT")
install.packages("gpclib")
install.packages("geosphere")
install.packages("deldir")
# install.packages("devtools")
devtools::install_github("tidyverse/ggplot2")
library(spdep)
library(sp)
library(maptools)
library(rgdal)
library(sf)
library(classInt)
library(RColorBrewer)
library(readr)
library(magrittr)
library(dplyr)
library(kokudosuuchi)
library(estatapi)
library(purrr)
library(needs)
library(DT)
library(gpclib)
library(geosphere)
library(deldir)
# library(devtools)
library(ggplot2)

以下の方法でパイプ %>% の優先順位を高めることができます。パイプの使い方は授業内などで説明します。

needs::prioritize(magrittr)

空間データの基本構造

空間オブジェクトを作成する

ランダムなポイントデータを作成してみよう

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)
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")))
plot(line_spdf, lty=1, 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")

グリッドデータを作成する

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) 
plot(sp_grd, lwd=3)
text(getSpPPolygonsLabptSlots(sp_grd), 
getSpPPolygonsIDSlots(sp_grd)) 
## Warning: use coordinates method
## Warning: use *apply and slot directly

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)

既存データを読み込む

横浜市のシェープファイル（ポリゴン）を読み込み描画する

yoko <- sf::st_read("yoko.shp")
yoko %>%
  st_geometry() %>%
  plot(border="white", col="grey")

単一レイヤの操作

セントロイドの作成

セントロイドはいくつかのパッケージで作成することができますが、ここではgeosphereパッケージのcentroid()関数を使う方法を紹介します。

ggplot2::ggplot()+
  ggplot2::geom_sf(data = yoko)+
  ggplot2::geom_sf(data = st_centroid(yoko))

距離計算（ユークリッド距離とマンハッタン距離）

dist()関数を使うと、距離を計算できます。ここでは、ユークリッド距離とマンハッタン距離の計算例を示します。

# ポイントデータの座標
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, pch=19)
abline(h=0:10,v=0:10,lty=2)

# ユークリッド距離
dist(xy, method="euclidean")
##           1         2         3         4
## 2  3.605551                              
## 3  4.123106  7.211103                    
## 4  6.708204  7.211103 10.198039          
## 5  5.099020  8.544004  5.000000  7.280110
# マンハッタン距離
dist(xy, method="manhattan")
##    1  2  3  4
## 2  5         
## 3  5 10      
## 4  9 10 12   
## 5  6 11  7  9

バッファリング

ここでは、緯度経度で座標が与えられた点オブジェクトに対して、バッファを作成する例を示します。

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)

ボロノイ分割

ボロノイ分割は、deldirパッケージのdeldier()関数を使って作成できます。

library(deldir)
## deldir 0.1-16
x <- c(5, 2, 6, 8,10)
y <- c(6, 4,10, 0, 7)
plot(deldir(x,y), axes=FALSE, cex=3)
## 
##      PLEASE NOTE:  The components "delsgs" and "summary" of the
##  object returned by deldir() are now DATA FRAMES rather than
##  matrices (as they were prior to release 0.0-18).
##  See help("deldir").
##  
##      PLEASE NOTE: The process that deldir() uses for determining
##  duplicated points has changed from that used in version
##  0.0-9 of this package (and previously). See help("deldir").