📈
自由な作図を提供するggplot2パッケージのレイヤー定義についていろいろ触ってみたまとめ
はじめに
これまでRの基礎的な話を延々としてきましたが、今回はパッケージ ggplot2 の話です。
ggplot2パッケージは自由にレイヤーを追加して、綺麗なグラフを簡単に作ることができる、Rでは人気のプロット用パッケージです。
ggplot2のようにレイヤー操作を可能にしているのが、 geomサブクラス による作図の設計です。ggplot2の代表的な関数を以下に見てみます。
ggplot2::geom_point
function (mapping = NULL, data = NULL, stat = "identity", position = "identity",
..., na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
{
layer(data = data, mapping = mapping, stat = stat, geom = GeomPoint,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...))
}
<bytecode: 0x55afe91ef088>
<environment: namespace:ggplot2>
geom_point 関数の内部では、レイヤーを作成する layer() 関数が呼び出されており、その引数には geom = GeomPoint と stat = "identity" が指定されています。これがサブクラスとして設定されている部分で、作図方法(ポイント、ライン、ポリゴン等)とデータの変換を決定しているところです。
それぞれ別の内部処理がされており、関数名で geom_hogehoge としているところでサブクラス名もわかります。
どのようなサブクラスが使用できるかは、公式のページで確認することができます。
多くのサブクラスが設計されていることがわかります。
普段はこれを組み合わせて自由に作っていきますが、凝った図や手間がかかる作図では、関数化してフォーマットにしたくなってきます。
一連の作業をまとめて関数化する方法も公式のページで紹介されています。
いつもやっている作業を無名関数 function() で囲って関数化すればいいよ、ということです。
今回はもうすこし踏み込んで、サブクラスを触りつつ、レイヤー操作は残したまま処理を簡素化する方法を探っていきます。
作図の確認
普段は使わないような図を作った方が説明として分かりやすいと思いますので、今回は、化学的な考察で使用するイオンバランスを図示する hexadiagram と piperdiagram を定義してみます。
それぞれどのようなものかは、以下のgoogle検索の結果をご確認ください。
- hexadiagram
- piperdiagram
元データは当然、イオン種別に列が作成され、地点別に整理されているはずです。
Name NaK Ca Mg SO4 NO3 HCO3 Cl
1 point.1 1.229 3.074 3.2 0.9 0.5 6.4 1.4
2 point.2 1.300 4.000 3.5 0.8 0.7 6.0 1.5
3 point.3 0.900 3.500 2.9 1.0 0.4 5.9 1.2
作図自体は難しいものではありませんが、プロットを作成するには少し面倒な作業です。
このような場合は関数化するメリットが見えてきます。
データの準備
まずは今回使用するサンプルデータを作成しておきます。
data.frameの列名は、地点名以外はそれぞれイオン名です。
データは meq/L 単位で作成しています。
ggplotで使うときは、データを x, y が指定できるように tidyr::pivot_longer 等を使用して変換する必要がありますので、データを準備したら変換します。
library(dplyr)
library(ggplot2)
data <- data.frame(Name = c("point.1", "point.2", "point.3"),
NaK = c(1.229, 1.3, 0.9),
Ca = c(3.074, 4, 3.5),
Mg = c(3.2, 3.5, 2.9),
SO4 = c(0.9, 0.8, 1.0),
NO3 = c(0.5, 0.7, 0.4),
HCO3 = c(6.4, 6, 5.9),
Cl = c(1.4, 1.5, 1.2))
str(data)
'data.frame': 3 obs. of 8 variables:
$ Name: chr "point.1" "point.2" "point.3"
$ NaK : num 1.23 1.3 0.9
$ Ca : num 3.07 4 3.5
$ Mg : num 3.2 3.5 2.9
$ SO4 : num 0.9 0.8 1
$ NO3 : num 0.5 0.7 0.4
$ HCO3: num 6.4 6 5.9
$ Cl : num 1.4 1.5 1.2
data <- tidyr::pivot_longer(data = within(data, {NO3 <- NO3 + SO4}),
cols = !Name,
names_to = "z",
values_to = "x") %>%
mutate(x = if_else(z %in% c("NaK", "Ca", "Mg"), -x, x)) %>%
mutate(y = do.call(c, lapply(as.list(z), function(m) {
switch(m,
"NaK" = 3,
"Ca" = 2,
"Mg" = 1,
"SO4" = 1,
"NO3" = 1,
"HCO3" = 2,
"Cl" = 3,
"SO4+NO3" = 1
)})))
str(data)
tibble [21 × 4] (S3: tbl_df/tbl/data.frame)
$ Name: chr [1:21] "point.1" "point.1" "point.1" "point.1" ...
$ z : chr [1:21] "NaK" "Ca" "Mg" "SO4" ...
$ x : num [1:21] -1.23 -3.07 -3.2 0.9 1.4 ...
$ y : num [1:21] 3 2 1 1 1 2 3 3 2 1 ...
hexadiagramの作図例
まずは、通常ggplot2を使うとどのようになるか確認します。
イオン名はy軸、x軸はイオンの当量濃度です。
## ポリゴンの作成
poly_1 <- within(subset(data, !z %in% "NO3"), {
group <- "Main Group"
})
poly_2 <- within(subset(data, z %in% c("SO4", "NO3", "HCO3")), {
group <- "NO3 Group"
z <- if_else(z == "NO3", "SO4+NO3", z)
})
d <- rbind(poly_1, poly_2)
## 作図
p <- ggplot(d) +
geom_polygon(aes(x = x, y = y, group = group, fill = group), col = "black") +
geom_segment(data = data.frame(x = 0, y = 1, xend = 0, yend = 3),
aes(x = x, y = y, xend = xend, yend = yend),
col = "black") +
scale_fill_manual(values = c("NA", "gray")) +
facet_wrap(vars(Name)) +
labs(x = "meq/L",
y = "",
title = "HexaDiagram") +
coord_equal() +
theme(legend.position = "none",
axis.text.y = element_blank())
p
このような図になる関数を作成していきます。
関数化の注意点(ggplot2パッケージ)
ggplot2パッケージを使って図を作成する際、 + 演算子を使ってレイヤーを重ねるようにします。
自前で関数を定義するときは、この + 演算子の使い方は要注意です。
既存の ggplot2関数 を使って + 演算子を使って関数化することはできますが、関数内部では + 演算子を使って処理をまとめる操作ができる場合とできない場合があります。
その理由が、ggplot2ではクラスの継承が内部で行われているためであり、 + 演算子を使う場合は、前から順番にしか処理ができない仕様になっています。
つまり、 ggplot() から処理をひとまとめにした関数は実行できますが、 geom_hogehoge() だけを関数化し、なおかつ + 演算子で処理を続行する場合は失敗します。
## 成功する例
geom_hogehoge <- function(...) {
geom_point(...)
}
ggplot(data) +
geom_hogehoge(...) +
theme(...)
## 成功する例
createGraph <- function(...) {
ggolot(data) +
geom_point(...) +
theme(...)
}
createGraph()
## 失敗する例
geom_hogehoge <- function(...) {
geom_point(...) +
theme(...)
}
ggplot(data) +
geom_hogehoge(...)
レイヤー1つを追加するラッパー関数を作成することは難しくありませんが、複雑な操作を実行する場合はサブクラスの作成等が必要です。
ggplot2公式ホームページにどのような記載があるかと言うと
ggproto() を使ってサブクラスを作ってね、と書いてあります。
筆者がいろいろ触ってみましたが、何かしらのエラーに当たり、時間ばかりが取られる結果になりました。そのため、以下を参考に関数定義のレベルを整理すると良いと思います。
| 定義レベル | やりたいこと |
|---|---|
| 1. 前処理でデータを作成 | 軸(x, y)が定まっていないデータをプロットする |
| 2. 関数内でggproto(デフォルトの読み込みを変更)上書き定義 | 既存のサブクラスをもとに、軸やポイント形状等、aes()等で定義していることを初期値化する |
| 3. ggproto定義(stat) | 軸(x, y)は決まっているが、処理を施して軸の値を変更する等 |
| 4. ggproto定義(geom)::gridパッケージを使用して作図を定義 | 作図方法を定義する場合 |
はじめから 3., 4. のサブクラスを定義するとかなりの労力が必要となるため、あまりおすすめはしません。
それよりも、既存のサブクラスを上書きする 2. の方法が簡単で便利でした。
今回の例のように、座標軸(x, y)を計算で定義しなければなならない場合は 1. の前処理も定義します。
データの継承について
ggplot() で指定したデータは、サブクラス内でのみ引用することができます。
そのため、クラスや環境の継承がされていない場所では引用できません。
データを ggplot() 関数から継承して使用する場合は stat のサブクラスを定義する必要があります。
geom_hogehoge(data = NULL, ...) {
## 前処理
...
## 作図
layer(data = data, ...)
}
## 失敗例
ggplot(data) +
geom_hogehoge(aes(x, y))
## 成功例
ggplot() +
geom_hogehoge(data = data, aes(x, y))
このような使用上の注意点もあるため、ggplot2公式で紹介されているように、 ggplot() 関数から、まるごと関数化する方法が安全で簡単でしょう。
それでも、レイヤーごとで操作できると、ggplot2の1操作1レイヤーの仕様を残すことができるため、実施にプロットするときの操作上の不都合は起きにくくなります。
また、複数のデータを重ね合わせることができる利点ありますので、事前にデータをまとめる必要がなく、ある程度の柔軟性が残ります。
特にデータ解析を進める上ではこのような柔軟性は必要だと筆者は思いますので、今回はレイヤーを個別に生成する関数として作成していきます。
HexaDiagramのレイヤー作図関数
前処理はめんどくさい作業なので、なるべく簡単に使うことができるようにします。
最初に geom_point で確認したとおり、レイヤー呼び出しの関数内部ではレイヤーを作成する layer() 関数が呼び出されています。
前処理はレイヤー呼び出し関数内で、 layer() が実行される前の段階で行う方法場最も簡単です。
レイヤ-呼び出し関数を geom_Hexadiag() として定義すると、以下のようになります。
library(ggplot2)
library(dplyr)
## 使用するデータの準備
data <- data.frame(point = c("point.1", "point.2", "point.3"),
NaK = c(1.229, 1.3, 0.9),
Ca = c(3.074, 4, 3.5),
Mg = c(3.2, 3.5, 2.9),
SO4 = c(0.9, 0.8, 1.0),
NO3 = c(0.5, 0.7, 0.4),
HCO3 = c(6.4, 6, 5.9),
Cl = c(1.4, 1.5, 1.2))
## レイヤー呼び出し関数定義
geom_Hexadiag <- function(mapping = NULL, data = NULL, pname = NULL,
stat = "identity", position = "identity",
rule = "evenodd", ..., na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE) {
## データの前処理 ==================================
## データの読み込み
force(data)
## 地点の列名を個別に設定
pname <- substitute(pname)
## 前処理関数
hexadiag <- function(data, pname) {
library(dplyr)
library(tidyr)
force(data)
pname <- deparse(substitute(pname))
names(data) <- ifelse(names(data) == pname, "Name", names(data))
ion_name <- c("NaK", "Ca", "Mg", "SO3", "NO3", "HCO3", "Cl")
data <- data %>%
mutate(NO3 = NO3 + SO4) %>%
tidyr::pivot_longer(cols = !Name,
names_to = "ion",
values_to = "x") %>%
mutate(x = if_else(ion %in% c("NaK", "Ca", "Mg"), -x, x)) %>%
#mutate(z = factor(z, levels = ion_name)) %>%
mutate(y = do.call(c, lapply(as.list(ion), function(m) {
switch(m,
"NaK" = 3,
"Ca" = 2,
"Mg" = 1,
"SO4" = 1,
"NO3" = 1,
"HCO3" = 2,
"Cl" = 3,
1)})))
poly_1 <- subset(data, !ion %in% "NO3") %>%
mutate(group = "Main Group") %>%
mutate(fill = "NA")
poly_2 <- subset(data, ion %in% c("SO4", "NO3", "HCO3")) %>%
mutate(group = "NO3 Group") %>%
mutate(ion = if_else(ion == "NO3", "SO4+NO3", ion)) %>%
mutate(fill = "gray")
rbind(poly_1, poly_2)
}
## データを作成
grid <- eval(bquote(hexadiag(data, .(pname))))
## サブクラスを修正(GeomPolygonクラスを上書き修正) ========================
geomHexaDiag <- ggproto("geomHexaDiag", GeomPolygon,
default_aes = aes(colour = "black", fill = grid$fill, size = 0.5,
linetype = 1, alpha = NA, subgroup = NULL))
## レイヤー定義 ======================================================
layer(data = grid,
mapping = aes(x = x, y = y, group = group),
stat = stat, geom = geomHexaDiag,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, rule = rule, ...))
}
## 実行
ggplot() +
geom_Hexadiag(data = data, pname = point) +
#ylim(-0.5, 4) +
facet_wrap(vars(Name)) +
labs(x = "meq/L",
y = "",
title = "HexaDiagram") +
coord_equal() +
theme_bw() +
theme(legend.position = "none",
axis.text.y = element_blank())
geom_HexaDiag() は六角形のポリゴンを描画する関数ですが、内部でデータの変換、サブクラスの定義、レイヤー定義を行っています。
データの変換の後はサブクラスを定義せずに geom_point() を呼び出しても同じ結果が得られますが、柔軟さを出すためにこのようにしています。
piper-diagram（トリリニアダイアグラム）の作図関数
今度は作図領域の作成からです。
piper-diagramは、2軸プロットではなく、イオン種3要素の割合によるプロットがされている点が特殊です。更に、3要素プロットされたポイント(カチオン、アニオン)を元に、中央領域に混合割合のプロットがされる必要があります。
これらの計算はめんどくさいので、やはり前処理に入れ込んでしまいます。
作図領域の描画は一度きりのデータ定義とした方がスマートですので、レイヤー定義関数を定義したら、内部で local() を使ってデータを作成します。
続いてサブクラスを定義して、レイヤー定義を行っており、流れは先程の geom_HexaDiag() と同じです。
しかしながら、今回は作図領域描画、ラベル作成、データをプロットする3段階があるため、それぞれ1レイヤーとして関数を作成します。
library(ggplot2)
library(tidyr)
data <- data.frame(point = c("point.1", "point.2", "point.3"),
NaK = c(1.229, 1.3, 0.9),
Ca = c(3.074, 4, 3.5),
Mg = c(3.2, 3.5, 2.9),
SO4 = c(0.9, 0.8, 1.0),
NO3 = c(0.5, 0.7, 0.4),
HCO3 = c(6.4, 6, 5.9),
Cl = c(1.4, 1.5, 1.2))
## 作図領域枠を作成する関数
geom_PiperArea <- function(stat = "identity",
position = "identity", ..., arrow = NULL, arrow.fill = NULL,
lineend = "butt", linejoin = "round", na.rm = FALSE,
show.legend = NA, inherit.aes = TRUE) {
## 作図領域描画データ作成
grid <- local({
## カチオン側
t1 <- data.frame(x = c(0, 0, 50),
xend = c(100, 50, 100),
y = c(0, 0, 50 * sqrt(3)),
yend = c(0, 50 * sqrt(3), 0),
size = 1,
group = "t1")
## アニオン側
t2 <- within(t1, {
x <- x + 120
xend <- xend + 120
group <- "t2"
})
## 中央領域
dia <- data.frame(x = rep(110, 4),
xend = c(60, 160, 60, 160),
y = c(10 * sqrt(3), 10 * sqrt(3),
rep(10 * sqrt(3) + 50 * sqrt(3) * 2, 2)),
yend = rep(10 * sqrt(3) + 50 * sqrt(3), 4),
size = 1,
group = "dia")
## カチオン側細軸
l1 <- data.frame(x = c(seq(20, 80, length = 4), seq(20, 80, length = 4), seq(10, 40, length = 4)),
xend = c(seq(10, 40, length = 4), seq(60, 90, length = 4),
rev(seq(60, 90, length = 4))),
y = c(rep(0, 8), seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
yend = c(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)],
rev(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
size = 0.25,
group = "l1")
## カチオン側細軸
l2 <- within(l1, {
x <- x + 120
xend <- xend + 120
group <- "l2"
})
## 中央領域細軸
ld <- data.frame(x = c(seq(70, 100, length = 4), seq(120, 150, length = 4)),
xend = c(seq(120, 150, length = 4), seq(70, 100, length = 4)),
y = c(rev(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)]),
seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)]),
yend = c(rev(seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
size = 0.25,
group = "ld")
## データを結合してreturn
rbind(t1, t2, dia, l1, l2, ld)
})
## サブクラスの定義(geom_segmentを上書き修正)
geomPiperArea <- ggproto("geomArea", GeomSegment,
default_aes = aes(colour = "black",
size = grid$size,
linetype = 1, alpha = NA))
## レイヤー定義
layer(data = grid,
mapping = aes(x = x, y = y, xend = xend, yend = yend),
stat = stat, geom = geomPiperArea,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(arrow = arrow, arrow.fill = arrow.fill,
lineend = lineend, linejoin = linejoin, na.rm = na.rm,
...))
}
## 描画領域の軸ラベルを作成する関数
geom_PiperAreaLabel <- function(stat = "identity", position = "identity",
..., nudge_x = 0, nudge_y = 0, check_overlap = FALSE,
na.rm = FALSE, show.legend = NA, inherit.aes = TRUE) {
## ラベルデータの作成
grid_label <- local({
## カチオン軸ラベル
pl1 <- data.frame(
x = c(seq(80, 20, length = 4),
(seq(10, 40, length = 4) - 5),
(seq(60, 90, length = 4) + 5)),
y = c(rep(-5, 4),
seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)],
rev(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)])),
label = rep(seq(20, 80, length = 4), 3),
size = 4,
angle = 0)
## アニオン軸ラベル
pl2 <- within(pl1, {
x <- c(seq(20, 80, length = 4),
(seq(90, 60, length = 4) + 5),
(seq(40, 10, length = 4) - 5)) + 120
})
## 中央領域軸ラベル
pl3 <- data.frame(
x = c((seq(100, 70, length = 4) - 5),
(seq(70, 100, length = 4) - 5),
(seq(120, 150, length = 4) + 5),
(seq(150, 120, length = 4) + 5)),
y = c(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)],
seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)],
rev(seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
rev(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)])),
label = c(seq(80, 20, length = 4),
seq(20, 80, length = 4),
seq(80, 20, length = 4),
seq(20, 80, length = 4)),
size = 4,
angle = 0)
## イオン名ラベル
l <- data.frame(
x = c(50, 19, 81,
(c(50, 19, 81) + 120),
72.5, 147.5),
y = c(-10, 50, 50,
-10, 50, 50,
150, 150),
label = c("Ca^'2+'", "Mg^'2+'", "Na^'+'~+~K^'+'",
"Cl^'-'", "HCO[3]^'-'", "SO[4]^'2-'",
"SO[4]^'2-'~+~Cl^'-'", "Ca^'2+'~+~Mg^'2+'"),
size = 5,
angle = c(0, 60, -60, 0, 60, -60, 60, -60))
## データを結合してretuen
rbind(pl1, pl2, pl3, l)
})
## サブクラスの作成(geomTextを上書き修正)
geomPiperAreaLabel <- ggproto("geomPiperAreaLabel", GeomText,
default_aes = aes(colour = "black",
size = grid_label$size,
angle = grid_label$angle,
hjust = 0.5, vjust = 0.5, alpha = NA,
family = "", fontface = 1,
lineheight = 1.2
))
## レイヤー定義
layer(data = grid_label,
mapping = aes(x = x, y = y, label = label),
stat = stat, geom = geomPiperAreaLabel,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(parse = TRUE, check_overlap = check_overlap,
na.rm = na.rm, ...))
}
## イオンデータをプロットする関数
geom_PiperDiag <- function(mapping = NULL, data = NULL,
pname = NULL, NaK = NULL, Ca = NULL, Mg = NULL,
SO4 = NULL, NO3 = NULL, HCO3 = NULL, Cl = NULL,
stat = "identity", position = "identity",
..., na.rm = FALSE, show.legend = NA, inherit.aes = TRUE) {
force(data)
NaK <- substitute(NaK)
Ca <- substitute(Ca)
Mg <- substitute(Mg)
SO4 <- substitute(SO4)
NO3 <- substitute(NO3)
HCO3 <- substitute(HCO3)
Cl <- substitute(Cl)
pname <- substitute(pname)
## ポイントデータをx, y軸のデータに修正
grid <- with(data, {
## 割合計算
t.cation <- eval(sum(NaK, Ca, Mg))
t.anion <- eval(sum(SO4, NO3, HCO3, Cl))
## 当量濃度を割合(%)に変換
base <- data.frame(Ca = eval(Ca) / t.cation,
Mg = eval(Mg) / t.cation,
Cl = eval(Cl) / t.anion,
SO4 = eval(SO4 / t.anion)) * 100
## 座標値(x, y)に変換
calc.cation <- data.frame(x = (100 - base$Ca) + cos(pi * 2 / 3) * base$Mg,
y = sin(pi * 1 / 3) * base$Mg)
calc.anion <- data.frame(x = ((base$Cl) + cos(pi * 1 / 3) * base$SO4) + 120,
y = sin(pi * 1 / 3) * base$SO4)
## カチオン側方程式 傾き: sqrt(3), 切片b1
## アニオン側方程式 傾き: -sqrt(3), 切片b2
## 2つの直線が交わる点=中央領域のプロット位置として計算
## 切片の計算
calc.cation$intercept <- calc.cation$y - sqrt(3) * calc.cation$x
calc.anion$intercept <- calc.anion$y + sqrt(3) * calc.anion$x
## 中央エリアの座標位置計算
## カチオン、アニオンそれぞれの方程式から解を求める連立方程式を組む
calc.centre <- as.data.frame(do.call(rbind,
lapply(data.frame(t(cbind(calc.cation$intercept,
calc.anion$intercept))),
function(m) {
## カチオン: -ax + y = b
## アニオン: ax + y = b
A <- matrix(c(-sqrt(3), sqrt(3), 1, 1), ncol = 2)
## 切片
b <- m
## 連立方程式の解を計算
solve(A, b)
})))
names(calc.centre) <- c("x", "y")
out <- data.frame(rbind(calc.cation[c("x", "y")],
calc.anion[c("x", "y")],
calc.centre))
out$group <- eval(pname)
out
})
## サブクラスの定義(geomPointを上書き修正)
geomPiperPlot <- ggproto("geomPiperPlot", GeomPoint,
default_aes = aes(shape = 21, colour = "black",
size = 5, fill = NA, alpha = 0.7,
stroke = 0.5))
## レイヤー定義
layer(data = grid,
mapping = aes(x = x, y = y, group = group, fill = group),
stat = stat, geom = geomPiperPlot,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...))
}
ggplot() +
geom_PiperArea() +
geom_PiperAreaLabel() +
geom_PiperDiag(data = data, pname = point, NaK = NaK, Ca = Ca, Mg = Mg,
SO4 = SO4, NO3 = NO3, HCO3 = HCO3, Cl = Cl) +
coord_equal() +
theme_void()
無事に描画できました。
終わりに
複雑な描画でも、関数にして処理を簡素化できました。
残念ながら、定義した関数をそのまま ggplotly() でplotlyの図に変換することはできませんでした。
plotlyにすると、Javascriptによりインタラクティブな操作が可能になりますが、自作関数では無理なようです。
これには、必ずggplot2パッケージ定義されている関数(geom_point(), geom_polygon(), ...)を使用する必要があるようです。
ggplotly() にも対応させるには、上記 3., 4. のように ggproto() での処理に任せる必要があります。
おまけ サブクラスを作成してプロットする方法
library(ggplot2)
library(tidyr)
data <- data.frame(point = c("point.1", "point.2", "point.3"),
NaK = c(1.229, 1.3, 0.9),
Ca = c(3.074, 4, 3.5),
Mg = c(3.2, 3.5, 2.9),
SO4 = c(0.9, 0.8, 1.0),
NO3 = c(0.5, 0.7, 0.4),
HCO3 = c(6.4, 6, 5.9),
Cl = c(1.4, 1.5, 1.2))
statPiperArea <- ggproto("statPiperArea", Stat,
#default_aes = aes(x = stat(x), y = stat(y), xend = stat(xend), yend = stat(yend)),
compute_group = function(data, scales) {
local({
## カチオン側
t1 <- data.frame(x = c(0, 0, 50),
xend = c(100, 50, 100),
y = c(0, 0, 50 * sqrt(3)),
yend = c(0, 50 * sqrt(3), 0),
size = 1,
group = "t1")
## アニオン側
t2 <- within(t1, {
x <- x + 120
xend <- xend + 120
group <- "t2"
})
## 中央領域
dia <- data.frame(x = rep(110, 4),
xend = c(60, 160, 60, 160),
y = c(10 * sqrt(3), 10 * sqrt(3),
rep(10 * sqrt(3) + 50 * sqrt(3) * 2, 2)),
yend = rep(10 * sqrt(3) + 50 * sqrt(3), 4),
size = 1,
group = "dia")
## カチオン側細軸
l1 <- data.frame(x = c(seq(20, 80, length = 4), seq(20, 80, length = 4), seq(10, 40, length = 4)),
xend = c(seq(10, 40, length = 4), seq(60, 90, length = 4),
rev(seq(60, 90, length = 4))),
y = c(rep(0, 8), seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
yend = c(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)],
rev(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)]),
size = 0.25,
group = "l1")
## カチオン側細軸
l2 <- within(l1, {
x <- x + 120
xend <- xend + 120
group <- "l2"
})
## 中央領域細軸
ld <- data.frame(x = c(seq(70, 100, length = 4), seq(120, 150, length = 4)),
xend = c(seq(120, 150, length = 4), seq(70, 100, length = 4)),
y = c(rev(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)]),
seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)]),
yend = c(rev(seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
size = 0.25,
group = "ld")
## データを結合してreturn
rbind(t1, t2, dia, l1, l2, ld)
})
}
)
statPiperAreaLabel <- ggproto("statPiperAreaLabel", Stat,
default_aes = aes(
#x = stat(x), y = stat(y),
colour = "black", hjust = 0.5, vjust = 0.5, alpha = NA,
family = "", fontface = 1,
lineheight = 1.2
),
compute_group = function(data, scales) {
local({
## カチオン軸ラベル
pl1 <- data.frame(
x = c(seq(80, 20, length = 4),
(seq(10, 40, length = 4) - 5),
(seq(60, 90, length = 4) + 5)),
y = c(rep(-5, 4),
seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)],
rev(seq(0, 50 * sqrt(3), length = 6)[-c(1, 6)])),
label = rep(seq(20, 80, length = 4), 3),
size = 4,
angle = 0)
## アニオン軸ラベル
pl2 <- within(pl1, {
x <- c(seq(20, 80, length = 4),
(seq(90, 60, length = 4) + 5),
(seq(40, 10, length = 4) - 5)) + 120
})
## 中央領域軸ラベル
pl3 <- data.frame(
x = c((seq(100, 70, length = 4) - 5),
(seq(70, 100, length = 4) - 5),
(seq(120, 150, length = 4) + 5),
(seq(150, 120, length = 4) + 5)),
y = c(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)],
seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)],
rev(seq(10 * sqrt(3) + 50 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3) * 2,
length = 6)[-c(1, 6)]),
rev(seq(10 * sqrt(3), 10 * sqrt(3) + 50 * sqrt(3), length = 6)[-c(1, 6)])),
label = c(seq(80, 20, length = 4),
seq(20, 80, length = 4),
seq(80, 20, length = 4),
seq(20, 80, length = 4)),
size = 4,
angle = 0)
## イオン名ラベル
l <- data.frame(
x = c(50, 19, 81,
(c(50, 19, 81) + 120),
72.5, 147.5),
y = c(-10, 50, 50,
-10, 50, 50,
150, 150),
label = c("Ca\u00B2\u207A", "Mg\u00B2\u207A", "Na\u207A+K\u207A",
"Cl\u207B", "HCO\u2083\u207B", "SO\u2084\u00B2\u207B",
"SO\u2084\u00B2\u207B+Cl\u207B", "Ca\u00B2\u207A+Mg\u00B2\u207A"),
size = 5,
angle = c(0, 60, -60, 0, 60, -60, 60, -60))
## データを結合してretuen
rbind(pl1, pl2, pl3, l)
})
})
table2Piperdata <- function(data = NULL,
pname = NULL, NaK = NULL, Ca = NULL, Mg = NULL,
SO4 = NULL, NO3 = NULL, HCO3 = NULL, Cl = NULL) {
force(data)
NaK <- substitute(NaK)
Ca <- substitute(Ca)
Mg <- substitute(Mg)
SO4 <- substitute(SO4)
NO3 <- substitute(NO3)
HCO3 <- substitute(HCO3)
Cl <- substitute(Cl)
pname <- substitute(pname)
## ポイントデータをx, y軸のデータに修正
grid <- with(data, {
## 割合計算
t.cation <- eval(sum(NaK, Ca, Mg))
t.anion <- eval(sum(SO4, NO3, HCO3, Cl))
## 当量濃度を割合(%)に変換
base <- data.frame(Ca = eval(Ca) / t.cation,
Mg = eval(Mg) / t.cation,
Cl = eval(Cl) / t.anion,
SO4 = eval(SO4 / t.anion)) * 100
## 座標値(x, y)に変換
calc.cation <- data.frame(x = (100 - base$Ca) + cos(pi * 2 / 3) * base$Mg,
y = sin(pi * 1 / 3) * base$Mg)
calc.anion <- data.frame(x = ((base$Cl) + cos(pi * 1 / 3) * base$SO4) + 120,
y = sin(pi * 1 / 3) * base$SO4)
## カチオン側方程式 傾き: sqrt(3), 切片b1
## アニオン側方程式 傾き: -sqrt(3), 切片b2
## 2つの直線が交わる点=中央領域のプロット位置として計算
## 切片の計算
calc.cation$intercept <- calc.cation$y - sqrt(3) * calc.cation$x
calc.anion$intercept <- calc.anion$y + sqrt(3) * calc.anion$x
## 中央エリアの座標位置計算
## カチオン、アニオンそれぞれの方程式から解を求める連立方程式を組む
calc.centre <- as.data.frame(do.call(rbind,
lapply(data.frame(t(cbind(calc.cation$intercept,
calc.anion$intercept))),
function(m) {
## カチオン: -ax + y = b
## アニオン: ax + y = b
A <- matrix(c(-sqrt(3), sqrt(3), 1, 1), ncol = 2)
## 切片
b <- m
## 連立方程式の解を計算
solve(A, b)
})))
names(calc.centre) <- c("x", "y")
out <- data.frame(rbind(calc.cation[c("x", "y")],
calc.anion[c("x", "y")],
calc.centre))
out$group <- eval(pname)
out
})
grid
}
p <- ggplot(table2Piperdata(data = data, pname = point, NaK = Nak, Ca = Ca, Mg = Mg,
SO4 = SO4, NO3 = NO3, HCO3 = HCO3, Cl = Cl)) +
geom_segment(stat = statPiperArea) +
geom_text(stat = statPiperAreaLabel) +
geom_point(aes(x = x, y = y, fill = group), shape = 21, size = 5, col = "black") +
coord_equal() +
theme_void()
p
## plotlyでインタラクティブな作図
## zennへの載せ方がわからなかったので、皆さんでお試しください。
plotly::ggplotly(p)
Discussion