「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

コードは原則として書籍内のコードに準じた記載としている。

まず、本書の付属パッケージであるGISToolsパッケージをロードする。

本パッケージには地理データのサンプルおよびその他GISのデータ処理に便利な関数が含まれている。

library("GISTools")

data()によってtornadosをロードするとtorn, torn2, us_states, us_states2の4つのデータが読み込まれる。 今回用いるus_statesデータはSpatialPolygonsDataFrameクラスのデータである。 またtornデータはSpatialPointsDataFrameクラスのデータである。

data(tornados)

構造やクラスを確認する際はstr()を用いると表示が冗長になるので代わりにsummary()を用いる。 summary()の結果もかなり長いため、ここでは割愛する。

summary(us_states)
summary(torn)

SpatialPolygonsDataFrameクラスおよびSpatialPointsDataFrameクラスのオブジェクトはS4クラスシステムで記述されている。

isS4(us_states)

## [1] TRUE

isS4(torn)

## [1] TRUE

アメリカ全土における竜巻の発生地点を描画する

par(mar=c(0,0,0,0))
plot(us_states) # アメリカの地図が表示される
plot(torn, add=TRUE, pch=1, col="orange", cex=0.4) # 上記表示範囲内の竜巻の発生地点が描画される
plot(us_states, add=TRUE) # 竜巻の発生地点を描画した際に州の境界が消えてしまうので、あらためて描画する

４州における竜巻の発生地点を描画する

まずは４州を指定するインデックスを作成する

index <- us_states$STATE_NAME == "Texas"|us_states$STATE_NAME == "New Mexico"|us_states$STATE_NAME == "Oklahoma"|us_states$STATE_NAME == "Arkansas"

これは以下のようにも指定できる

index <- us_states$STATE_NAME %in% c("Texas", "New Mexico", "Oklahoma", "Arkansas")

４州をArea of Interest（AoI）として抽出する。

AoI <- us_states[index, ]

いったん４州を描画してみよう。

par(mar=c(0,0,0,0))
plot(AoI)

さて、ここでtornの描画を重ねるとAoI以外の範囲の竜巻の発生地点も描画されてしまう。

したがって、AoIの範囲内にtornを限定したい。

この際、rgeosパッケージのgIntersection()を用いる。

この関数は地理情報版のintersection()であり、2つのオブジェクトの共通範囲のデータを抽出する。

AoI.torn <- gIntersection(AoI, torn)

これでAoIの範囲内の竜巻発生地点が格納されたAoI.tornが作成できた。 AoIで境界を描画した上で、AoI.tornで竜巻の発生地点を重ねる。

par(mar=c(0,0,0,0))
plot(AoI)
plot(AoI.torn, add=TRUE, pch=1, col="orange", cex=0.4)

ggplot2を用いた描画

ここまでplot()を用いて描画してきたが、ggplot2による描画も試してみる。

アメリカ全土における竜巻の発生地点を描画する

ggplot2で描画するには原則としてデータフレームである必要がある。

SpatialPolygonDataFrameは受け入れてくれるが、SpatialPointsDataFrameは受け入れてくれないため、後者からデータフレームを抜き出す必要がある。

S4クラスシステムのオブジェクトの内部構造にアクセスする際は@を用いる。

@dataでデータフレームとして保持されているデータにアクセスできるので、SpatialPointsDataFrameからデータフレームを抜き出す。

library("ggplot2")
ggplot(us_states) + geom_polygon(fill="white", color="black", aes(x=long, y=lat, group=group)) + 
  geom_point(data=torn@data, aes(x=coords.x1, y=coords.x2), size=0.1, color="orange", alpha=0.3)

４州における竜巻の発生地点を描画する

gIntersectionを適用して作られたAoI.tornはSpatialPointsDataFrameクラスからSpatialPointsクラスに変換されている。

ggplot2で描画するにはいったんデータフレームに変換する必要がある。

broomパッケージのtidyで変換したいところだがSpatialPointsに対してはメソッドが用意されていないため、as.data.frameを用いる。

変換したAoI.torn_dfにおいては緯度経度を表す列名がx, yなので注意。

AoI.torn_df <- as.data.frame(AoI.torn)
ggplot(AoI) + geom_polygon(fill="white", color="black", aes(x=long, y=lat, group=group)) + 
  geom_point(data=AoI.torn_df, aes(x=x, y=y), size=0.1, color="orange", alpha=0.3)

sfパッケージによる記述

余談だが、最近開発が進んでいるsfパッケージを用いるとdplyrの関数で抽出ができる。

普段dplyrパッケージによるデータ操作に慣れ親しんでいる場合はsfパッケージの利用も考えてみてほしい。

library("sf")
library("dplyr")
us_state_sf <- st_as_sf(us_states)
AoI.us_state_sf <- us_state_sf %>% filter(STATE_NAME %in% c("Texas", "New Mexico", "Oklahoma", "Arkansas"))
plot(AoI.us_state_sf$geometry) # 境界のみを表示したい際は、plot()時にgeometryを指定する。

本当はsfパッケージを用いて全て書き換えたかったのだが、gIntersectionに相当するst_intersectionの実行速度があまりに遅いので今回は見送った。