[R]対数軸グラフを描く

R

グラフを描くときに、データがべき乗則に乗っている場合などでは対数軸を使用することが多いかと思います。例えば地震の規模(マグニチュード)別の回数ではグーテンベルグ・リヒター則が知られています。気象庁の震源データを使って、マグニチュード別の地震数を描いてみます。

setwd("D:/WorkingDir/path/")
par(mar=c(3, 3, 1, 1)) #余白設定
par(mgp=c(2, 0.7, 0))  #グラフ枠と軸との間隔設定
par(xaxs = "i")        #データ範囲とグラフの範囲(X軸)
par(yaxs = "i")        #データ範囲とグラフの範囲(Y軸)

#データ読み込み
dfEQ <- read.csv("h2017.csv", header=T)

#プロット
plot(dfEQ$M, dfEQ$N,
 log="y", #片対数
 xlim=c(0, 8), ylim=c(1, 1e+6), # xとyの範囲
 xlab="M", ylab="N",
 pch=20, type="o",
 lty=1, lwd=1
)

8行目で読み取っている「h2017.csv」は気象庁のデータからマグニチュードごとに地震数を集計したもので、下記のようなデータになっています。

元のデータ(h2017.csv)

元のデータ(h2017.csv)

Rで描いた片対数グラフ①

Rで描いた片対数グラフ①

13行目のように「log=”x”」とか「log=”y”」と指定すれば、指定した軸を対数軸にすることが出来ますが、上図のように10の累乗にしか目盛り線が描画されません。これでは見てくれが悪い、というか値を読みづらい場合があるのでなんとかならないかと検索していたところ、三重大学・奥村先生のサイトに「グラフの例：物理論文風グラフ」というそのものズバリなページがありました。
結果的には、下記のようにRでゴリゴリ書いて補助目盛り線を描くしかなさそうです。
ただ、軸目盛りを描くときに「expression(10^{0})」とすることで、乗数を右肩の小文字で書くことが出来るようです(37-39行目)。そして、48・50行目のように、sapply関数を使うことによってFor文を使わずに繰り返して補助目盛り線を描いています。

setwd("D:/WorkingDir/path/")

# 余白・隙間設定
par(mar=c(3, 3, 1, 1)) #余白設定
par(mgp=c(2, 0.7, 0))  #グラフ枠と軸との間隔設定
par(xaxs = "i")        #データ範囲とグラフの範囲(X軸)
par(yaxs = "i")        #データ範囲とグラフの範囲(Y軸)

# プロットするデータ
# 出典(気象庁) https://www.data.jma.go.jp/svd/eqev/data/bulletin/hypo.html
dfEQ <- read.csv("h2017.csv", header=T)

# グラフの範囲
iLogL <- 0
iLogU <- 6

#まずは空白グラフ
plot(NULL,
 log="y", #片対数
 xlim=c(-1, 8), ylim=c(10^iLogL, 10^iLogU), # xとyの範囲
 xaxt="n", yaxt="n",  # x軸もy軸も描かない（後で描く）
 xlab="", ylab=""
)
#軸ラベル
mtext("M", side=1, line=1.5)
mtext("N/n", side=2, line=2)

#X軸:ここは普通の線形軸
axis(side=1,        #side1:下
 at=seq(-1,8,1),     #0から8まで1ずつ
 tck=0.03,          #長さ0.03のティック
 labels=seq(-1,8,1), #0から8まで1刻みでラベル出力
 mgp=c(1,0.5,0)
)

#Y軸:目盛
sLab <- c(
 expression(10^{0}),expression(10^{1}),expression(10^{2}),expression(10^{3}),expression(10^{4}),expression(10^{5}),expression(10^{6})
)
axis(side=2,          #side2:左
 at=10^(iLogL:iLogU), #0から8まで1ずつ
 tck=0.03,            #長さ0.03のティック
 labels=sLab,
 mgp=c(1,0.5,0)
)

#Y軸:目盛
fPow <- function(x){(2:9)*10^x}
axis(side=2,        #side2:左
 at=sapply(iLogL:iLogU, fPow), #繰り返し(2:9)×10^(iLogL:iLogU)
 tck=0.01,          #長さ0.01のティック
　labels=FALSE,      #ラベル出力なし
 mgp=c(1,0.5,0)
)

#肝心のプロット
points(dfEQ$M, dfEQ$N,
 pch=20, type="o",
 lty=1, lwd=1
)
points(dfEQ$M, dfEQ$n,
 pch=21, type="o",
 lty=1, lwd=1, col=gray(0.6)
)

#凡例
legend("topright", 
 c("N(積算地震数)", "n(M区間地震数)"),
 pch=c(20, 21),
 lty=c(1, 1),
 col=c("black", gray(0.6))
)
Rで描いた片対数グラフ②

Rで描いた片対数グラフ②

ちょっと欲をかいて、折れ線を増やしたりしていますけど、やっていることは簡単ながら長ったらしいコードになっています。
ということで、仕事中とかで自分でコピペする用のメモでした。

この記事の作成には、下記のサイトを参考(というよりほとんどそのまま)にさせて頂きました。
グラフの例：物理論文風グラフ(三重大学奥村先生のサイト)