「An Introduction to R for Spatial Analysis and Mapping」の第5章をベースにRにおける地理情報の操作を学んでいく。

https://us.sagepub.com/en-us/nam/an-introduction-to-r-for-spatial-analysis-and-mapping/book241031

前回はこちらから。

http://rpubs.com/dichika/rsa_5_1

まずは前回の一連のオブジェクトを作成するところから再開する。

library("GISTools")
data(tornados)
index <- us_states$STATE_NAME %in% c("Texas", "New Mexico", "Oklahoma", "Arkansas")
AoI <- us_states[index, ]

前回の操作ではtornをAoIの範囲に限定したAoI.tornを作成し描画した。

AoI.tornはSpatialPointsクラスとなっており、座標の情報のみを格納したオブジェクトである。 このオブジェクトを属性情報をデータフレームの形で付加したSpatialPointsDataFrameに戻したい。

この操作を下記手順で進めることにする。

AoIとtornが対応づけられたrownameを取得し、AoIとtornのrownameに分割する

ここでgIntersection()においてbyid=TRUEを指定すると、AoI.tornのrownameにAoIとtornからそれぞれ抽出されたrownameがスペースで区切られて結合された形で格納される。 なお、byid=TRUEを指定すると実行にかなり長い時間がかかるので注意する。

system.time(AoI.torn <- gIntersection(AoI, torn, byid=TRUE))
##    user  system elapsed 
## 270.705   1.077 277.030
class(AoI.torn) # SpatialPointsクラス
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
head(data.frame(AoI.torn)) # rownameにAoIとtornのrownameが結合された形で格納されているのが確認できる
##             x     y
## 37 139 -97.60 35.55
## 37 140 -95.75 34.85
## 37 141 -97.02 35.82
## 37 142 -95.83 36.13
## 37 143 -99.28 34.88
## 37 144 -96.40 35.08

このrownameをstrsplit()を用いてAoIとtornのrownameに分割する。

なお、AoIのrownameはus_statesのrownameを示していることに注意する。

tmp <- rownames(data.frame(AoI.torn))
tmp <- strsplit(tmp, " ")

us_statesとtornからデータを抽出する

分割して得たAoIとtornのrownameを用いて、us_statesとtornからデータを抽出する。 抽出したデータはcbindで結合する。

# AoI.tornのrownamesを用いてus_statesからデータを抽出する
state.id <- sapply(tmp, "[[", 1)
df1 <- data.frame(us_states)[as.numeric(state.id), ]

# AoI.tornのrownamesを用いてtornからデータを抽出する
torn.id <- sapply(tmp, "[[", 2)
df2 <- data.frame(torn)[as.numeric(torn.id), ]

df <- cbind(df1, df2)

AoI.tornにAoIとtornの情報を付与しSpatialPointsDataFrameに変換する

SpatialPointsクラスであるAoI.tornにdata属性としてdfを与えてSpatialPointsDataFrameクラスに変換する。

AoI.torn_ <- SpatialPointsDataFrame(AoI.torn, data=df)
class(AoI.torn_)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"