地球儀と楕円体の描画
1. 地球儀を描く
「１．地球楕円体上の点」で、緯度経度（φ、λ）から（x,y,z）を求めたので地球儀を描いてみます。
RGL は、マウスで回転やズームができます。
- source("lib/DrawEarth.r")
- EARTH <- DrawEarth() # 地球の描画
- # mapdataパッケージのworld,riversを3Dで表示
- library(mapdata)
- d <- map("world",plot=F) # 行政区画(x:経度,y:緯度)
- d <- degXYZ(d$y,d$x,EARTH$b) # XYZ座標系に変換
- lines3d(d[,1],d[,2],d[,3],col="gray40",add=T)
- d <- map("rivers",plot=F) # 河川
- d <- degXYZ(d$y,d$x,EARTH$b) # XYZ座標系に変換
- lines3d(d[,1],d[,2],d[,3],col="paleturquoise",add=T)
|
地球の外観を描く方法はさておき、国境線は、パッケージ mapdata から、"world"と"rivers"という2組の座標データを参照して、lines3d() で結ぶだけです。
mapdata の図形は、緯度：y、経度：x で、単位は「度」です。
- r$> d <- map("world",plot=F) # 行政区画(x:経度,y:緯度)
- r$> str(d)
- List of 4
- $ x : num [1:82403] -69.9 -69.9 -69.9 -70 -70.1 ...
- $ y : num [1:82403] 12.5 12.4 12.4 12.5 12.5 ...
- $ range: num [1:4] -180 190.3 -85.2 83.6
- $ names: chr [1:1627] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda" ...
- - attr(*, "class")= chr "map"
- r$> range(na.omit(d$x))
- [1] -180.0000 190.2708
- r$> range(na.omit(d$y))
- [1] -85.19218 83.59961
|
x の範囲は [-180,190.2708 ] で、西経を負で表しています。y は、南緯が負になっているようです。
- # 弧度で与えられた経度と緯度（地理緯度）からXYZ座標を計算
- radXYZ <- function(φ,λ,b)
- {
- β <- atan(b*tan(φ))
- cbind(cos(β)*cos(λ),cos(β)*sin(λ),b*sin(β))
- }
- #------------------------------------------------------
- # 度で与えられた経度と緯度（地理緯度）からXYZ座標を計算
- degXYZ <- function(latitude,longitude, b) #!!! 経度、緯度の順
- { radXYZ(latitude*pi/180, longitude*pi/180, b) }
|
地球の外観がないと下図のように描画されます。地球の裏側を隠すために塗りつぶしています。
|
- library(mapdata)
- b <- 1 - 1/298.257222101
- d <- map("world",plot=F) # 行政区画(x:経度,y:緯度)
- d <- degXYZ(d$y,d$x,b) # XYZ座標系に変換
- lines3d(d[,1],d[,2],d[,3],col="gray40",add=T)
- d <- map("rivers",plot=F) # 河川
- d <- degXYZ(d$y,d$x,b) # XYZ座標系に変換
- lines3d(d[,1],d[,2],d[,3],col="paleturquoise",add=T)
|
|
2. 楕円体の描画
不透明な楕円体を描いて、地球の裏側の線を消します。この楕円を描く方法が分かりません。
始めた頃は、ellipsoid() で描けていたと思います。そのうちに未定義となったので、Rquarkパッケージをinstallして Rquark::rgl.ellipsoid() としていました。現在は、これも不可です。ellipsoid() の実態は cdaパッケージのものらしいのですが、cdaパッケージは、「使用しているバージョンのRにインストールできない」となります。
別途、ellipse3d() という関数があり、これも楕円体を描きます。
|
- clear3d()
- ee <- ellipse3d(diag(3),scale = c(1, 1, 1),
- centre = c(0, 0, 0),
- col="cyan",alpha=0.5)
- shade3d(ee)
- t <- seq(0,2*pi,length.out=360)
- lines3d(cos(t),sin(t),rep(0,360),col="red")
|
|
しかし、scale = c(1, 1, 1) は、半径1ではありません。scale = c(0.36, 0.36, 0.36) とすると半径が概ね1になるので代用できると思います。
しかし、理屈が分かりません。
結局、cdaパッケージから ellipsoid() のコードを借用することにします。そのまま写して使いますが、理解したいので、最小を求めてみます。
cdaパッケージから ellipsoid() の処理は概ね以下のようです。
| １．原点を中心に一辺が長さが2の立方体の図形を作る。8頂点は(-1,-1,-1)、(1,-1,-1)、(-1,1,-1)、・・・ r$> cube <- rgl::cube3d() shade3d(cube) str(cube) List of 6 $ vb : num [1:4, 1:8] -1 -1 -1 1 1 -1 -1 1 -1 1 ... $ material : list() $ normals : NULL $ texcoords: NULL $ meshColor: chr "vertices" $ ib : num [1:4, 1:6] 1 3 4 2 3 7 8 4 2 4 ... - attr(*, "class")= chr [1:2] "mesh3d" "shape3d" |
|
| ２．各辺の2分割を3回繰り返して頂点を増やす r$> sphere <- rgl::subdivision3d(cube, 3) shade3d(sphere) str(sphere) List of 3 $ material: list() $ ib : num [1:4, 1:384] 1 195 99 196 51 197 99 195 27 198 ... $ vb : num [1:4, 1:386] -15627 -15627 -15627 27171 15627 ... - attr(*, "class")= chr "mesh3d" |
|
| ３．各頂点の座標を正規化（原点との距離が1になる） r$> norm <- sqrt( sphere$vb[1, ]^2 + sphere$vb[2, ]^2 + sphere$vb[3, ]^2 ) for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm sphere$vb[4, ] <- 1 r$> str(sphere) List of 3 $ material: list() $ ib : num [1:4, 1:384] 1 195 99 196 51 197 99 195 27 198 ... $ vb : num [1:4, 1:386] -0.577 -0.577 -0.577 1 0.577 ... - attr(*, "class")= chr "mesh3d" |
|
|
|
原点から1の386頂点の多角形が描かれます。以下がXY平面を赤道面とした楕円体を描くのに必要なコードのようです。
- a <- 1; b <- 1; c <- 0.6
- sphere <- rgl::subdivision3d(rgl::cube3d(col="cyan"), 3)
- norm <- sqrt(sphere$vb[1, ]^2 + sphere$vb[2, ]^2 + sphere$vb[3, ]^2 )
- for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm
- sphere$vb[4, ] <- 1
- sphere$normals <- sphere$vb
- shade3d(rgl::scale3d(sphere, a,b,c))
|
3. 楕円体の回転、平行移動
cdaパッケージの ellipsoid() の描画は、中心が原点以外だったり、座標平面に対して傾いている楕円体も描きます。
cdaパッケージの ellipsoid() には以下のような続きがあります。
- rotM <- cpp_euler_passive(phi,theta,psi)
- result <- rgl::rotate3d(result,matrix=rotM)
- result <- rgl::translate3d(result, x,y,z)
- invisible(result)
|
3軸、それぞれの回転角 phi,theta,psi から、回転行列を求め回転します。その上で、原点が指定された x,y,z になるように平行移動しています。
4. 回転について
cdaパッケージはインストールできなかったのですが、utils.cpp の cpp_euler_passive()を見ることはできます。
各軸の周りの回転は、以下のように表されます。
cpp_euler_passive()の計算は、Z軸に回転(ψ)、Y軸に回転(θ)、Z軸に回転(φ)の順に回転することを示すと思います。
これは、私には理解できません。引数に、phi、theta、psi を指定するわけですが、どんな場合に、それが既知なのか想像できません。また、３軸ではなく、２軸の回転んで表すことの利点も分かりません。
私が必要性を感じるのは任意の直線を軸に回転することです。直線は２点で決まるので、２点を通る直線の周りの回転です。
