こんにちは、エンヂニアです。
今日は初のR関連の記事を投稿してみたいと思います。
第一弾はVoronoi grid（ヴォロノイグリッド）の作成方法とそのプロットの仕方を自分自身の備忘録を兼ねて投稿したいと思います。
Rはversion 3.4.3を使用しており、RStudio (Version 1.0.136)を使ってスクリプトを作成・実行しています。
背景
まずは簡単に本投稿の背景から。
資源・エネルギー開発Engineerであるエンヂニアは、あるエリア・地域に広く・しばしば疎に散らばった空間データを扱うことが頻繁にあります。
そこで、それぞれのデータポイントが地域的にどのような分布をしているのか、色々なVisualisationをしてみたくなるわけです。
今回取り上げるVoronoi gridはそのようなVisualisationの一つとなります。
ここでは、以下のステップを踏んで行きたいと思います。
- 適当にデータを作成
- ヴォロノイグリッドプロット作成
- シェープファイルの作成・出力
1. データの作成
まずは適当に空間データを作成します。
### 必要となるパッケージをインストール
require("latticeExtra")
require("deldir")
require("colorspace")
require("tripack")
require("shapefiles")
### 必要となるパッケージを読み込み
library(latticeExtra)
library(deldir)
library(colorspace)
library(tripack)
library(shapefiles)
### アウトプットファイル名の指定
Output.file.voro = "Voronoi Coords.csv"
Output.file.area = "Voronoi Areas.csv"
Output.file.shape = "Voronoi"
### ランダムなデータポイントを100点用意
### コンテキストは空間データであるため、以下のような意味合い
### x: 東西位置, y: 南北位置, z: 各地点におけるある物理量の測定値
num.data <- 100
x <- rnorm(num.data)
y <- rnorm(num.data)
z <- y*2 + rnorm(num.data)
### x, y, zをデータフレームdata.dfに統合
data.df <- data.frame(data.num = seq(1, 100),
data.id = paste("D", seq(1, 100), sep = ""),
X = x, Y = y, Z = z)
### まずはプロットしてデータを理解
plot(data.df)
アウトプットファイルは、後で使いますが、カレントフォルダに直接出力するようにしています。RStudioでの、カレントフォルダの変更方法は以下を参考にしてみて下さい。
2. ヴォロノイグリッドプロット作成
それでは早速ヴォロノイグリッドのプロットを作成してみます。まずは、latticeExtraパッケージのtileplot関数を使ってみます。
x及びyを座標としデータポイントを平面図にプロットし、その周辺にヴォロノイグリッドを作成、各データポイントでのz値によってグリッドの各セルに色をつけます。
### latticeExtra::tilepolot関数によるヴォロノイグリッドプロット
mypal <- diverge_hcl(n = 50, h = c(168, 330), c = 100, l = c(19, 92),
power = 1.5, fixup = TRUE, gamma = NULL, alpha = 1)
tileplot(z ~ x * y, data.df, pch = 16, cex = 1, border = 1,
col = "black", col.regions = mypal)
一行目のカラーパレットの設定に関しては、こちらの記事を参考にしてみて下さい。
続いて、latticeパッケージのlevelplot関数を使ってみます。
### lattice::levelplot関数によるヴォロノイグリッドプロット
levelplot(z ~ x * y, data.df, panel = panel.voronoi, points = TRUE, border = 1,
col = "black", col.regions = mypal, pch = 16, cex = 1)
全く同じであることが分かります。
3. シェープファイルの作成・出力
最後に、shapefilesパッケージを使った、シェープファイル (ESRI Shapefile)の作成・出力にも触れてみたいと思います。
シェイプファイルは、ESRI社が考案したファイル形式で、様々なGISソフトウェアにおいて利用されている、地図や図形を平面図にプロットする為の形式です。
エンヂニアはSpotfireでよく利用しています。詳しくはESRI社のウェブサイトをご覧下さい。
まずは、tripackパッケージを使って、ヴォロノイグリッドのポリゴンを作成し、プロットしてみます。
### tripackパッケージを用いたヴォロノイグリッドポリゴン作成及びそのプロット tri.vm <- voronoi.mosaic(x, y) plot(x, y, pch = 16, col = "black") plot.voronoi(tri.vm, pch = 1, cex = 0, add = TRUE)
色以外、再度全く結果が得られました。
続いて、いくつか結果を見る為、また、その他ソフトウェアパッケージへの読み込みの為のcsvファイルへの出力と共に、シェープファイルを出力してみます。
### voronoiクラスであるtri.vmをvoronoi.polygonsクラスに変換
voro.polygon <- voronoi.polygons(tri.vm)
num.voro.polygon <- nrow(summary(voro.polygon))
### voronoi polygonsの面積を計算し、元のデータフレームにデータ行として追加
voro.area <- data.frame(area = voronoi.area(tri.vm))
data.df.withArea <- cbind(data.df, voro.area)
head(data.df.withArea)
### 面積データを追加したデータセットをcsvファイルに出力
write.csv(data.df.withArea, Output.file.area, row.names = FALSE)
### voronoi gridの座標: "voro.coords"を生成
voro.coords <- NULL
j <- 1
i <- 1
if(is.na(voro.area[i, 1])){
#Do nothing
}else{
tmp.df <- data.frame(voro.polygon[j])
tmp.df <- cbind(tmp.df, poly.id = c(i))
voro.coords <- rbind(voro.coords, tmp.df)
j <- j + 1
}
for(i in 2:num.data){
if(is.na(voro.area[i, 1])){
#Do nothing
}else{
tmp.df <- data.frame(voro.polygon[j])
tmp.df <- cbind(tmp.df, poly.id = c(i))
voro.coords <- rbind(voro.coords, tmp.df)
j <- j + 1
}
}
voro.coords <- voro.coords[c("poly.id", "x", "y")] #行順の入れ替え
##### voronoi grid座標をcsvファイルに出力
write.csv(voro.coords, Output.file.voro, row.names = FALSE)
####################################################
############# Generate ESRI Shapefiles #############
####################################################
##### data.frame型オブジェクトをshapefile型オブジェクトに変換
ddTable <- data.df.withArea[!(is.na(data.df.withArea$area)), c("data.num", "data.id")]
voro.shapefile <- convert.to.shapefile(voro.coords, ddTable , "data.num", 5)
##### shapefileを出力
write.shapefile(voro.shapefile, Output.file.shape, arcgis = T)
参考サイト
- https://www.rdocumentation.org/packages/latticeExtra/versions/0.6-28/topics/tileplot
- https://www.rdocumentation.org/packages/latticeExtra/versions/0.6-28/topics/panel.voronoi
今日も最後までお付き合い頂きありがとうございます。↓ポチッとクリックして下さったら尚感謝です。
