R

2013年7月21日 (日)

『文字列に含まれる改行文字(\n)を除く』 などのRメモ

●文字列に含まれる改行文字(\n)を除く
例1
> text1 <- "あいうえを\nかきくけこ\nさしすせそ"
> text2  <- gsub('\\n', " ", text1)
> text3  <- gsub('\\n', "", text1)
> text1
[1] "あいうえを\nかきくけこ\nさしすせそ"
> text2
[1] "あいうえを かきくけこ さしすせそ"
> text3
[1] "あいうえをかきくけこさしすせそ"

例2
> text1="あいうえを\nかきくけこ\nさしすせそ"
> text2=unlist(strsplit(text1,"\n"))
> text3<- paste(text2,collapse = " ") #ベクトルの各要素を一つの文字列にする
> text4<- paste(text2,collapse = "")  #ベクトルの各要素を一つの文字列にする
> text1
[1] "あいうえを\nかきくけこ\nさしすせそ"
> text2
[1] "あいうえを" "かきくけこ" "さしすせそ"
> text3
[1] "あいうえを かきくけこ さしすせそ"
> text4
[1] "あいうえをかきくけこさしすせそ"

● ベクトルの各要素を一つの文字列にする
> a <- c(1,2,3,4)
> a
[1] 1 2 3 4
> paste(a,collapse=" ")
[1] "1 2 3 4"
> paste(a,collapse="")
[1] "1234"

> b <- paste(1:4)
> paste(b,collapse=" ")
[1] "1 2 3 4"

●参考
『11. データファイルの順次処理:ファイル名のいろいろな決め方』
http://takenaka-akio.org/doc/r_auto/chapter_11.html
『17. 文字列を操作する』
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/17.html
『R Programming/Text Processing』
http://en.wikibooks.org/wiki/R_Programming/Text_Processing
『R における正規表現』
http://www.okada.jp.org/RWiki/?R%20%A4%CB%A4%AA%A4%B1%A4%EB%C0%B5%B5%AC%C9%BD%B8%BD

●区切り文字で文字列を分割する
> mojiretu <- "$y=ax^2+bx+c$"
> s <- strsplit(mojiretu,"")[[1]]
> s
[1] "$" "y" "=" "a" "x" "^" "2" "+" "b" "x" "+" "c" "$"

●文字列から空白(半角スペース)を除く
> st <- "$\\abc   \\def  {10}  {20}  ghi$"
> sub(" ","",st)
[1] "$\\abc  \\def  {10}  {20}  ghi$"
> gsub(" ","",st)
[1] "$\\abc\\def{10}{20}ghi$"

●文字列の最初が『\\』にマッチする要素の位置番号とその部分の実際の文字列を返す.
> a=c("123","\\abc","\\def")
> grep("^\\\\",a)
[1] 2 3

> grep("^\\\\",a,value=TRUE)
[1] "\\abc" "\\def"

> b=c("123","abc","def")
> grep("^\\\\",b)
integer(0)

●文字列の最後が『\\』にマッチする要素の位置番号と
その部分の実際の文字列を返す.
> s1=c("123\\","abc\\","\\def")
> grep("\\\\$",s1)
[1] 1 2
> grep("\\\\$",s1,value=TRUE)
[1] "123\\" "abc\\"

●文字列の最初が『$』にマッチする要素の位置番号と
その部分の実際の文字列を返す.
> s2=c("$sin(x)$" , "$cos(x)" , "$tan(x)$" , "$(x+y)/z$")
> grep("^\\$",s2)
[1] 1 2 3 4
> grep("^\\$",s2,value=TRUE)
[1] "$sin(x)$"  "$cos(x)"   "$tan(x)$"  "$(x+y)/z$"

●文字列の最後が『$』にマッチする要素の位置番号と
その部分の実際の文字列を返す.
> s2=c("$sin(x)$" , "$cos(x)" , "$tan(x)$" , "$(x+y)/z$")
> grep("\\$$",s2)
[1] 1 3 4
> grep("\\$$",s2,value=TRUE)
[1] "$sin(x)$"  "$tan(x)$"  "$(x+y)/z$"

●文字列の最初と最後が『$』にマッチする要素の位置番号と
その部分の実際の文字列を返す.
> s2=c("$sin(x)$" , "$cos(x)" , "$tan(x)$" , "$(x+y)/z$")
> pattern="^\\$.*\\$$"
> grep(pattern,s2) 
[1] 1 3 4
> grep(pattern,s2,value=TRUE) 
[1] "$sin(x)$"  "$tan(x)$"  "$(x+y)/z$"

【例題1】
> s3=c("$sin(x)$" , "$cos(x)" , "$tan(x)$" , "$(x+y)/z$","$$")
> pattern="^\\$.*\\$$"
> grep(pattern,s3,value=TRUE) 
[1] "$sin(x)$"  "$tan(x)$"  "$(x+y)/z$" "$$"      

【例題2】
> s3=c("$sin(x)$" , "$cos(x)" , "$tan(x)$" , "$(x+y)/z$","$$")
> pattern="^\\$.+\\$$"
> grep(pattern,s3,value=TRUE) 
[1] "$sin(x)$"  "$tan(x)$"  "$(x+y)/z$"

●参考
・Pattern Matching and Replacement(R Documentation)
http://stat.ethz.ch/R-manual/R-devel/library/base/html/grep.html
・Quick-R
http://www.statmethods.net/management/functions.html

●all() と any()関数の使用例
> x<- 0:5
> any(x>3)
[1] TRUE
> x>3
[1] FALSE FALSE FALSE FALSE  TRUE  TRUE
> all(x>6)
[1] FALSE
> x>6
[1] FALSE FALSE FALSE FALSE FALSE FALSE
> all(x>=0)
[1] TRUE
> (x>=0)
[1] TRUE TRUE TRUE TRUE TRUE TRUE

【例題】 x<- 0:5とする.
・ all(x>0)
[1] FALSE

・ any(x<0)
[1] FALSE

・ x<4のtrueの数をsum(x<4)で求めてみよう.
> sum(x<4)
[1] 4
> x<4
[1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE

【メモ】
> b=x<4
> b
[1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE
> as.numeric(b)
[1] 1 1 1 1 0 0
> count=as.numeric(b)
> sum(count)
[1] 4

●データ型の変換
> x <- c(TRUE, FALSE)
> x
[1]  TRUE FALSE
> class(x)
[1] "logical"
> as.numeric(x)
[1] 1 0
> as.logical(c(-3 , -2 , -1 , 0 , 1 ,  2 ,  3 ,  4.2 , 2.57))
[1]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE

● 『&&、|| 』  『かつ、または』
> FALSE && TRUE #二つ目の式(TRUE)は評価されない
[1] FALSE
> FALSE &&  (print(1))  #二つ目の式(print(1))は評価されない
[1] FALSE
> TRUE &&  (print(1))  #二つ目の式(print(1))が評価されて初めて条件式の真偽の評価ができる
[1] 1
[1] TRUE

> TRUE || FALSE #二つ目の式FALSE)は評価されない
[1] TRUE
> TRUE  ||  print(1) #二つ目の式(print(0))は評価されない
[1] TRUE
> TRUE  ||  print(0) #二つ目の式(print(0))は評価されない
[1] TRUE
> FALSE ||  print(0) #二つ目の式(print(0))が評価されて初めて条件式の真偽の評価ができる
[1] 0
[1] FALSE
> FALSE ||  print(1) #二つ目の式(print(0))が評価されて初めて条件式の真偽の評価ができる
[1] 1
[1] TRUE

●ベクトルの要素ごとの演算結果
> c(TRUE, TRUE) & c(FALSE, TRUE) #ベクトルの要素ごとの演算
[1] FALSE  TRUE

> c(TRUE, TRUE) && c(FALSE, TRUE)
[1] FALSE
#最初のベクトルの第1要素(TRUE)と次のベクトルの第1要素(FALSE)の演算だけが行われる

> c(FALSE, FALSE) | c(FALSE, TRUE) #ベクトルの要素ごとの演算
[1] FALSE  TRUE

> c(TRUE, FALSE) | c(FALSE, FALSE) #ベクトルの要素ごとの演算
[1]  TRUE FALSE

> c(FALSE, FALSE) || c(FALSE, TRUE) #第1要素の演算だけが行われる
[1] FALSE
#最初のベクトルの第1要素(FALSE)と次のベクトルの第1要素(FALSE)の演算だけが行われる

> c(FALSE, FALSE) || c(TRUE, FALSE)
[1] TRUE
#最初のベクトルの第1要素(FALSE)と次のベクトルの第1要素(TRUE)の演算だけが行われる


2013年4月15日 (月)

日時の取得/時間をファイル名に/対話的入力/ディレクトリ内にあるファイル読込み/など

■日時の取得
・Sys.time {base}  Get Current Date and Time
http://stat.ethz.ch/R-manual/R-patched/library/base/html/Sys.time.html
・trptime {base}   Date-time Conversion Functions to and from Character
http://stat.ethz.ch/R-manual/R-patched/library/base/html/strptime.html

%x  Date. Locale-specific on output, "%y/%m/%d" on input.
%X Time. Locale-specific on output, "%H:%M:%S" on input.

■時間をファイル名に利用

> Sys.time()
[1] "2013-03-02 16:53:40 JST"
> Sys.Date()
[1] "2013-03-02"

> nichiji=paste(Sys.Date(), gsub(":", "-", format(Sys.time(), "%X")),sep = "_")
> nichiji
[1] "2013-03-02_16-27-15"
> fname=paste("file",nichiji,".csv",sep="")
> fname
[1] "file2013-03-02_16-27-15.csv"

■対話的入力の例
赤文字はキー入力

例1
> data = readline("input=> ")
input=> 12
> data
[1] "12"

Input_2
カーソルの箇所でキー入力してEnterを押すとdataに文字列として代入される.

例2
> data=readline("input=> ")
input=> kisarazu3
> data
[1] "kisarazu3"

例3 ディレクトリの入力例
> workDIR=readline("Enter Dir=> ")
Enter Dir=> C:/test_tmp
> workDIR
[1] "C:/test_tmp"
> setwd(workDIR)
> getwd()
[1] "C:/test_tmp"

■作業ディレクトリ内にある特定の拡張子をもつファイルの処理例(複数ファイルの処理)

infiles=list.files()   #現在の作業ディレクトリのファイル一覧を代入
file_type="\\.jpeg$|\\.jpg$\\|.JPEG$|\\.JPG$"
for (file.name in infiles) { #ファイルごとの繰り返し
  if (regexpr(file_type, file.name) < 0){
    print("next")
    next
  }
  print("ファイルの内容を読み込む")
  print("読み込んだデータを使った処理")
}


"\\.jpeg$|\\.jpg$\\|.JPEG$|\\.JPG$" の説明は参考サイト(1)
拡張子が.jpeg .jpg .JPEG .JPG のファイルを処理対象とする.

実行例
[1] "next"
[1] "ファイルの内容を読み込む"
[1] "読み込んだデータを使った処理"
> infiles
[1] "test.txt"  "test.jpeg"

■「対話的(interactively)なファイル選択」による
ファイル読み込み(テキストファイルの読み込み)
> read.table(file=file.choose())

■インターネットからのファイル読み込み(テキストファイルの読み込み)
【記述例】
> site="http://www.****"    ここに具体的なurlを記述する
> read.table(url(site),header=TRUE)   または  >  read.table(site,header=TRUE)

■対話的(interactively)なファイル選択(複数選択が可能のようです)
ヘルプで利用方法を調べて使いましょう.
選択結果を変数に代入して各種処理に利用します.
> choose.files()

■対話的(interactively)なディレクトリ選択(一つのディレクトリを選択)
ヘルプで利用方法を調べて使いましょう.
選択結果を変数に代入して各種処理に利用します.
> choose.dir()

■パスを除いたファイル名を取得
・最後のパス分離記号(もし存在すれば)までの全てのパスを削除する.
> basename()
ファイル操作 を参考にさせて頂きました.

■ファイルのパスを取得
・最後のパス分離記号まで(しかし最後は除く)の path の部分を返す。 もしパス分離記号が無ければ "." を返す。チルダ記号の展開が行われる.
>dirname()
ファイル操作 を参考にさせて頂きました.
http://www.is.titech.ac.jp/~mase/mase/html.jp/temp/files.jp.html

■文字列を1文字ごとに分離する
> mojiretu="1234567890123456789012345678"
> bunri=strsplit(mojiretu, split=character(0))  #strsplit(mojiretu, split="")
> bunri
[[1]]
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "0" "1" "2" "3" "4" "5" "6" "7" "8"
[19] "9" "0" "1" "2" "3" "4" "5" "6" "7" "8"
> bunri2=unlist(bunri)
> bunri2
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "0" "1" "2" "3" "4" "5" "6" "7" "8"
[19] "9" "0" "1" "2" "3" "4" "5" "6" "7" "8"

■二進数の文字列を十進整数に変換(入力ミスの処理は無し)
base2=readline(prompt="二進数を入力して下さい ")
ketasu=nchar(base2)
bunri2=unlist(strsplit(base2,split=character(0)))
bunri3=rev(bunri2)
bunri3=as.numeric(bunri3)
keisu=c()
for (i in (0:(ketasu-1))){
   tmp= 2^i
   keisu=c(keisu,tmp)
}
dec=sum(keisu*bunri3)
print(dec)

・使用例
28桁の場合
> source("C:\\R2013\\nishin2dec.r")
二進数を入力して下さい 1111111111111111111111111111
[1] 268435455

■二進数の文字列を十進整数に変換(上記と同じ)
base2="0101000100000010001001000010" #28桁のデータ例
ketasu=nchar(base2)
bunri2=unlist(strsplit(base2, split=character(0)))
bunri3=rev(bunri2)
bunri3=as.numeric(bunri3)
keisu=c()
for (i in (0:(ketasu-1))){
   tmp= 2^i
   keisu=c(keisu,tmp)
}
dec=sum(keisu*bunri3)
print(dec)

■十進数を二進数に変換
dec2baseはパッケージ ‘oro.dicom’ を使うようです.
> dec2base(268435455,2)
[1] "1111111111111111111111111111"     (28桁)

■文字列の結合例
所定のn桁の二進数の左に,0を(28-n)個結合して28桁の二進数の文字列を作る.
nは28以下.

ketasu=28
nishin="10011" #データ例
n=nchar(nishin)
Tmp=rep("0",ketasu-n)
n0=paste(Tmp,collapse="")
nishin2=paste(n0,nishin,sep="")
print(nishin2)

結果
[1] "0000000000000000000000010011"

■参考サイト
(1) データファイルの順次処理:ファイル名のいろいろな決め方
http://takenaka-akio.org/doc/r_auto/chapter_11.html
(2) R における正規表現
http://www.okada.jp.org/RWiki/?R%20%A4%CB%A4%AA%A4%B1%A4%EB%C0%B5%B5%AC%C9%BD%B8%BD
(3) RプログラミングTips大全
http://www.okada.jp.org/RWiki/?R%A5%D7%A5%ED%A5%B0%A5%E9%A5%DF%A5%F3%A5%B0Tips%C2%E7%C1%B4#f3455df4
(4) grep {base}    R Documentation Pattern Matching and Replacement
http://127.0.0.1:11012/library/base/html/grep.html
(5) regex {base}    R Documentation Regular Expressions as used in R
http://stat.ethz.ch/R-manual/R-patched/library/base/html/regex.html
(6) Index of /R-manual/R-patched/library
http://stat.ethz.ch/R-manual/R-patched/library/
(7) The R Manuals
http://cran.r-project.org/manuals.html


2013年1月29日 (火)

■『Rによる空間データの統計分析』で空間データの可視化手法を学ぶ(その3)

■『R による空間データの統計分析』で空間データの可視化手法を学習する(3章の一部)

●書籍
『R による空間データの統計分析』
古谷知之 (著)
朝倉書店 (2011/6/10)
※具体的にコードを入力しながら学習することができると思います.

●参照サイト
空間データの統計分析
サイト内のデータを利用させて頂きました.ここに記して感謝の意を表します.
http://web.sfc.keio.ac.jp/~maunz/wiki/index.php?%B6%F5%B4%D6%A5%C7%A1%BC%A5%BF%A4%CE%C5%FD%B7%D7%CA%AC%C0%CF

※Rのバージョンは2.15.1で行いました.

library(sp)
library(spatstat)
# 正規分布
px <- rnorm(500, mean=0.5, sd=0.15)
py <- rnorm(500, mean=0.5, sd=0.15)
pz <- as.data.frame(rep(1, 500))
colnames(pz) <- c("pz")
plot(pnt, type="p")

Zu_3_1_a_w

plot(density(pnt), 0.1)

Zu_3_1_a_density_w

contour(density(pnt), add=TRUE)
NULL

Zu_3_1_a_density_contour_w

plot(pnt, type="p", add=TRUE)

Zu_3_1_a_density_contour_pnt_w

pnt_xy <- cbind(px, py)
head(pnt_xy)
            px        py
[1,] 0.6703488 0.6826212
[2,] 0.3866471 0.6146397
[3,] 0.5537766 0.2806091
[4,] 0.4628794 0.5317109
[5,] 0.6450791 0.6735637
[6,] 0.6081204 0.5018943

pnt_sp <- SpatialPoints(data.frame(px, py))
head(pnt_sp)
SpatialPoints:
            px        py
[1,] 0.6703488 0.6826212
[2,] 0.3866471 0.6146397
[3,] 0.5537766 0.2806091
[4,] 0.4628794 0.5317109
[5,] 0.6450791 0.6735637
[6,] 0.6081204 0.5018943
Coordinate Reference System (CRS) arguments: NA

class(pnt_sp)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"

pnt_spdf <- SpatialPointsDataFrame(pnt_sp, pz)

class(pnt_spdf)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

pnt_spdf[1:6,]
           coordinates pz
1 (0.670349, 0.682621)  1
2  (0.386647, 0.61464)  1
3 (0.553777, 0.280609)  1
4 (0.462879, 0.531711)  1
5 (0.645079, 0.673564)  1
6  (0.60812, 0.501894)  1

plot(pnt_spdf)

Zu_pnt_spdf_w

quadratcount(pnt, nx=2, ny=2)
         x
y         [0,0.5] (0.5,1]
  (0.5,1]     113     139
  [0,0.5]     126     122
quadratcount(pnt, nx=2, ny=2)/(0.5^2)
         x
y         [0,0.5] (0.5,1]
  (0.5,1]     452     556
  [0,0.5]     504     488

plot(pnt, type="p")
plot(quadratcount(pnt, nx=2, ny=2), add=TRUE, col="red")
Quadratcount_pnt_w

quadratcount(pnt, nx=4, ny=4)
            x
y            [0,0.25] (0.25,0.5] (0.5,0.75] (0.75,1]
  (0.75,1]          1         14         12        0
  (0.5,0.75]        9         89        113       14
  (0.25,0.5]       10        107        106        4
  [0,0.25]          1          8         12        0
quadratcount(pnt, nx=4, ny=4)/(0.25^2)
            x
y            [0,0.25] (0.25,0.5] (0.5,0.75] (0.75,1]
  (0.75,1]         16        224        192        0
  (0.5,0.75]      144       1424       1808      224
  (0.25,0.5]      160       1712       1696       64
  [0,0.25]         16        128        192        0

plot(pnt, type="p")
plot(quadratcount(pnt, nx=4, ny=4), add=TRUE, col="red")

Quadratcount4_4_pnt_w

poly1 <- cbind(
  c(0, 0.75, 0.75, 0.5, 0.5, 0.25, 0.25, 0, 0),
c(0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 0))
poly2 <- cbind(
  c(0.25, 0.75, 0.75, 0.25, 0.25),
  c(0.5, 0.5, 0.75, 0.75, 0.5))
poly3 <- cbind(
  c(0.75, 1, 1, 0.75, 0.75, 0.5, 0.5, 0.75, 0.75),
  c(0, 0, 0.75, 0.75, 0.5, 0.5, 0.25, 0.25, 0))
poly4 <- cbind(
  c(0, 1, 1, 0, 0),
  c(0.75, 0.75, 1, 1, 0.75))
poly1
      [,1] [,2]
[1,] 0.00 0.00
[2,] 0.75 0.00
[3,] 0.75 0.25
[4,] 0.50 0.25
[5,] 0.50 0.50
[6,] 0.25 0.50
[7,] 0.25 0.75
[8,] 0.00 0.75
[9,] 0.00 0.00
poly2
     [,1] [,2]
[1,] 0.25 0.50
[2,] 0.75 0.50
[3,] 0.75 0.75
[4,] 0.25 0.75
[5,] 0.25 0.50
poly3
      [,1] [,2]
[1,] 0.75 0.00
[2,] 1.00 0.00
[3,] 1.00 0.75
[4,] 0.75 0.75
[5,] 0.75 0.50
[6,] 0.50 0.50
[7,] 0.50 0.25
[8,] 0.75 0.25
[9,] 0.75 0.00
poly4
     [,1] [,2]
[1,]    0 0.75
[2,]    1 0.75
[3,]    1 1.00
[4,]    0 1.00
[5,]    0 0.75

poly1_pl <- Polygons(list(Polygon(poly1)), "poly1")

poly1_pl
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 0.2916667 0.2916667

Slot "area":
[1] 0.375

Slot "hole":
[1] FALSE

Slot "ringDir":
[1] 1

Slot "coords":
      [,1] [,2]
[1,] 0.00 0.00
[2,] 0.00 0.75
[3,] 0.25 0.75
[4,] 0.25 0.50
[5,] 0.50 0.50
[6,] 0.50 0.25
[7,] 0.75 0.25
[8,] 0.75 0.00
[9,] 0.00 0.00

Slot "plotOrder":
[1] 1

Slot "labpt":
[1] 0.2916667 0.2916667

Slot "ID":
[1] "poly1"

Slot "area":
[1] 0.375


summary(poly1_pl)
  Length    Class     Mode
       1 Polygons       S4


poly2_pl <- Polygons(list(Polygon(poly2)), "poly2")
poly3_pl <- Polygons(list(Polygon(poly3)), "poly3")
poly4_pl <- Polygons(list(Polygon(poly4)), "poly4")

poly_sp <- SpatialPolygons(list(poly1_pl, poly2_pl, poly3_pl, poly4_pl))
poly_spdf <- SpatialPolygonsDataFrame(poly_sp, data.frame(c(1:4),
row.names=c("poly1", "poly2", "poly3", "poly4")))

class(poly_sp)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"

poly_spdf <- SpatialPolygonsDataFrame(poly_sp, data.frame(c(1:4),
  row.names=c("poly1", "poly2", "poly3", "poly4")))

class(poly_spdf)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

plot(poly_spdf)
plot(pnt, type="p", add=TRUE)

Zu_3_1_d_w

overlay(pnt_spdf, poly_spdf, fn=sum)
   pz
1 147
2 202
3 124
4  27

ov1=overlay(pnt_spdf, poly_spdf, fn=sum)
ov1
   pz
1 147
2 202
3 124
4  27

P1= poly_spdf@polygons[[1]]@labpt
P2= poly_spdf@polygons[[2]]@labpt
P3= poly_spdf@polygons[[3]]@labpt
P4= poly_spdf@polygons[[4]]@labpt
text(P1[1],P1[2],ov1[1],col="red",cex=1.4)
text(P2[1],P2[2],ov1[2],col="red",cex=1.4)
text(P3[1],P3[2],ov1[3],col="red",cex=1.4)
text(P4[1],P4[2],ov1[4],col="red",cex=1.4)

Zu_3_1_d2_w

area.poly_spdf <- c(0.25^2*6, 0.25^2*2, 0.25^2*4, 0.25^2*4)
area.poly_spdf
[1] 0.375 0.125 0.250 0.250

ov2=overlay(pnt_spdf, poly_spdf, fn=sum)/area.poly_spdf
ov2
    pz
1  392
2 1616
3  496
4  108

col=gray.colors(4, start = 0.6, end = 0.8)
plot(poly_spdf,col=col)
plot(pnt, type="p", col="blue",cex=1.5,add=TRUE)

P1= poly_spdf@polygons[[1]]@labpt
P2= poly_spdf@polygons[[2]]@labpt
P3= poly_spdf@polygons[[3]]@labpt
P4= poly_spdf@polygons[[4]]@labpt
text(P1[1],P1[2],ov2[1],col="red",cex=1.4)
text(P2[1],P2[2],ov2[2],col="red",cex=1.4)
text(P3[1],P3[2],ov2[3],col="red",cex=1.4)
text(P4[1],P4[2],ov2[4],col="red",cex=1.4)

Zu_3_2_w

■参考
Polygon(poly1)
An object of class "Polygon"
Slot "labpt":
[1] 0.2916667 0.2916667

Slot "area":
[1] 0.375

Slot "hole":
[1] TRUE

Slot "ringDir":
[1] -1

Slot "coords":
      [,1] [,2]
[1,] 0.00 0.00
[2,] 0.75 0.00
[3,] 0.75 0.25
[4,] 0.50 0.25
[5,] 0.50 0.50
[6,] 0.25 0.50
[7,] 0.25 0.75
[8,] 0.00 0.75
[9,] 0.00 0.00

list(Polygon(poly1))
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 0.2916667 0.2916667

Slot "area":
[1] 0.375

Slot "hole":
[1] TRUE

Slot "ringDir":
[1] -1

Slot "coords":
      [,1] [,2]
[1,] 0.00 0.00
[2,] 0.75 0.00
[3,] 0.75 0.25
[4,] 0.50 0.25
[5,] 0.50 0.50
[6,] 0.25 0.50
[7,] 0.25 0.75
[8,] 0.00 0.75
[9,] 0.00 0.00


■参考サイト
・Package ‘sp’December 21, 2012
http://cran.r-project.org/web/packages/sp/sp.pdf

Package ‘spatstat’January 17, 2013
http://cran.r-project.org/web/packages/spatstat/spatstat.pdf

2013年1月14日 (月)

■『Rによる空間データの統計分析』で空間データの可視化手法を学ぶ(その2)

■『R による空間データの統計分析』で空間データの可視化手法を学習する(4章の一部)

●書籍
『R による空間データの統計分析』
古谷知之 (著)
朝倉書店 (2011/6/10)
※具体的にコードを入力しながら学習することができると思います.

●参照サイト
空間データの統計分析
サイト内のデータを利用させて頂きました.ここに記して感謝の意を表します.
http://web.sfc.keio.ac.jp/~maunz/wiki/index.php?%B6%F5%B4%D6%A5%C7%A1%BC%A5%BF%A4%CE%C5%FD%B7%D7%CA%AC%C0%CF

※Rのバージョンは2.15.1で行いました.

4.1 等量分類  p.44

library(maptools)
library(classInt)

jpn_pref <- readShapePoly("C:/rkukan_work/jpn_pref.shp",IDvar="PREF_CODE")
names(jpn_pref)
[1] "PREF"      "PREF_CODE" "Area"      "Habitable" "Pop"       "Pop_Dens"
summary(jpn_pref)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min       max
x 122.9339 153.98243
y  24.0448  45.55691
Is projected: NA
proj4string : [NA]
Data attributes:
      PREF      PREF_CODE         Area         Habitable            Pop              Pop_Dens    
Aichi  : 1   Min.   : 1.0   Min.   : 1876   Min.   :  850.5   Min.   :  607012   Min.   : 257.0 
Akita  : 1   1st Qu.:12.5   1st Qu.: 4165   1st Qu.: 1302.2   1st Qu.: 1163534   1st Qu.: 654.9 
Aomori : 1   Median :24.0   Median : 6096   Median : 2022.0   Median : 1753179   Median : 852.6 
Chiba  : 1   Mean   :24.0   Mean   : 8039   Mean   : 2582.6   Mean   : 2718468   Mean   :1372.1 
Ehime  : 1   3rd Qu.:35.5   3rd Qu.: 8087   3rd Qu.: 2952.9   3rd Qu.: 2762151   3rd Qu.:1223.3 
Fukui  : 1   Max.   :47.0   Max.   :83456   Max.   :21900.7   Max.   :12576601   Max.   :9009.5 
(Other):41 

jpn_COD <- read.table("C:/rkukan_work/COD.csv",sep=",",header=TRUE)
names(jpn_COD)
[1] "PREF_CODE"    "geriatric"    "malignant"    "diabetes"     "hypertensive"
summary(jpn_COD)
   PREF_CODE      geriatric       malignant        diabetes      hypertensive 
Min.   : 1.0   Min.   :351.9   Min.   :185.8   Min.   : 7.80   Min.   :2.400 
1st Qu.:12.5   1st Qu.:514.6   1st Qu.:255.9   1st Qu.:10.40   1st Qu.:4.000 
Median :24.0   Median :561.4   Median :271.8   Median :11.50   Median :4.800 
Mean   :24.0   Mean   :560.7   Mean   :273.9   Mean   :11.51   Mean   :4.943 
3rd Qu.:35.5   3rd Qu.:619.4   3rd Qu.:293.6   3rd Qu.:12.35   3rd Qu.:5.700 
Max.   :47.0   Max.   :688.1   Max.   :337.9   Max.   :18.00   Max.   :8.700 

# 地図属性テーブルと空間データテーブルとの結合
ID.match <- match(jpn_pref$PREF_CODE, jpn_COD$PREF_CODE)
ID.match
[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
[40] 40 41 42 43 44 45 46 47

jpn_COD1 <- jpn_COD[ID.match,]
names(jpn_COD1)
[1] "PREF_CODE"    "geriatric"    "malignant"    "diabetes"     "hypertensive"
summary(jpn_COD1)
   PREF_CODE      geriatric       malignant        diabetes      hypertensive 
Min.   : 1.0   Min.   :351.9   Min.   :185.8   Min.   : 7.80   Min.   :2.400 
1st Qu.:12.5   1st Qu.:514.6   1st Qu.:255.9   1st Qu.:10.40   1st Qu.:4.000 
Median :24.0   Median :561.4   Median :271.8   Median :11.50   Median :4.800 
Mean   :24.0   Mean   :560.7   Mean   :273.9   Mean   :11.51   Mean   :4.943 
3rd Qu.:35.5   3rd Qu.:619.4   3rd Qu.:293.6   3rd Qu.:12.35   3rd Qu.:5.700 
Max.   :47.0   Max.   :688.1   Max.   :337.9   Max.   :18.00   Max.   :8.700 

jpn_pref_COD <- spCbind(jpn_pref,jpn_COD1)
names(jpn_pref_COD)
[1] "PREF"         "PREF_CODE"    "Area"         "Habitable"    "Pop"          "Pop_Dens"     "PREF_CODE.1"
[8] "geriatric"    "malignant"    "diabetes"     "hypertensive"

summary(jpn_pref_COD)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min       max
x 122.9339 153.98243
y  24.0448  45.55691
Is projected: NA
proj4string : [NA]
Data attributes:
      PREF      PREF_CODE         Area         Habitable            Pop              Pop_Dens       PREF_CODE.1 
Aichi  : 1   Min.   : 1.0   Min.   : 1876   Min.   :  850.5   Min.   :  607012   Min.   : 257.0   Min.   : 1.0 
Akita  : 1   1st Qu.:12.5   1st Qu.: 4165   1st Qu.: 1302.2   1st Qu.: 1163534   1st Qu.: 654.9   1st Qu.:12.5 
Aomori : 1   Median :24.0   Median : 6096   Median : 2022.0   Median : 1753179   Median : 852.6   Median :24.0 
Chiba  : 1   Mean   :24.0   Mean   : 8039   Mean   : 2582.6   Mean   : 2718468   Mean   :1372.1   Mean   :24.0 
Ehime  : 1   3rd Qu.:35.5   3rd Qu.: 8087   3rd Qu.: 2952.9   3rd Qu.: 2762151   3rd Qu.:1223.3   3rd Qu.:35.5 
Fukui  : 1   Max.   :47.0   Max.   :83456   Max.   :21900.7   Max.   :12576601   Max.   :9009.5   Max.   :47.0 
(Other):41                                                                                                      
   geriatric       malignant        diabetes      hypertensive 
Min.   :351.9   Min.   :185.8   Min.   : 7.80   Min.   :2.400 
1st Qu.:514.6   1st Qu.:255.9   1st Qu.:10.40   1st Qu.:4.000 
Median :561.4   Median :271.8   Median :11.50   Median :4.800 
Mean   :560.7   Mean   :273.9   Mean   :11.51   Mean   :4.943 
3rd Qu.:619.4   3rd Qu.:293.6   3rd Qu.:12.35   3rd Qu.:5.700 
Max.   :688.1   Max.   :337.9   Max.   :18.00   Max.   :8.700


# カラーパレットの作成
#pal0 <- c("white","grey","grey2")
#pal0 <- c("white",grey(1),grey(0.5))
pal0 <- c("antiquewhite","darkorange1","darkorange3")

# 4.1 等量分類
q_pref <- classIntervals(round(jpn_pref_COD$malignant,2), n=5, style="quantile")
q_pref
style: quantile
  one of 163,185 possible partitions of this variable into 5 classes
[185.8,251.68) [251.68,266.42)  [266.42,278.5)   [278.5,304.6)
             10               9               9               9
  [304.6,337.9]
             10


plot(q_pref, pal=pal0, cex.axis=1.3, cex.lab=1.2, lwd=2, main="",col="blue",
xlab="悪性新生物による死亡者数(人口10万人あたり)")

Zu_4_1a

q_pref_Col <- findColours(q_pref,pal0)
plot(jpn_pref_COD,col=q_pref_Col)

Zu_4_2a

title("悪性新生物による死亡者数(人口10万人あたり) (等量分類)", cex=1.4)
legend("topleft",fill=attr(q_pref_Col,"palette"), cex=1.4,
legend=names(attr(q_pref_Col,"table")),bty="n")

Zu_4_3a

■Global Administrative AreasからDLした日本地図による作図
下記のサイトから入手した都道府県境界マップを使用した場合のメモ.
Download | Global Administrative Areas
http://www.gadm.org/country
Global Administrative Areas | Boundaries without limits
http://www.gadm.org/

JPN_adm0:日本マップ
JPN_adm1:都道府県境界マップ
JPN_adm2:市町村境界マップ(あらいようです)
kmzファイル
などが入手可能のようです.

今回は,「JPN_adm1:都道府県境界マップ」を利用させて頂きました.

都道府県別コード(ID_1)はアルファベット順位となっているようです.

QGISのTable ManagerによるJPN_adm1の属性表示
Jpn_adm1_zokusei

そこで,『COD.csv』(「Rによる空間データの統計分析」のサイトから入手可)ファイル
にアルファベット順のコード(PREF_CODE2)を1列目に追加し,別名保存し利用(下図).
Zu_cod_b

・作図結果
Zu_4_adm1_a

Zu_4_adm2_a

※島の一部が表示されない.

以下は,この
Global Administrative AreasからDLした日本地図を用いた.
(都道府県名はアルファベットとなっている)

Global Administrative Areasからのマップを使っているため
分かりにくいと思います.すみません.

■前作業
・以下の作業のために,都道府県名をアルファベット順に整理した
CODe_1.csvを作成

Code_1

4.8 ドットマップ   p.52

# 人口5万人に1ドット
#tohoku_dots <- dotsInPolys(tohoku_COD, as.integer(tohoku_COD$Pop/50000))
エラー:  関数 "sample.Polygons" を見つけることができませんでした

となり作図できなかった.

4.9 シンボルマップ   p.53

jpn_pref_b <- readShapePoly("C:/rkukan_work/PN_adm1.shp",IDvar="ID_1")
jpn_COD_e <- read.table("C:/rkukan_work/CODe_1.csv",sep=",",header=TRUE)
IIDe.match <- match(jpn_pref_b$ID_1, jpn_COD_e$PREF_CODE2)
IDe.match
[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
[27] 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

jpn_COD1_e  <- jpn_COD_e[IDe.match,]
jpn_pref_b_COD_e <- spCbind(jpn_pref_b,jpn_COD1_e)
names(jpn_pref_b_COD_e)
[1] "ID_0"         "ISO"          "NAME_0"       "ID_1"         "NAME_1"       "NL_NAME_1"   
[7] "VARNAME_1"    "TYPE_1"       "ENGTYPE_1"    "PREF_CODE2"   "PREF_CODE"    "geriatric"   
[13] "malignant"    "diabetes"     "hypertensive" "PREF_name"    "Area"         "Habitable"   
[19] "Pop"          "Pop_Dens"   

# 東北地方を切り出す
tohoku_COD <- jpn_pref_b_COD_e[jpn_pref_b_COD_e$PREF_CODE>=2,]
tohoku_COD <- tohoku_COD[tohoku_COD$PREF_CODE<=7,]
#names(tohoku_COD)
#summary(tohoku_COD)
# 座標の抽出
tohoku_coord <- coordinates(tohoku_COD)
#names(tohoku_coord)
#summary(tohoku_coord)
# 図4.9  p.53
# ポリゴンデータを表示
plot(tohoku_COD, col="grey", border="white", lwd=2)
symbols(x=tohoku_coord[,1], y=tohoku_coord[,2],
circles=tohoku_COD$hypertensive/35, inch=FALSE, bg="black", add=TRUE)
text(x=tohoku_coord[,1], y=tohoku_coord[,2]+0.3, cex=1.3, col="black",
c("秋田県","青森県","福島県","岩手県", "宮城県","山形県"))

#県名をアルファベット順に設定していることにより,
#"秋田県","青森県","福島県","岩手県", "宮城県","山形県"の順にするようです.

Zu_4_9_c_1

4.9 複数の属性データの表示   p.54

# 図4.10  p.55
# パッケージTeachingDemosを使用
library(TeachingDemos)
# 棒グラフで表示するデータテーブルを作成
tohoku_COD2 <- cbind(tohoku_COD$Pop_Dens,tohoku_COD$malignant)
# ポリゴンデータを表示
plot(tohoku_COD, lwd=2,col="antiquewhite",border="ivory4")
# 棒グラフを表示
for(i in 1:nrow(tohoku_COD)){
subplot(barplot(tohoku_COD2[i,], yaxt="n", col=c("cyan4","indianred")),
x=tohoku_coord[i,1],y=tohoku_coord[i,2], vadj=0, size=c(0.4,0.6))
}
# 凡例を表示
legend(138, 41.5, c("人口密度", paste("悪性新生物","死亡者数",sep="\n")),
cex=1.3, fill=c("cyan4","indianred"), bty="n")

Zu_4_10_1

■色に関するサイトなどの参考URL
統計グラフの色
http://oku.edu.mie-u.ac.jp/~okumura/stat/colors.html
Colorbrewer: Color Advice for Maps
http://colorbrewer2.org/
Colors in R
http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
Download | Global Administrative Areas
http://www.gadm.org/country
Global Administrative Areas | Boundaries without limits
http://www.gadm.org/

■世界のCO2排出量マップの作図
●データ入手先
CO2 emissions from fuel combustion
http://www.oecd-ilibrary.org/environment/co2-emissions-from-fuel-combustion_2075826x-table1

●作図例
【CO2 emissions from fuel combustion(2009) [Million tonnes]】
データ:40カ国
Australia Austria Belgium Brazil Canada Switzerland Chile
China Czech Republic Germany Denmark Spain Estonia
Finland France United Kingdom Greece Hungary Indonesia
India Ireland Iceland Israel Italy Japan South_Korea
Luxembourg Mexico Netherlands Norway New Zealand Poland
Portugal Russia Slovakia Slovenia Sweden Turkey
United States South_Africa

Co2_fn_b
 
Co2_world_map2009_b1

Fix_map_co2_b

Co2_world_map2009_iso_b

●上図に対するコード(参考,試行錯誤で描いた図)
library(maptools)
library(classInt)
#world <- readShapePoly("C:/R2013/rkukan_work/co2/all_co2/all_co2.shp",IDvar="ISO")
world <- readShapePoly("C:/R2013/rkukan_work/chizu_world.shp",IDvar="ISO")
world_CO2 <- read.table("C:/R2013/rkukan_work/co2/CO2_data_2.csv",sep=",",header=TRUE)
ID.match <- match(world$ISO, world_CO2$ISO)
world_CO2_1 <- world_CO2[ID.match,]
row.names(world_CO2_1)<- world$ISO
world_county_CO2<- spCbind(world,world_CO2_1)
pal <- cm.colors(n=10,alpha=1)
fix_pref <- classIntervals(world_county_CO2$X2009, n=10, style="fixed",
fixedBreaks=c(2, 10,50, 100, 500, 1000,2000,3000,4000,5000,6832))

plot(fix_pref, pal=pal, cex.axis=1.3, cex.lab=1.2, lwd=2, main="",col="blue",
xlab="CO2 emissions from fuel combustion(2009) [Million tonnes]")

fix_pref_Col <- findColours(fix_pref,pal)
plot(world_county_CO2,col=fix_pref_Col,border="gray90",cex=0.05)
警告メッセージ:
1: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)
2: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)
3: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)
4: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)
5: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)
6: In polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  :
  Reached total allocation of 1535Mb: see help(memory.size)

world_coord<- coordinates(world_county_CO2)
world_coord
           [,1]        [,2]
AUS  134.355123 -25.5767961
AUT   14.144274  47.5881231
BEL    4.661211  50.6409251
BRA  -53.080320 -10.7425745
CAN -102.336873  57.9433211
CHE    8.241408  46.8104133
CHL  -71.231826 -35.3089075
CHN  103.832182  36.6197849
CZE   15.332121  49.7391324
DEU   10.386771  51.0967695
DNK    9.314695  56.0332699
ESP   -3.554164  40.3905257
EST   25.843485  58.6830018
FIN   26.287923  64.5199965
FRA    2.459853  46.6250364
GBR   -2.522990  53.9440593
GRC   22.576521  39.4722305
HUN   19.411116  47.1676406
IDN  114.001523  -0.2046302
IND   79.583090  22.9112341
IRL   -8.146708  53.1799074
ISL  -18.594070  64.9906731
ISR   35.023004  31.4742777
ITA   12.162522  43.5248900
JPN  137.977928  36.6522662
KOR  127.849389  36.4545149
LUX    6.093842  49.7748812
MEX -102.527156  23.9438089
NLD    5.621097  52.2725774
NOR   13.991355  64.2386706
NZL  170.480865 -43.9854127
POL   19.408963  52.1215081
PRT   -7.963010  39.6861197
RUS   99.132607  61.6644109
SVK   19.483510  48.7094553
SVN   14.818129  46.1244741
SWE   16.743484  62.8337959
TUR   35.436876  38.9880081
USA  -99.152262  39.5294320
ZAF   25.166140 -29.0107588

legend("top",fill=attr(fix_pref_Col,"palette"), cex=0.8,
+ legend=names(attr(fix_pref_Col,"table")),bty="n")
title("CO2 emissions from fuel combustion(2009) [Million tonnes]", cex=1)
cp<- rep(1,40)
symbols(x=world_coord[,1],y=world_coord[,2],circles=cp/2,inch=FALSE,add=TRUE,bg="gray70")
text(x=world_coord[,1],y=world_coord[,2],world$ISO,cex=0.4,col="gray82")

メモ
1) chizu_world.shpの作成
各国のshapeファイルをDownload | Global Administrative AreasからDLし,
DLした40個のshapeファイル(JPN_adm0.shp など)を結合する.

結合には,『GeoMerge』を使用させて頂きました.
GeoMerge - VDS TECHNOLOGIES
http://www.vdstech.com/geomerge.aspx

作成されたdbfファイルのイメージ
Shape_table

2) CO2_data_2.csvの作成
(ファイル名をCO2_sample.csvとしたものを参考として置きます).
「CO2_sample.zip」 (2KB)

※1)と2)のファイルのレコードの国別順番は同じにする.
具体的には,1)に合わせた.

3) plot(world_county_CO2,col=fix_pref_Col,border="gray90",cex=0.05)
でメモリに関する警告が表示された・・・.対処法分かりません.

4) 参考パッケージなど
Package ‘maptools’January 16, 2013
http://cran.r-project.org/web/packages/maptools/maptools.pdf
Package ‘classInt’November 5, 2012
http://cran.r-project.org/web/packages/classInt/classInt.pdf
Introduction to spatstat
http://www.stats.uwo.ca/faculty/kulperger/S9934a/Papers/SpatStatIntro.pdf

統計局ホームページ
http://www.stat.go.jp/index.htm
統計局ホームページ/世界の統計
http://www.stat.go.jp/data/sekai/index.htm
統計局ホームページ/世界の統計 第16章 環境
http://www.stat.go.jp/data/sekai/16.htm

●GeoMergeのイメージ
Geomerge_2

●Rの履歴の補足
world_coord<- coordinates(world_county_CO2)
world_coord
           [,1]        [,2]
AUS  134.355123 -25.5767961
AUT   14.144274  47.5881231
BEL    4.661211  50.6409251
BRA  -53.080320 -10.7425745
CAN -102.336873  57.9433211
CHE    8.241408  46.8104133
CHL  -71.231826 -35.3089075
CHN  103.832182  36.6197849
CZE   15.332121  49.7391324
DEU   10.386771  51.0967695
DNK    9.314695  56.0332699
ESP   -3.554164  40.3905257
EST   25.843485  58.6830018
FIN   26.287923  64.5199965
FRA    2.459853  46.6250364
GBR   -2.522990  53.9440593
GRC   22.576521  39.4722305
HUN   19.411116  47.1676406
IDN  114.001523  -0.2046302
IND   79.583090  22.9112341
IRL   -8.146708  53.1799074
ISL  -18.594070  64.9906731
ISR   35.023004  31.4742777
ITA   12.162522  43.5248900
JPN  137.977928  36.6522662
KOR  127.849389  36.4545149
LUX    6.093842  49.7748812
MEX -102.527156  23.9438089
NLD    5.621097  52.2725774
NOR   13.991355  64.2386706
NZL  170.480865 -43.9854127
POL   19.408963  52.1215081
PRT   -7.963010  39.6861197
RUS   99.132607  61.6644109
SVK   19.483510  48.7094553
SVN   14.818129  46.1244741
SWE   16.743484  62.8337959
TUR   35.436876  38.9880081
USA  -99.152262  39.5294320
ZAF   25.166140 -29.0107588

names(world_county_CO2)
[1] "ID_0"         "ISO"          "NAME_ENGLI"   "NAME_ISO"   
[5] "NAME_FAO"     "NAME_LOCAL"   "NAME_OBSOL"   "NAME_VARIA" 
[9] "NAME_NONLA"   "NAME_FRENC"   "NAME_SPANI"   "NAME_RUSSI" 
[13] "NAME_ARABI"   "NAME_CHINE"   "WASPARTOF"    "CONTAINS"   
[17] "SOVEREIGN"    "ISO2"         "WWW"          "FIPS"       
[21] "ISON"         "VALIDFR"      "VALIDTO"      "EUmember"   
[25] "ID_0.1"       "ISO.1"        "NAME_ENGLI.1" "X2002"      
[29] "X2003"        "X2004"        "X2005"        "X2006"      
[33] "X2007"        "X2008"        "X2009"      



2013年1月 5日 (土)

『Rによる空間データの統計分析』で空間データの可視化手法を学ぶ(その1)

■『Rによる空間データの統計分析』で空間データの可視化手法を学習する(2章の一部)

●書籍
『Rによる空間データの統計分析』
古谷知之 (著)
朝倉書店 (2011/6/10)

※具体的にコードを入力しながら学習することができると思います.

●サイト
空間データの統計分析
サイト内のデータを利用させて頂きました.ここに記して感謝の意を表します.
http://web.sfc.keio.ac.jp/~maunz/wiki/index.php?%B6%F5%B4%D6%A5%C7%A1%BC%A5%BF%A4%CE%C5%FD%B7%D7%CA%AC%C0%CF

●利用するRのバージョン:version2.12.2 または 2.15,1
R_version

R_ver2_15_1

●関連サイト
GISA2012HandsOnSpatStat
http://web.sfc.keio.ac.jp/~maunz/wiki/index.php?GISA2012HandsOnSpatStat

■学習履歴メモ
●【Rの分析例 p17】
二次元平面上の座標にプロットされた5つのポイントデータについての
直線距離とマンハッタン距離

# ポイントデータの座標
x <- c(1, 3, 4, 8, 9)
y <- c(4, 2, 7, 4, 10)
xy <- t(rbind(x,y))

#ベクトルの結合
rbind(x,y)

【出力】
  [,1] [,2] [,3] [,4] [,5]
x    1    3    4    8    9
y    4    2    7    4   10

#転置行列
t(rbind(x,y))

【出力】
     x  y
[1,] 1  4
[2,] 3  2
[3,] 4  7
[4,] 8  4
[5,] 9 10

#作図
plot(xy,xlim=c(0,10),ylim=c(0,10),cex=2, pch=19)

【出力】
Fig_p17_1_2
#垂直線,水平線を引く
abline(h=0:10,v=0:10,lty=2) #lty=2 線プロットのスタイルを破線に指定

【出力】
Fig_p17_2
# ユークリッド距離
dist(xy, method="euclidean")

【出力】
          1         2         3         4
2  2.828427                              
3  4.242641  5.099020                   
4  7.000000  5.385165  5.000000         
5 10.000000 10.000000  5.830952  6.082763

# マンハッタン距離
dist(xy, method="manhattan")

【出力】
   1  2  3  4
2  4         
3  6  6      
4  7  7  7   
5 14 14  8  7

■全国市区町村界データのシェープファイルは
ESRIジャパン株式会社のホームページからDLさせて頂きました.
ここに記して感謝の意を表します.
※『使用規定で許諾された範囲内でのみ』使用できるとなっています.
http://www.esri.com/
http://www.esrij.com/products/data/japan-shp/
(ここでは使用していない)

●【R分析例 p23】
【例1】住んでいる地域を学習対象とするとよいと思います.
・学習に使うshapeファイルを準備する
 ここでは,『地図で見る統計(統計GIS)』から入手した木更津市の境界データを使用.

・入手先:『地図で見る統計(統計GIS)』
国勢調査データ【平成22年国勢調査(小地域)2010/10/01】のダウンロード
サービスからDLした境界データ.
世界測地系緯度経度・Shape形式 木更津市(319KB)

参考サイト
http://kiyomidais.cocolog-nifty.com/blog/2012/12/20101001dl-29e1.html

・参考サイト
国土数値情報 行政区域データ の入手先(ここでは使用していない)
国土数値情報ダウンロードサービス
http://nlftp.mlit.go.jp/ksj/index.html

kisa <- readShapeSpatial("C:/rkukan_work/h22ka12206.shp")
kisa_coord<- coordinates(kisa)
names(kisa)
[1] "AREA"       "PERIMETER"  "H22KA07_"   "H22KA07_ID" "KEN"      
[6] "CITY"       "KEN_NAME"   "SITYO_NAME" "GST_NAME"   "CSS_NAME" 
[11] "HCODE"      "KIHON1"     "DUMMY1"     "KIHON2"     "KEYCODE1" 
[16] "KEYCODE2"   "AREA_MAX_F" "KIGO_D"     "N_KEN"      "N_CITY"   
[21] "N_C1"       "KIGO_E"     "KIGO_I"     "TATE"       "DIR"      
[26] "HIGHT"      "JIKAKU"     "NMOJI"      "MOJI"       "SEQ_NO2"   
[31] "KSUM"       "CSUM"       "JINKO"      "SETAI"      "X_CODE"   
[36] "Y_CODE"     "KCODE1"     "KEY_CODE"   "H22KA08_"   "H22KA08_ID"
[41] "H22KA09_"   "H22KA09_ID" "H22KA10_"   "H22KA10_ID" "H22KA11_" 
[46] "H22KA11_ID" "H22KA12_"   "H22KA12_ID" "H22KA13_"   "H22KA13_ID"
[51] "H22KA14_"   "H22KA14_ID"


plot(kisa, col="grey", border="red")

Kisa_a

text(kisa_coord[,1], kisa_coord[,2], kisa$KEY_CODE, cex=0.7, col="black")

Kisa__b

この図には海域が含まれている.

次に,海域を除く処理を行う.

which(kisa$KEY_CODE=="122069990")  #KEY_CODE=="122069990"のある行を取得
[1] 7
kisa_b<- kisa[-7,]  #7行を除く
kisa_b_coord<- coordinates(kisa_b)
names(kisa_b)
[1] "AREA"       "PERIMETER"  "H22KA07_"   "H22KA07_ID" "KEN"        "CITY"      
[7] "KEN_NAME"   "SITYO_NAME" "GST_NAME"   "CSS_NAME"   "HCODE"      "KIHON1"   
[13] "DUMMY1"     "KIHON2"     "KEYCODE1"   "KEYCODE2"   "AREA_MAX_F" "KIGO_D"   
[19] "N_KEN"      "N_CITY"     "N_C1"       "KIGO_E"     "KIGO_I"     "TATE"      
[25] "DIR"        "HIGHT"      "JIKAKU"     "NMOJI"      "MOJI"       "SEQ_NO2"   
[31] "KSUM"       "CSUM"       "JINKO"      "SETAI"      "X_CODE"     "Y_CODE"   
[37] "KCODE1"     "KEY_CODE"   "H22KA08_"   "H22KA08_ID" "H22KA09_"   "H22KA09_ID"
[43] "H22KA10_"   "H22KA10_ID" "H22KA11_"   "H22KA11_ID" "H22KA12_"   "H22KA12_ID"
[49] "H22KA13_"   "H22KA13_ID" "H22KA14_"   "H22KA14_ID"


plot(kisa_b, col="grey", border="red")
Kisa_01_a 
text(kisa_b_coord[,1], kisa_b_coord[,2], kisa_b$KEY_CODE, cex=0.7, col="black")
Kisa_02_a

左上が「アクアラインと海ほたる」です.

・属性テーブルの名前の確認(確認だけです.編集などは未経験です)
dbfファイルをRで読み,属性名を確認する.
リスクを避けるため,もとのdbfファイルを複製して作業を行った.
【h22ka12206.dbf】を複製し,【copy_h22ka12206.dbf】とした.

library(foreign)
tmp_dbf <- read.dbf("C:/rkukan_work/copy_h22ka12206.dbf")
names(tmp_dbf)
[1] "dbf"    "header"
head(tmp_dbf$dbf,n=3L)
       AREA PERIMETER H22KA07_ H22KA07_ID KEN CITY KEN_NAME SITYO_NAME GST_NAME
1  161879.2 10425.770        0          0  12  206   千葉県       <NA> 木更津市
2 1095848.0  5863.181        0          0  12  206   千葉県       <NA> 木更津市
3 3363032.0 18400.620        0          0  12  206   千葉県       <NA> 木更津市
  CSS_NAME HCODE KIHON1 DUMMY1 KIHON2  KEYCODE1 KEYCODE2 AREA_MAX_F KIGO_D N_KEN
1     <NA>  8101   8150      -     00 206815000  2068150          M   <NA>  <NA>
2     <NA>  8101   6650      -     00 206665000  2066650          M   <NA>  <NA>
3     <NA>  8101   6620      -     00 206662000  2066620          M   <NA>  <NA>
  N_CITY N_C1 KIGO_E KIGO_I TATE DIR HIGHT JIKAKU NMOJI     MOJI SEQ_NO2 KSUM CSUM
1   <NA>    0   <NA>   <NA>    0   0    50     10     4 中島地先    5007    1    1
2   <NA>    0   <NA>   <NA>    0   0    50     10     2     牛込    5115    3    3
3   <NA>    0   <NA>   <NA>    0   0    50     10     2     中島    5136   12   12
  JINKO SETAI   X_CODE   Y_CODE  KCODE1 
KEY_CODE H22KA08_ H22KA08_ID H22KA09_
1     0     0 139.8842 35.45755 8150-00 122068150        0          0        0
2   552   178 139.9408 35.43908 6650-00 122066650        0          0        0
3  2076   599 139.9240 35.43185 6620-00 122066620        0          0        0
  H22KA09_ID H22KA10_ H22KA10_ID H22KA11_ H22KA11_ID H22KA12_ H22KA12_ID H22KA13_
1          0        0          0        0          0     5008       5007        0
2          0        0          0        0          0     5116       5115        0
3          0        0          0        0          0     5137       5136        0
  H22KA13_ID H22KA14_ H22KA14_ID
1          0        0          0
2          0        0          0
3          0        0          0


tmp_dbf$dbf$KEY_CODE
  [1] 122068150   122066650   122066620   122066660   122066630   122066640 
  [7]
122069990   122064550   122064560   12206453000 122064570   122067760 
[13] 12206450004 122064570   12206453003 12206450003 12206451003 12206450001
[19] 12206452003 122067770   12206451002 122067780   12206450002 12206452002
[25] 12206451001 12206453002 122064580   12206452001 122067750   12206453004
[31] 122067720   12206453000 12206453001 122064600   12206104000 122064540 
[37] 122067700   122067740   122067710   12206459000 122067730   12206459002
[43] 122063300   12206340000 122063420   12206459001 122063411   12206104002
[49] 12206104000 12206332000 12206340002 12206340003 12206104001 12206340001
[55] 12206340004 12206332002 12206341201 12206332001 122063460   12206341202
[61] 122063430   12206104000 12206333001 122068080   122061030   122063450 
[67] 122063440   12206350002 12206106000 12206334001 12206334003 122061250 
[73] 12206334002 12206105003 12206105002 12206119001 12206105001 12206333002
[79] 12206102001 12206102002 122061200   12206102003 12206106002 12206101003
[85] 12206106001 12206335001 12206101001 12206106003 12206119002 12206335004
[91] 12206107001 12206107002 12206350003 12206333003 12206101002 122068010 
[97] 12206335005 12206107003 12206108002 12206108001 12206335002 122068020 
[103] 12206121001 12206350001 12206119004 12206111004 12206335003 12206350004
[109] 12206100003 12206100002 12206100001 12206108003 12206118001 12206118003
[115] 12206109001 12206119003 12206118002 12206118003 12206109002 12206109003
[121] 122061220   12206116001 12206110001 12206109006 12206111001 12206109005
[127] 122068010   12206116002 12206109004 12206121002 12206110002 12206118004
[133] 12206111002 12206116003 12206116004 12206116000 12206110003 122063470 
[139] 12206112001 12206112002 12206116000 12206111003 12206112003 12206130002
[145] 12206130003 12206130001 122061230   12206116000 12206110004 122068040 
[151] 122068030   12206113001 12206113002 12206121003 122065600   122061240 
[157] 12206130004 12206114203 122068020   12206114202 122068040   12206130005
[163] 122068020   122068130   12206114201 122068040   12206111005 12206117001
[169] 122061141   122061141   12206117005 122068020   12206130006 12206114205
[175] 12206114204 122068030   12206130007 12206111006 122068050   12206116000
[181] 122068050   122068020   12206140001 122068110   12206130008 12206140002
[187] 12206117002 12206111007 122068120   12206117004 12206140005 12206140003
[193] 12206222001 12206140004 12206111008 122062230   12206117003 122062290 
[199] 12206224000 12206111009 122062270   122062250   12206222003 12206220000
[205] 12206222002 122068140   12206222004 122068060   12206220001 122065610 
[211] 12206230007 12206221001 12206230001 12206220002 12206222005 12206220003
[217] 12206226002 12206226001 12206220000 12206230005 12206230006 122062230 
[223] 12206220000 122062230   12206562001 12206220004 12206224001 12206224005
[229] 12206230002 12206226003 12206224002 122068090   12206221002 12206224003
[235] 122068110   12206221003 122068100   12206224004 12206224006 122068050 
[241] 12206226004 12206221006 12206230003 12206224000 12206226006 12206226007
[247] 12206221004 12206230004 12206221005 12206226005 122068070   12206562002
[253] 12206562003
228 Levels: 12206100001 12206100002 12206100003 12206101001 ... 122069990


122069990 が海域と判断したコード番号

※QGISでshapeファイルを読み「レイヤ」→「属性テーブルを開く」で
確認する方が簡単かもしれませんが,ここではRの学習として行います.

・参考サイト
Package ‘foreign’January 6, 2013』
http://cran.r-project.org/web/packages/foreign/foreign.pdf

【例2】
library(sp)
px <- runif(100,0,10)
py <- runif(100,0,10)
pz <- as.data.frame(px+py)
colnames(pz) <- c("pz")
pnt_xy <- cbind(px, py)
pnt_sp <- SpatialPoints(data.frame(px, py))
pnt_spdf <- SpatialPointsDataFrame(pnt_sp, pz)
plot(pnt_spdf, pch=20, cex=1.5, axes=TRUE, col="cyan")

Zu_ex2_1_a
・各オブジェクトの確認
px[1:5]
[1] 5.5363970 0.4466756 1.5218465 3.8153368 1.2822201

py[1:5]
[1] 6.839750 6.765926 4.789548 8.114962 4.502171

head(pz,n=4L)
         pz
1 12.376147
2  7.212601
3  6.311395
4 11.930298

head(pnt_xy,n=4L)
            px       py
[1,] 5.5363970 6.839750
[2,] 0.4466756 6.765926
[3,] 1.5218465 4.789548
[4,] 3.8153368 8.114962

head(pnt_sp,n=4L)
SpatialPoints:
            px       py
[1,] 5.5363970 6.839750
[2,] 0.4466756 6.765926
[3,] 1.5218465 4.789548
[4,] 3.8153368 8.114962
Coordinate Reference System (CRS) arguments: NA

pnt_spdf[1:4,]
          coordinates        pz
1   (5.5364, 6.83975) 12.376147
2 (0.446676, 6.76593)  7.212601
3  (1.52185, 4.78955)  6.311395
4  (3.81534, 8.11496) 11.930298


・pnt_spdfをspplotで表示
(spplotの使い方は理解していませんが覚書として記載)
pzの値を色分けで分類して表示.

spplot(pnt_spdf,xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y",scales=list(draw = TRUE),
col.regions=cm.colors(5))


Zu_ex2_2a
 spplotのパラメタは理解できていませんので,誤りがあると思います.
・参考サイト
Package ‘sp’December 21, 2012
http://cran.r-project.org/web/packages/sp/sp.pdf

======
以下メモ
cols <- grey.colors(5, 0.95, 0.55, 2.2)
spplot(pnt_spdf,xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y",
scales=list(draw = TRUE),col.regions=cols,cuts=5)

Gray_a

summary(px)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
0.009392 3.314000 5.502000 5.302000 7.426000 9.773000

summary(py)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.0684  2.2140  4.7650  4.8440  7.6640  9.9230

summary(pz)
       pz         
Min.   : 0.2976 
1st Qu.: 8.0831 
Median :10.3355 
Mean   :10.1461 
3rd Qu.:12.6310 
Max.   :16.7444 

summary(pnt_sp)
Object of class SpatialPoints
Coordinates:
           min      max
px 0.009391883 9.772706
py 0.068397804 9.922550
Is projected: NA
proj4string : [NA]
Number of points: 100

summary(pnt_spdf)
Object of class SpatialPointsDataFrame
Coordinates:
           min      max
px 0.009391883 9.772706
py 0.068397804 9.922550
Is projected: NA
proj4string : [NA]
Number of points: 100
Data attributes:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.2976  8.0830 10.3400 10.1500 12.6300 16.7400

●階級を指定すると
cuts <- c(0,3,6,9,12,15,17)
grys <- grey.colors(length(cuts), 0.90, 0.50, 2.2)
spplot(pnt_spdf,xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y",
scales=list(draw = TRUE),col.regions=grys,cuts=cuts)

Gray_b_1

・凡例を図の左に設定する方法は?

●plotでpzを円のサイズで表現する
plot(pnt_spdf,axes=TRUE,col="red",pch=1,cex  = sqrt(pnt_spdf$pz))
plot(pnt_spdf,col="green",pch=3,cex=1.5,add=TRUE)
legVals <- c(3, 6,10,14,17)
legend("topleft", legend = legVals, pch = 1, pt.cex = sqrt(legVals),
bty = "n", title  = "pz値")

Zu_en_1_a_2
・凡例の円が重なっているが,調整の方法?


2012年12月18日 (火)

RのreadLinesの使用例

■『次のテキストファイル(in.txt)から行番号とピリオドを除く』を例として
1.<!DOCTYPE html>
2.<html lang='en'>
3.<head>
4.    <meta charset='utf-8' />
5.    <title>My OpenLayers Map</title>
6.    <script type='text/javascript' src='OpenLayers.js'></script>
7.    <script type='text/javascript'>

●コード例
fin <- file("in.txt","r")      # open a input file connection
fout <- file("out.txt","w") # open an output file connection
a <- readLines(fin) #aはファイルの1行が1要素の文字列となるベクトル.ファイルを最後まで読む
for (i in 1:length(a)){
  Tmp <- a[i]
  start <- regexpr('\\.', Tmp)  # ピリオドの位置を返す.\\. はピリオドを表す.
  stop <- nchar(Tmp)
  Tmp1 <- substr(Tmp, start+1, stop)
  cat(paste(Tmp1,"\n",sep=""), file=fout, append=TRUE)
}
close(fin)
close(fout)

※参考
a <- readLines(fin)
"in.html"を最後まで読む.
ファイルの1行が1要素の文字列となるベクトルとしてa に代入する.

start <- regexpr('\\.', Tmp)
ピリオドの位置を返す.\\. はピリオドを表す.
1行にピリオドが複数ある場合は最初の位置を返す.

start <- regexpr('.', Tmp, fixed=TRUE) でもOKかも(未確認).
fixed=TRUEで ピリオドを正規表現と解釈しないことを指示.
デフォルトでは正規表現が有効になっているよう.

See Also: gregexpr

※参考サイト
データファイルの順次処理:ファイル名のいろいろな決め方
http://takenaka-akio.org/doc/r_auto/chapter_11.html
生物統計学 正規表現
http://r.livedocs.net/r/regex.html

●出力ファイル【out.txt】(行番号とピリオドを除いたテキストファイル)
<!DOCTYPE html>
<html lang='en'>
<head>
    <meta charset='utf-8' />
    <title>My OpenLayers Map</title>
    <script type='text/javascript' src='OpenLayers.js'></script>
    <script type='text/javascript'>



2012年4月 8日 (日)

Rコードにawkコードを組み込む

■Rコードにawkコードを組み込む
Rからawkを利用する.
memo_utf_8N.txt → 「memo_utf_8N.zip」
次の3行をawk_test.Rに書く.文字コードは環境に合わせてください.

system(' echo \' hello awk! \'  ')
system(' awk \' BEGIN { print "木更津高専" } \' ')
system(' awk \' { print NR,":"$0} \'  memo_utf_8N.txt ')

テキストファイルに行番号を付けて表示するコード例を含む.

●実行結果
$ Rscript awk_test.R   #Unix コンマンドラインで実行.$ Rscript awk_test.R  > res もOK
hello awk!
木更津高専
1 :Rのコードの中にawkのコードを書くことができるようです.
2 :system('your awk command here') のように書くようです.
3 :例として次の3行をRコードにかいて確認してみよう.
4 :このテキストに行番号を付けて表示する.
5 :
6 :system(' echo \' hello awk! \'  ')
7 :system(' awk \' BEGIN { print "木更津高専" } \' ')
8 :system(' awk \' { print NR,":",$0} \' memo_utf_8N.txt ')

【参考】
Awk script in R
http://stackoverflow.com/questions/9075208/awk-script-in-r

■awkをコンマンドラインから実行すると
$ awk '{print NR,":",$0}'  memo_utf_8N.txt   #awkをUnix のコンマンドラインから実行する.
1 : Rのコードの中にawkのコードを書くことができるようです.
2 : system('your awk command here') のように書くようです.
3 : 例として次の3行をRコードにかいて確認してみよう.
4 : このテキストに行番号を付けて表示する.
5 :
6 : system(' echo \' hello awk! \'  ')
7 : system(' awk \' BEGIN { print "木更津高専" } \' ')
8 : system(' awk \' { print NR,":",$0} \' memo_utf_8N.txt ')

●MSDOSコマンドプロンプトでもできると思いますが,
'  " の使い方がUnixと異なっているようです.
また文字コードも混乱してしまいました.
Unixでの作業がよいと思います.

●以下はパスして下さい.
MSDOSで実行例

■awk_test.R,memo_utf_8N.txtとも文字コードUTF_8N

C:\mywork>Rscript awk_test.R  #RscriptをMSDOSで実行
hello awk!
木更津高専
1 :Rのコードの中にawkのコードを書くことができるようです.
2 :system('your awk command here') のように書くようです.
3 :例として次の3行をRコードにかいて確認してみよう.
4 :このテキストに行番号を付けて表示する.
5 :
6 :system(' echo \' hello awk! \'  ')
7 :system(' awk \' BEGIN { print "木更津高専" } \' ')
8 :system(' awk \' { print NR,":",$0} \' memo_utf_8N.txt ')

■memo_shift_jis.txtはSHIFT-JIS 「memo_shift_jis.zip」

C:\mywork>gawk "{print NR,\":\",$0}"  memo_shift_jis.txt     #awkをMSDOSで実行
1 : Rのコードの中に,awkのコードを書くことができるようです.
2 : system('your awk command here') のように書くようです.
3 : 例として次の3行をRコードにかいて確認してみよう.
4 :
5 : system(' echo \' hello guy! \'  ')
6 : system(' awk \' BEGIN { print "木更津高専" } \' ')
7 : system(' awk \' { print } \'  memo_shift_jis.txt ')


2010年11月20日 (土)

図形描画関数群(http://aoki2.si.gunma-u.ac.jp/R/plot.html)

図形描画のための関数群を提供する。(ソースコード青木 繁伸先生、一部変更)
(http://aoki2.si.gunma-u.ac.jp/R/plot.html)
使用法・引数

  • plot.start(x1=0, y1=0, x2=20, y2=20, ...)
    プロット関数群の開始。
    (x1, y1)-(x2, y2) 描画領域の左下隅と右上隅の座標を宣言する。
    普通に R の関数(hist や plot などなど)を使って描かれたグラフに以下の関数を使って図形を付加する場合には不要。
    縦横比を1にするためには,asp=1 を加えること。

  • plot.line(x1, y1, x2, y2, ...)
    (x1, y1)-(x2, y2) を結ぶ直線を描画する。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。

  • plot.hline(x1, x2, y, ...)
    (x1, y)-(x2, y) を結ぶ水平な直線を描画する。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。

  • plot.vline(x, y1, y2, ...)
    (x, y1)-(x, y2) を結ぶ垂直な直線を描画する。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。

  • plot.box(x1, y1, x2, y2, ...)
    (x1, y1)-(x2, y2) を対角頂点とする長方形(正方形)を描画する。
    実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。

  • plot.circle(ox, oy, r, start=0, end=360, ...)
    中心を (ox, oy) とする,半径 r の円を描く。
    start, end には,描き始めと描き終わりの位置を指定できる。水平線を基準として,角度(度単位)で指定できる(0, 360 とすると,完全な円を描くことになる。90, 270 とすると,左半分の円を描くことを指示することになる)。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。

  • plot.circlef(ox, oy, r, ...)
    中心を (ox, oy) とする,半径 r の円を描き,内部を塗りつぶす。
    実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。

  • plot.ellipse(ox, oy, ra, rb, phi=0, start=0, end=360, ...)
    中心を (ox, oy) とする,長径 ra,短径 rb の楕円を描く。
    phi は楕円の傾き(長径が水平線となす角度),start, end には,描き始めと描き終わりの位置を指定できる。長径を基準として,角度(度単位)で指定できる(0, 360 とすると,完全な楕円を描くことになる。90, 270 とすると,左半分の楕円を描くことを指示することになる)。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。

  • plot.ellipsef(ox, oy, ra, rb, phi=0, ...)
    中心を (ox, oy) とする,長径 ra,短径 rb の楕円を描き,内部を塗りつぶす。
    phi は楕円の傾き(長径が水平線となす角度)。
    実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。

  • plot.polygon(x1, y1, l, n, ...)
    (x1, y1) を描きはじめとして,右方向へ水平
    に l の長さの辺を描きその後反時計回りに辺を構成して n 正方角形を描く
    。 実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。

  • plot.polygon2(ox, oy, r, n, phi=90, ...)
    (ox, oy) 中心とする円に内接する正 n 角形を描画する。
    最初の頂点の位置は phi 引数(度)によって,反時計回りに再設定できる。
    実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。

  • plot.grid(x1, y1, x2, y2, wx, wy=NULL, ...)
    (x1, y1)-(x2, y2) を対角頂点とする長方形(正方形)を描画し,横方向の間隔が wx, 縦方向の間隔が wy となる格子を描く。
    実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。
    長方形内部を塗りつぶすには,前もって plot.box で塗りつぶした長方形を描いてから plot.grid 関数を使う。

  • radian(degree)
    度をラジアンに変換する
無料ブログはココログ
2017年8月
    1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31