LoginSignup
0
0

More than 3 years have passed since last update.

@kaz0414(Kaz)

Rで３次元の色付き等高線図を描く

Last updated at Posted at 2020-12-26

要約

RStudioで以下のような3次元の色付き等高線図の描き方を紹介します。
20201226-1.png

  • 等高線は、数式で与えられているものとします。上の図の場合は、点$(-1,-1,-1)$を中心とする球の式（$(x+1)^2+(y+1)^2+(z+1)^2$）で描画しています。

なお、一般的に「3次元の等高線図」というのは、下図のようなものを呼びますが、本記事では、上に示した図、強いて言い直せば「3次元空間内の等高面図」のことを指すこととします。
20201226-5.png
※上図の参考Webサイト

環境

  • OS: Windows10
  • R: R x64 4.0.3
  • RStudio: Version 1.3.1093

説明

まずは、ソースコード（Rファイル）とその結果をお見せします。後ほど、説明をいたします。

ソースコード

test1.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2+1)^2+(arg3+1)^2)
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)

valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}

lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)

xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

scatterplot3d(xyz,color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z")

結果

20201226-1.png

解説

3次元の等高線図を描く関数というものは、（私の知る限り）存在しません。
今回は、3次元の散布図を描画するscatterplot3d()という関数で、無理やり3次元等高線図を描いています。

scatterplot3d()は、n個の空間内の点の座標を示すn行3列のデータを渡すと、3次元の散布図を描画してくれる関数です。
こちらのWebサイトの例を用いると、以下のソースコード

example1.r
scatterplot3d(iris[,1:3])

で、以下の3次元散布図を描画できます。
20201226-2.png

ここで用いたirisとは、Rにはじめから組み込まれているデータセットで、中身は以下のような感じです。
20201226-3.png
この例では、このデータセットの最初の3列を用いて、3次元散布図を描いています。

今回は、このscatterplot3d()を用いて、数式の値を元に色を付けた点を、ひたすら空間内に敷き詰めた散布図を作ることで、3次元等高線図を作成します。

ところで、点を敷き詰めるのは、グラフで見えている3つの面だけで十分です。見えてないところに敷き詰めてもいいのですが、描画処理時間と、グラフを外部出力するなら出力ファイル容量の無駄です。今回のソースコードでも、見えている3つの面にのみ、色付きの点を敷き詰めています。

一方、色を付ける際には、敷き詰める点についてだけ考えるのでは不十分です。
等高線図の色をつけるためには、グラフ描画範囲における値の範囲、すなわち、最大値と最小値が必要です。値の範囲がわかることで、各値について、どの色で塗ればいいのかがわかります。
ところが、点を敷き詰める3つの面に、最大値と最小値が両方ともあるとは限りません。なので、値の範囲を調べるために、値についてはグラフ描画範囲全体について調べる必要があります。

以上を踏まえたうえで、ソースコードを見てみましょう。

test1-1.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2+1)^2+(arg3+1)^2)
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)

formuraは、等高線に用いる数式です。値の範囲を調べるときと、点に色をつけるときの2回、この数式を用いるので、はじめに関数化しておきます。

scatterplot3d()を使うためには、library(scatterplot3d)でライブラリを呼び出す必要があります。なお、scatterplot3dがインストールされていない場合は、別途、インストールする必要があります。インストールの際には、ユーザに管理者権限がないとエラーが出ますのでご注意ください。

n <- 301で、1行あたりに敷き詰める点の個数を設定します。ここでは301で設定しています。

x <- seq(-1,1,length=n)、y <- seq(-1,1,length=n)、z <- seq(-1,1,length=n)で、3軸それぞれの配列を作成します。

test1-2.r
valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}

さて、ここでは、値の範囲（最大値・最小値）を計算しています。グラフ描画範囲全体（n×n×n）を走査して、値の最大値と最小値を得ます。
はじめに作成したformula()に、x[i]、y[j]、z[k]を渡して、値を取得しています。
なお、最小値の初期値にInf（正の無限大）を、最大値の初期値に-Inf（負の無限大）を、それぞれ設定しています。Rでは無限大を変数に代入することが可能です。

test1-3.r
lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)

lvsには、等高線の水準（塗り分ける色の境界値）を代入しています。pretty()は、与えられた配列に対して、「いい感じの」水準配列を返してくれる関数です。このときのpretty(c(valmin,valmax),10)は、

example2.r
> pretty(c(valmin,valmax),10)
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12

のように、0から12までの13個の整数の配列を返しています。
なお、pretty(c(valmin,valmax),10)の第二引数10は、水準数の「目安」を指定しています。あくまで目安なので、この数ぴったりの水準数になるとは限らない、という点に注意です。（今回も水準数は13個になっている）
nlvsには、lvsの長さを代入しています。

test1-4.r
xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

xyzには、描画する点の座標をセットします。1面あたりの点の数はn×nで、それが3面分あるので、xyzは3*n^2行3列の行列としています。初期値は0です。

二重for文の中では、具体的に点の座標を代入しています。
最初の3行はxy平面の、次の3行はzx平面、最後の3行はyz平面の点をx座標、y座標、z座標を代入しています。
意味合い的には以下のような3つの二重for文に分けた方が分かりやすいですが、処理時間短縮のために、1つの二重for文で済ましています。

example3.r
#xy平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
  }
}

#zx平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
  }
}

#yz平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

さて、次は点につける色の配列を作成します。

test1-5.r
clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

scatterplot3d(xyz,color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z")

clrsは、scatterplot3d()の引数colorに渡す配列です。この配列は、描画する点の数と同じ長さの配列である必要があり、対応する点の色をこの配列で指定することができます。

具体例を示します。

example4.r
scatterplot3d(rbind(c(1,1,1),c(2,2,2),c(3,3,3)),color=c("#FF0000","#00FF00","#0000FF"),pch=19)

このようなコードを実行すると、以下のような3次元散布図が描画されます。
20201226-4.png
1つ目の点c(1,1,1)には#FF0000（赤）を、2つ目の点c(2,2,2)には#00FF00（緑）を、3つ目の点c(3,3,3)には#0000FF（青）を、それぞれ対応させています。

このようにして、各点の色を指定する配列を作成します。今回は、それがclrsです。

clrsを作成するにあたり、その色サンプルとしてclrを用意します。
clrには、cm.colors(nlvs-1)をセットします。

for文の中では、formula()関数から得た値が、各水準のどこに当てはまるのかをif文で判定し、Trueとなったところの色をclrsにセットしています。

最後のscatterplot3d()では、点の座標データxyzと、各点の色指定配列clrs、そして、マーカーとしてドット.を指定しています。また、各軸ラベルもセットしています。

マーカーとしてドット.を指定しているのは、最も小さな点だからです。pch=20が、もともと組み込まれているマーカーの中では最も小さな点ですが、これでもまだ大きいので、ドット.を指定しています。

おまけ（y軸反転）

「Rの３次元散布図（scatterplot3d()）での軸の反転方法」に記載の方法を使うことで、y軸を反転させて、右手系座標で描画することもできます。

ソースコード

test2.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2-1)^2+(arg3+1)^2) #中心を(-1,1,-1)に変更
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)

valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}

lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)

xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- max(y) #ここをmin(y)からmax(y)に変更
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

xyz[,2] <- -xyz[,2] #y成分反転

scatterplot3d(xyz[,1],xyz[,2],xyz[,3],color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z",
              y.ticklabs=seq(max(xyz[,2]), min(xyz[,2]), -0.5)
              )

結果

20201226-6.png

姉妹記事

Rで色付き等高線図を描く

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0