Rによるボロノイグリッド・ボロノイ図の作成、シェープファイルの出力

こんにちは、エンヂニアです。

今日は初のR関連の記事を投稿してみたいと思います。

第一弾はVoronoi grid(ヴォロノイグリッド)の作成方法とそのプロットの仕方を自分自身の備忘録を兼ねて投稿したいと思います。

Rはversion 3.4.3を使用しており、RStudio (Version 1.0.136)を使ってスクリプトを作成・実行しています。

背景

まずは簡単に本投稿の背景から。

資源・エネルギー開発Engineerであるエンヂニアは、あるエリア・地域に広く・しばしば疎に散らばった空間データを扱うことが頻繁にあります。

そこで、それぞれのデータポイントが地域的にどのような分布をしているのか、色々なVisualisationをしてみたくなるわけです。

今回取り上げるVoronoi gridはそのようなVisualisationの一つとなります。

ここでは、以下のステップを踏んで行きたいと思います。

  1. 適当にデータを作成
  2. ヴォロノイグリッドプロット作成
  3. シェープファイルの作成・出力

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での、カレントフォルダの変更方法は以下を参考にしてみて下さい。

Rの醍醐味は、大量のデータを簡単なコードで読み込ませ、計算させ、描画できちゃう点です。これには、未だに感動します。 そのためには、データを読みこませるステップが必要になります。さらにその前に、今Rのコンソールがどのディレクトリ(フォ...

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)

一行目のカラーパレットの設定に関しては、こちらの記事を参考にしてみて下さい。

こんにちは、エンヂニアです。 今日はRにおける、colorspaceパッケージを用いたインタラクティブなカラーの設定について簡...

続いて、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社のウェブサイトをご覧下さい。

シェープファイル(Shapefile)とは、Esri 社の提唱したベクトル形式の業界標準フォーマットです。Esri 製品はもちろん、多くの GIS ソフトウェアで利用が可能です。 シェープファイルは仕様が公開されています。シェープファイルの読み込み、書き込みに対応した独自のアプリケーションを開発されたい場合は仕様書をご...

まずは、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)

参考サイト

  1. https://www.rdocumentation.org/packages/latticeExtra/versions/0.6-28/topics/tileplot
  2. https://www.rdocumentation.org/packages/latticeExtra/versions/0.6-28/topics/panel.voronoi

今日も最後までお付き合い頂きありがとうございます。↓ポチッとクリックして下さったら尚感謝です。

にほんブログ村 海外生活ブログ オーストラリア情報へ

にほんブログ村

シェアする

  • このエントリーをはてなブックマークに追加

フォローする