mikeo_410
  1. ３．R
    1. １．Rの常套句
    2. ２．いろいろ（R）
    3. ３．座標の引数（ベクトル、matrix、bind、list）
    4. ４．plot
    5. ５．RGL
    6. ６．ラグランジュの未定乗数法(Rsolnp)
    7. ７．常微分方程式の汎用解法（deSolve）

座標の引数（ベクトル、matrix、bind、list）

points(x,y) は、points(1,2)のようにも、points(c(1,2),c(3,4)) のようにも記述できます。

引数のx,yは、それぞれ、ベクトルでも、行列でも、リストでも良いように出来ています。同じように関数を作りたいと考えました。

1. Rの配列、リスト

C++やC#のようなプログラミング言語では、型を作って、その配列を扱うことができます。しかし、Rでは配列として扱えるのは基本的な型だけのようです。逆に、C++やC# で型を作っても組み込みの演算子は作用しませんが、Rの演算子はベクトルを演算するように出来ています。

  1. > a <- as.numeric(1); b <- as.integer(2); c <- as.double(3)
  2. > d <- as.complex(4); e <- as.raw(5)
  3. > f <- as.logical(6); g <- as.character(7)
  4. > sprintf("typeof() a:%s b:%s c:%s d:%s e:%s f:%s g:%s",
  5. +     typeof(a),typeof(b),typeof(c),typeof(d),typeof(e),typeof(f),typeof(g))
  6. [1] "typeof() a:double b:integer c:double d:complex e:raw f:logical g:character"
  7. > sprintf("class() a:%s b:%s c:%s d:%s e:%s f:%s g:%s",
  8. +     class(a),class(b),class(c),class(d),class(e),class(f),class(g))
  9. [1] "class() a:numeric b:integer c:numeric d:complex e:raw f:logical g:character"
  10. > sprintf("mode() a:%s b:%s c:%s d:%s e:%s f:%s g:%s",
  11. +     mode(a),mode(b),mode(c),mode(d),mode(e),mode(f),mode(g))
  12. [1] "mode() a:numeric b:numeric c:numeric d:complex e:raw f:logical g:character"

配列として扱うとは、array[i] のようにインデクスで参照する構文が使えることを考えています。array()は、基本的な型の配列を作ります。基本的な型のオブジェクト以外は、as.vector()されるとあるので、array()が作る配列の要素はベクトルだと考えられ、ベクトルの要素は基本の型です。

ここではmode()が表す数値型（numeric）を主に考えます。

Arrayの訳語が「配列」だと考えて、インデクスで参照する構文で扱えるものの総称に使います。

Rのmatrix()が作るオブジェクトは配列で、他の方法で作成したベクトルや配列と共に、行列計算の演算子や関数の入力になります。行列は数学関数などの説明で転置行列のように使われているようです。

Rにある配列を作る関数は、以下のようです。

関数

説明

array

Multi-way Arrays

c

Combine Values into a Vector or List

cbind

Combine R Objects by Columns

rbind

Combine R Objects by Rows

matrix

Matrices

vector

Vectors - Creation, Coercion, etc

character

  1. > character(length=10)
  2.  [1] "" "" "" "" "" "" "" "" "" ""
  3. > logical(length=10)  
  4.  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  5. > numeric(length=10)
  6.  [1] 0 0 0 0 0 0 0 0 0 0

logical

numeric

ここでの、ベクトル、配列、行列には、単に同じ型のものが並んでいてインデクスで要素にアクセスすると言う以外の意味はありません。c(1,"a") のように書けば、最初の要素の型も character になって格納され、同じ型のものが並ぶことになります。

我流に解釈すれば、Rでは演算子が作用する対象はベクトルで、配列も準じたものなのです。ベクトルは主に基本型の列で、配列はベクトルをバインドしたもので、矩形に並びます。

list() や data.frame() もベクトルをバインドする機能だと思います。バインドするベクトルの型が異なっていても良い点で上述の配列とは異なっています。また、上述の配列と data.frame は矩形で、行、列の長さは、それぞれ全て同じですが、リストの要素は、それぞれ独立したベクトルで個別の長さです。

  1. > vector1 <- c(1,2,3)
  2. > vector2 <- c("4","5","6")
  3. > (u <- data.frame(vector1,vector2))
  4.   vector1 vector2
  5. 1       1       4
  6. 2       2       5
  7. 3       3       6
  8. > u[,1];u[,2]
  9. [1] 1 2 3
  10. [1] "4" "5" "6"
  11. > vector1 <- c(1,2)
  12. > vector2 <- c("4","5","6")
  13. > data.frame(vector1,vector2)
  14.  data.frame(vector1, vector2) でエラー:
  15.    引数に異なる列数のデータフレームが含まれています: 2, 3
  16. > list(vector1,vector2)
  17. [[1]]
  18. [1] 1 2
  19. [[2]]
  20. [1] "4" "5" "6"
2. vector について

Rの演算子が作用する対象がベクトルなのだと解釈すると、ベクトルの要素は同じ型で一律に処理できることを前提にしていると考えられます。

入れ物と中身と考えると、リストがコンテナで、中身はベクトルです。list() や data.frame() の行ベクトルか列ベクトルの一方が元のベクトルで、他方は要素ごとに型が異なる可能性があることになります。

まず、行と列を確認しておきます。

  1. > (u <- matrix(1:6, nrow=2,ncol=3))
  2.      [,1] [,2] [,3]
  3. [1,]    1    3    5
  4. [2,]    2    4    6
  5. > u[1,]
  6. [1] 1 3 5
  7. > u[,1]
  8. [1] 1 2

2行、3列の matrix で確認すると、[1,] が水平方向の3つの値を示し、行ベクトルです。[,1]は垂直方向の2つの値を示し、列ベクトルです。list や data.frame も同様に[行、列]とインデクスすると考えます。

  1. > vector1 <- c(1,2)
  2. > vector2 <- c("4","5","6")
  3. > u<-list(vector1,vector2)
  4. > str(u[1])
  5. List of 1
  6.  $ : num [1:2] 1 2
  7. > str(u[[1]])
  8.  num [1:2] 1 2
  9. > (vector <- u[[1]])
  10. [1] 1 2

試してみると、list は2次元の配列としてアクセスできないようです。u[1] は最初の要素を list で返し、u[[1]] が元のベクトルを返します。

data.frame の場合は、ベクトルは「列」として格納されています。「行」を要求すると data.frame 型で返されます。

  1. > x <- c(1,2,3)
  2. > y <- c("4","5","6")
  3. > (u <- data.frame(x,y))
  4.   x y
  5. 1 1 4
  6. 2 2 5
  7. 3 3 6
  8. > str(u[,1])
  9.  num [1:3] 1 2 3
  10. > str(u[,2])
  11.  chr [1:3] "4" "5" "6"
  12. > str(u[1,])
  13. 'data.frame':   1 obs. of  2 variables:
  14.  $ x: num 1
  15.  $ y: chr "4"
3. 配列のデータの並び

t <- (0:360)*(pi/180)

x <- cos(t); y <- sin(t)

plot(x,y,t=”l”)

と書けば円を描けます。

  1. > m <- cbind(x,y)
  2. > str(m)
  3.  num [1:361, 1:2] 1 1 0.999 0.999 0.998 ...
  4.  - attr(*, "dimnames")=List of 2
  5.   ..$ : NULL
  6.   ..$ : chr [1:2] "x" "y"
  7. > plot(m,t="l")
  8. > n <- rbind(x,y)
  9. > str(n)
  10.  num [1:2, 1:361] 1 0 0.9998 0.0175 0.9994 ...
  11.  - attr(*, "dimnames")=List of 2
  12.   ..$ : chr [1:2] "x" "y"
  13.   ..$ : NULL
  14. > plot(n,t="l")

cbind(x,y) は、xがm[,1]です。rbindだとxはn[1,]です。mをplotすると円、nをplotすると直線が描かれます。

  1. > u <- matrix(c(x,y),ncol=2)
  2. > str(u)
  3.  num [1:361, 1:2] 1 1 0.999 0.999 0.998 ...
  4. > plot(u,t="l")
  5. > v <- matrix(c(x,y),nrow=2)
  6. > str(v,t="l")
  7.  num [1:2, 1:361] 1 1 0.999 0.999 0.998 ...
  8. > plot(v,t="l")

matrix の場合、uをplotすると円、vをplotすると直線が描かれます。

plot などのグラフ描画関数は、xy座標を配列で受け取りますが、[,1] が x、[,2] が y と決まっています。

配列は「列」をバインドして出来ていると考えるのが普通のようです。

しかし、配列の要素は全て同じ型なのが普通と考えられるので、通常は、列ベクトルと行ベクトルは対等に扱われるものと思います。以下は、ソートの例です。

  1. > x <- c(1,3,2,4)
  2. > y <- c("a","c","b","d")
  3. > u <- cbind(x,y)
  4. > str(u)
  5.  chr [1:4, 1:2] "1" "3" "2" "4" "a" "c" "b" "d"
  6.  - attr(*, "dimnames")=List of 2
  7.   ..$ : NULL
  8.   ..$ : chr [1:2] "x" "y"
  9. > u[order(u[,1]),]
  10.      x   y
  11. [1,] "1" "a"
  12. [2,] "2" "b"
  13. [3,] "3" "c"
  14. [4,] "4" "d"
  15. > v <- rbind(x,y)
  16. > str(v)
  17.  chr [1:2, 1:4] "1" "a" "3" "c" "2" "b" "4" "d"
  18.  - attr(*, "dimnames")=List of 2
  19.   ..$ : chr [1:2] "x" "y"
  20.   ..$ : NULL
  21. > v[,order(v[1,])]
  22.   [,1] [,2] [,3] [,4]
  23. x "1"  "2"  "3"  "4"
  24. y "a"  "b"  "c"  "d"

1列目、あるいは1行目で配列をソートしています。

多くの場合、配列はランダムにアクセスするので、メモリ上での連続性を意識する必要はなく、どちらの並びでも間違えなければ良いことです。特別な理由がなければ「普通」にしておけば良さそうです。

4. matrixのbyrow

数列 c(1,2,3, 4,5,6) があって、意味的に x <- c(1,2,3) 、y <- c(4,5,6) と分けられるとします。

  1. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 2)
  2.      [,1] [,2] [,3]
  3. [1,]    1    3    5
  4. [2,]    2    4    6
  5. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 3)
  6.      [,1] [,2]
  7. [1,]    1    4
  8. [2,]    2    5
  9. [3,]    3    6
  10. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 2)
  11.      [,1] [,2]
  12. [1,]    1    4
  13. [2,]    2    5
  14. [3,]    3    6
  15. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 3)
  16.      [,1] [,2] [,3]
  17. [1,]    1    3    5
  18. [2,]    2    4    6

matrix()の nrow、ncol では、x、y が 行になるような配列（3行2列）しか作れません。

byrow=TRUE とすると、x、y が 列になるような配列（2行3列）しか作れません。

  1. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, byrow=TRUE)
  2.      [,1] [,2] [,3]
  3. [1,]    1    2    3
  4. [2,]    4    5    6
  5. > matrix(c(1, 2, 3, 4, 5, 6), nrow = 3, byrow=TRUE)
  6.      [,1] [,2]
  7. [1,]    1    2
  8. [2,]    3    4
  9. [3,]    5    6
  10. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 2, byrow=TRUE)
  11.      [,1] [,2]
  12. [1,]    1    2
  13. [2,]    3    4
  14. [3,]    5    6
  15. > matrix(c(1, 2, 3, 4, 5, 6), ncol = 3, byrow=TRUE)
  16.      [,1] [,2] [,3]
  17. [1,]    1    2    3
  18. [2,]    4    5    6

plot などのグラフ描画関数は、列ベクトルが、それぞれ x、y なので、2つのベクトル x、y を c(x,y) と書いて matrix を作るとすると byrow=TRUE は必須です。

5. 座標を表す配列引数

plot()やpoints()などの関数は、引数が配列である場合、[,1]をx、[,2]をyと見なすと考えられます。

配列は基本的な型で作られ、(x,y)のような組の配列を作るようには考えられていないようです。

[,1]はxの値のベクトルで、[,2]はyの値のベクトルです。

座標を表す配列の引数は、xのベクトルとyのベクトルをバインドしたものです。

  1. x<-0; y<-0
  2. plot(x,y,col="white",xlim=c(-5,5),ylim=c(-5,5))
  3. points(0,0,pch=16,col="red")
  4. points(c(1,2,3),c(3,2,1),pch=16,col="blue")
  5. u <- cbind(c(1,2,3),c(-3,-2,-1))
  6. points(u,pch=16,col="green")
  7. u <- matrix(c(-3,-2,-1, 1,2,3),ncol=2)
  8. points(u,pch=16,col="magenta")
  9. points(list(x=c(-1,-2,-3),y=c(1,2,3))) # OK
  10. points(list(y=c(-1,-2,-3),x=c(1,2,3))) # x,yを逆に書いてもOK

この例は、意図通り描画されます。

以下は、意図通りには描かれません。

  1. u <- rbind(c(1,2,3),c(-3,-2,-1))
  2. #(u[1,1]=1,u[1,2]=2),(u[2,1]=-3,u[2,2]=-2)の
  3. points(u,pch=16,col="green")    # 2点のみ描画
  4. u <- matrix(c(-3,-2,-1, 1,2,3),nrow=2,byrow=TRUE)
  5. points(u,pch=16,col="green")    # 2点のみ描画

上手くいかないのは、行と列が逆になっているからです。plot()やpoints()は、作成した側の意図とは無関係に決まった順でベクトルを参照します。

しかし、list()の場合は、ラベルには意味があります。

6. リストとラベル

リストの場合は、要素のベクトルの名前x、yが有効に見えます。

  1. > points(list(a=c(-1,-2,-3),b=c(1,2,3)))
  2.  xy.coords(x, y) でエラー: 
  3.    'x' はリストですが，成分 'x' と 'y' を持ちません

確かに、リストのラベルはplot()やpoints()関数で認識されるラベルです。

配列にも行名、列名を付けられます。

  1. > x<-0; y<-0
  2. > plot(x,y,col="white",xlim=c(-5,5),ylim=c(-5,5))
  3. > x<-c(-1,-2,-3); y<-c(1,2,3)
  4. > u <- cbind(x,y)
  5. > u
  6.       x y
  7. [1,] -1 1
  8. [2,] -2 2
  9. [3,] -3 3
  10. > points(u,pch=16,col="red")
  11. > u <- cbind(y,x)
  12. > u
  13.      y  x
  14. [1,] 1 -1
  15. [2,] 2 -2
  16. [3,] 3 -3
  17. > points(u,pch=16,col="blue")

cbind() の引数の順番を変えれば、値と共にラベルも入れ替わります。異なる図になるのは、ラベルに無関係に順番によって図の x、y が決まっているからだと考えられます。

matrixを使って、値が入れ替わらないように、ラベルを y、x の順に付与してみます。

  1. > x<-0; y<-0
  2. > plot(x,y,col="white",xlim=c(-5,5),ylim=c(-5,5))
  3. > x<-c(-1,-2,-3); y<-c(1,2,3)
  4. > u <- matrix(c(x,y),ncol=2)
  5. > colnames(u) <- c("x","y")
  6. > u
  7.       x y
  8. [1,] -1 1
  9. [2,] -2 2
  10. [3,] -3 3
  11. > points(u,pch=16,col="red")
  12. > colnames(u) <- c("y","x")
  13. > u
  14.       y x
  15. [1,] -1 1
  16. [2,] -2 2
  17. [3,] -3 3
  18. > points(u,col="blue")

ラベルを入れ替えても、同じ位置に表示され、ラベルは使われないようです。

7. 座標を示す引数の扱い
7.1. 座標を示す引数

原点に点を表示するには、points(0,0,pch=20) などと書きます。

points()などのplot関連の関数の座標は、plot(x,y) のように指定します。x、y は、それぞれ同じ長さのベクトルです。

R は、x<-0 のように単一の値の場合も x[1] と記述して、xだけを書いた場合と同じに扱われます。

plot(0) のように、引数x、yは省略可能です。

引数が1つの場合、yの座標値が与えられたものと解釈されます。xについては、1からの連番が生成されます。x軸はyの値のインデクスになります。したがって、plot(c(0,0)) とすると、(1,0)、(2,0) と2つの点が打たれることになります。

ただし、引数がベクトルではない配列やリストのときには、1つの配列やリストがx、yの両方を含んでいると見なされます。

これは、c(x,y) のようには考えないことを示しています。1つのベクトルに x と y のペアを記録するようには考えず、それぞれ個別のベクトルにあるようになっています。

7.2. x、y の省略

plot() などの引数 x、y の省略は便利です。しかし、自身が作成する関数で同じようにするのが便利なケースはないと思います。省略が有効なケースがあっても、それが 1 からの連番であることはないと思います。

7.3. 1点の座標

lines(c(1,2),c(3,4)) は、(1,3) と (2,4) の間に線を描きます。c(1,2) と c(3,4) の間ではありません。

p1 <- c(1,2)

p2 <- c(3,4)

lines(p1,p2)とは行かず、

lines(c(p1[1],p2[1]), c(p1[2],p2[2])

と記す必要があります。

p1 <- cbind(1,2) のように書けば、x、yの2つのベクトルとなって一貫性がありますが、lines() が lines(p1,p2) と書けるようになる訳ではありません。

自身が作成するプログラムでは、2つの値を持つベクトルを (x,y) と見なすことは避けようがないと思います。

8. 引数や戻り値を揃える

作成する関数の戻り値については、list(x=x,y=y,・・・) のように決めます。そのままplot()などの入力になり、座標以外の値を含めることができます。

引数も、これを受容することになります。

XYZ座標から（緯度、経度）を計算するプログラムで考えてみます。緯度は南緯を負で表し、西経を360°から180°で表します。ただし、弧度です。

すべて、赤道半径を1として計算しています。

  1. XYZ2φλ <- function(x,y,z, f=0){
  2.   b <- 1-f
  3.   r <- sqrt(x^2+y^2)
  4.   λ <- atan2(y,x)
  5.   λ <- ifelse(abs(λ)<(16*.Machine$double.eps), 0, λ)
  6.   λ <- ifelse(λ<0, 2*pi+λ, λ)
  7.   θ <- atan2(z,r) # 中心角
  8.   φ <- atan(tan(θ)/b^2)
  9.   list(φ=φ, λ=λ)
  10. }

このプログラムは、引数x、y、z が、それぞれ1つの値のときと、それぞれが同じ長さのベクトルのときに機能します。

配列やリストを引数には出来ないので、呼び出し元で、配列やリストの中の対応するベクトルに分けて引数にしています。

配列やリストを引数にするには、以下のようなことが必要です。

  1. 配列やリストのときは引数の y、z が省略されていなければならない。
  2. 配列なら、x<-[,1]、y<-[,2]、x<-[,3] とベクトルを取り出す。
  3. リスト、data.frame なら、x<-$x、y<-$y、x<-$z とベクトルを取り出す。

これらを考慮して書き直してみます。

  1. XYZ2φλ <- function(x,y,z, f=0){
  2.   b <- 1-f
  3.   if(missing(y) | missing(z)){
  4.     if(missing(y) & missing(z)){
  5.         if(is.list(x)){ y <- x$y; z <- x$z; x <- x$x
  6.         } else if(is.array(x)){ y <- x[,2]; z <- x[,3]; x <- x[,1] }
  7.         else {warning("引数のy,zは省略されているがxはリストでも配列でもない")}
  8.     } else {warning("引数が揃っていない")}
  9.   }
  10.   r <- sqrt(x^2+y^2)
  11.   λ <- atan2(y,x)
  12.   λ <- ifelse(abs(λ)<(16*.Machine$double.eps), 0, λ)
  13.   λ <- ifelse(λ<0, 2*pi+λ, λ)
  14.   θ <- atan2(z,r) # 中心角
  15.   φ <- atan(tan(θ)/b^2)
  16.   list(φ=φ, λ=λ)
  17. }

これで、いろいろな形式の引数に対応できました。

  1. > x <- c(cos(pi/6),cos(pi/4),cos(pi/3));
  2. > y <- c(0,0,0);
  3. > z <- c(sin(pi/6),sin(pi/4),sin(pi/3))
  4. > str(XYZ2φλ(x[1],y[1],z[1]))
  5. List of 2
  6.  $ φ: num 0.524
  7.  $ λ: num 0
  8. > str(XYZ2φλ(x,y,z))
  9. List of 2
  10.  $ φ: num [1:3] 0.524 0.785 1.047
  11.  $ λ: num [1:3] 0 0 0
  12. > str(XYZ2φλ(cbind(x,y,z)))
  13. List of 2
  14.  $ φ: num [1:3] 0.524 0.785 1.047
  15.  $ λ: num [1:3] 0 0 0
  16. > str(XYZ2φλ(matrix(c(x,y,z),ncol=3)))
  17. List of 2
  18.  $ φ: num [1:3] 0.524 0.785 1.047
  19.  $ λ: num [1:3] 0 0 0
  20. > str(XYZ2φλ(list(x=x,y=y,z=z)))
  21. List of 2
  22.  $ φ: num [1:3] 0.524 0.785 1.047
  23.  $ λ: num [1:3] 0 0 0
  24. > str(XYZ2φλ(data.frame(x,y,z)))
  25. List of 2
  26.  $ φ: num [1:3] 0.524 0.785 1.047
  27.  $ λ: num [1:3] 0 0 0

要点は、data.frame も is.list() にTRUEが返り、matrix()やcbind()で作った配列もis.array()にTRUEを返すことです。

また、エラーが起きるケースのうち、この関数で決めた引数のルール以外は標準のエラーになることで良いと考えました。

逆向きの変換は以下のように書き直しました。

  1. φλ2XYZ <- function(φ,λ, f=0){
  2.   b <- 1-f
  3.   if(missing(λ)){
  4.     if(is.list(φ)){
  5.         if(exists("φ",where=φ) & exists("λ",where=φ)){
  6.             λ <- φ$λ; φ <- φ$φ
  7.         } else if(exists("x",where=φ) & exists("y",where=φ)){
  8.             λ <- φ$x; φ <- φ$y
  9.         } else {warning("リストには(φ,λ)か(x,y)が必要")}
  10.     } else if(is.array(x)){ λ <- φ[,2]; φ <- φ[,1] }
  11.     else {warning("引数のλは省略されているがφはリストでも配列でもない")}
  12.   }
  13.   β <- atan2(b*sin(φ),cos(φ))
  14.   list(x=cos(β)*cos(λ),y=cos(β)*sin(λ),z=b*sin(β))
  15. }
9. 常用系との変換

常用系と言っても文字表記ではなく、経度を180°、緯度を90°で表すと言うことです。

経度は東経を正、西経を負で0°から180°で表します。緯度は、北緯を正、南緯を負で 0°から90°で表します。

内部は弧度を使用し、経度は0から2πで東回りで表します。緯度は、赤道が0で北にπ/2、南にーπ/2 まで表します。

  1. lon2λ <- function(lon){
  2.   (lon %% 360) /180*pi
  3. }

  1. λ2lon <- function(λ){
  2.   u <- (λ %% (2*pi)) /pi*180
  3.   i <- which(u > 180)
  4.   u[i] <- u[i] - 360
  5. }

  1. lat2φ <- function(lat){
  2.   lat /180*pi
  3. }

  1. φ2lat <- function(φ){
  2.   φ /pi*180
  3. }

10. 航程線を描いて見る

いくらか簡潔になったか航程線を引いてみます。経線は極で1点に集まりますが、メルカトル図法の地図は経線を平行に描きます。したがって、極に近づくほど左右に拡大されます。その率は緯度をφとして 1/cosφ です。これと同じだけ上下にも拡大したものがメルカトル図法の地図のようです。このメルカトル図法の地図で2点間に直線を引いたものを航程線と考え球面に移します。

まず、球面と平面に地図を描いて、2地点に赤い点を描きます。

  1. clear3d(); rgl.light();
  2. ee <- ellipsoid(col="white",alpha=0.9); shade3d(ee)
  3. d <- map("world2Hires",plot=FALSE)
  4. # 3D
  5. u <- φλ2XYZ(lat2φ(d$y), lon2λ(d$x))
  6. lines3d(u)
  7. YO <- cbind(lat2φ(35.45033),  lon2λ(139.63422))
  8. VA <- cbind(lat2φ(49.266667), lon2λ(-123.116667))
  9. yo <- φλ2XYZ(YO)
  10. va <- φλ2XYZ(VA)
  11. p <- ellipsoid(yo, a=0.02,b=0.02,c=0.02,col="red"); shade3d(p)
  12. p <- ellipsoid(va, a=0.02,b=0.02,c=0.02,col="red"); shade3d(p)
  13. # 2D
  14. MercatorY <- function(φ){ log(tan(pi/4+φ/2)) }
  15. Mercatorφ <- function(y){ 2*atan(exp(y))-pi/2 }
  16. x <- d$x/180*pi
  17. y <- MercatorY(lat2φ(d$y))
  18. par(mar=c(3,3,1,1),mgp=c(1.5, 0.5, 0))
  19. plot(x,y,t="l")
  20. x1 <- YO[2]; y1 <- MercatorY(YO[1])
  21. x2 <- VA[2]; y2 <- MercatorY(VA[1])
  22. points(x1,y1,pch=16,col="red")
  23. points(x2,y2,pch=16,col="red")

平面に2点を通る直線を描きます。直線の描画に使った点列の座標値を緯度、経度に直して、球面に描きます。

  1. # 2点を通る直線の傾きと切片を求めて直線を引く
  2. a <- (y2-y1)/(x2-x1); b <- (x2*y1-x1*y2)/(x2-x1)
  3. x <- seq(0,2*pi,length.out=361)
  4. y <- a*x+b
  5. lines(x,y,col="blue")
  6. # 直線を球面に描く
  7. λ <- x
  8. φ <- Mercatorφ(y)
  9. u <- φλ2XYZ(φ, λ)
  10. lines3d(u,col="blue")

map()が返す経度は、最初から360°表示でした。ただし、0°付近は負の値や360°以上になっている箇所があります。

角度は0°から360°未満に丸めたい訳ですが、Rの %% 演算子による剰余はこの目的に都合良く使えます。

剰余はプログラミング言語によって値が異なっています。

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