要約

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

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

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

※上図の参考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")

結果

解説

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

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

example1.r

scatterplot3d(iris[,1:3])

で、以下の3次元散布図を描画できます。

ここで用いたirisとは、Rにはじめから組み込まれているデータセットで、中身は以下のような感じです。

この例では、このデータセットの最初の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次元散布図が描画されます。

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)
              )

結果

姉妹記事

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