mikeo_410
  1. 2.地球
    1. 1.地球楕円体上の点
    2. 2.地球儀と楕円体の描画
    3. 3.3Dの直線と平面
    4. 4.卯酉線
    5. 5.楕円の標準形
    6. 6.5点で決まる楕円
    7. 7.曲率半径と平均曲率
    8. A1.地図、地球儀(Rスクリプト、データ)
    9. A2. 付録6.計算式集の公式

地球儀と楕円体の描画

1. 地球儀を描く

1.地球楕円体上の点」で、緯度経度(φ、λ)から(x,y,z)を求めたので地球儀を描いてみます。

RGL は、マウスで回転やズームができます。

  1. source("lib/DrawEarth.r")
  2. EARTH <- DrawEarth() # 地球の描画
  3. # mapdataパッケージのworld,riversを3Dで表示
  4. library(mapdata)
  5. d <- map("world",plot=F)       # 行政区画(x:経度,y:緯度)
  6. d <- degXYZ(d$y,d$x,EARTH$b)   # XYZ座標系に変換
  7. lines3d(d[,1],d[,2],d[,3],col="gray40",add=T)
  8. d <- map("rivers",plot=F)      # 河川
  9. d <- degXYZ(d$y,d$x,EARTH$b)   # XYZ座標系に変換
  10. lines3d(d[,1],d[,2],d[,3],col="paleturquoise",add=T)

地球の外観を描く方法はさておき、国境線は、パッケージ mapdata から、"world"と"rivers"という2組の座標データを参照して、lines3d() で結ぶだけです。

mapdata の図形は、緯度:y、経度:x で、単位は「度」です。

  1. r$> d <- map("world",plot=F) # 行政区画(x:経度,y:緯度)
  2. r$> str(d)
  3. List of 4
  4.  $ x    : num [1:82403] -69.9 -69.9 -69.9 -70 -70.1 ...
  5.  $ y    : num [1:82403] 12.5 12.4 12.4 12.5 12.5 ...
  6.  $ range: num [1:4] -180 190.3 -85.2 83.6
  7.  $ names: chr [1:1627] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda" ...
  8.  - attr(*, "class")= chr "map"
  9. r$> range(na.omit(d$x))
  10. [1] -180.0000  190.2708
  11. r$> range(na.omit(d$y))
  12. [1] -85.19218  83.59961

x の範囲は [-180,190.2708 ] で、西経を負で表しています。y は、南緯が負になっているようです。

  1. # 弧度で与えられた経度と緯度(地理緯度)からXYZ座標を計算
  2. radXYZ <- function(φ,λ,b)
  3. {
  4.     β <- atan(b*tan(φ))
  5.     cbind(cos(β)*cos(λ),cos(β)*sin(λ),b*sin(β))
  6. }
  7. #------------------------------------------------------
  8. # 度で与えられた経度と緯度(地理緯度)からXYZ座標を計算
  9. degXYZ <- function(latitude,longitude, b) #!!! 経度、緯度の順
  10. { radXYZ(latitude*pi/180, longitude*pi/180, b) }

地球の外観がないと下図のように描画されます。地球の裏側を隠すために塗りつぶしています。

  1. library(mapdata)
  2. b <- 1 - 1/298.257222101
  3. d <- map("world",plot=F)  # 行政区画(x:経度,y:緯度)
  4. d <- degXYZ(d$y,d$x,b)    # XYZ座標系に変換
  5. lines3d(d[,1],d[,2],d[,3],col="gray40",add=T)
  6. d <- map("rivers",plot=F) # 河川
  7. d <- degXYZ(d$y,d$x,b)    # XYZ座標系に変換
  8. lines3d(d[,1],d[,2],d[,3],col="paleturquoise",add=T)

2. 楕円体の描画

 不透明な楕円体を描いて、地球の裏側の線を消します。この楕円を描く方法が分かりません。

始めた頃は、ellipsoid() で描けていたと思います。そのうちに未定義となったので、Rquarkパッケージをinstallして Rquark::rgl.ellipsoid()  としていました。現在は、これも不可です。ellipsoid() の実態は cdaパッケージのものらしいのですが、cdaパッケージは、「使用しているバージョンのRにインストールできない」となります。

別途、ellipse3d() という関数があり、これも楕円体を描きます。

  1. clear3d()
  2. ee <- ellipse3d(diag(3),scale = c(1, 1, 1),
  3.                 centre = c(0, 0, 0),
  4.                 col="cyan",alpha=0.5)
  5. shade3d(ee)
  6. t <- seq(0,2*pi,length.out=360)
  7. 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() の処理は概ね以下のようです。

1.原点を中心に一辺が長さが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.各辺の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"

3.各頂点の座標を正規化(原点との距離が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平面を赤道面とした楕円体を描くのに必要なコードのようです。

  1. a <- 1; b <- 1; c <- 0.6
  2. sphere <- rgl::subdivision3d(rgl::cube3d(col="cyan"), 3)
  3. norm <- sqrt(sphere$vb[1, ]^2 + sphere$vb[2, ]^2 + sphere$vb[3, ]^2 )
  4. for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm
  5. sphere$vb[4, ] <- 1
  6. sphere$normals <- sphere$vb
  7. shade3d(rgl::scale3d(sphere, a,b,c))

3. 楕円体の回転、平行移動

cdaパッケージの ellipsoid() の描画は、中心が原点以外だったり、座標平面に対して傾いている楕円体も描きます。

cdaパッケージの ellipsoid() には以下のような続きがあります。

  1.   rotM <- cpp_euler_passive(phi,theta,psi)
  2.   result <- rgl::rotate3d(result,matrix=rotM)
  3.   result <- rgl::translate3d(result, x,y,z)
  4.   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 を指定するわけですが、どんな場合に、それが既知なのか想像できません。また、3軸ではなく、2軸の回転んで表すことの利点も分かりません。

私が必要性を感じるのは任意の直線を軸に回転することです。直線は2点で決まるので、2点を通る直線の周りの回転です。


題目一覧へmikeo_410@hotmail.com(updated: 2022/09/09)