Match Y dual axis values in R-CodePudding

I am aware that it is discouraged to plot 2 separate y-axes on a single plot, and that it can be confusing and misleading.

Is there a way to match the values of the Y axis ? Such as finding a way to shift the right Y axis data within the graph (?) For example, I would like 40 in the left Y axis to match -1 in the right Y axis, because those two values correspond to the beginning of drought conditions categories.

Sample:

structure(list(Time = structure(c(9862, 9893, 9921, 9952, 9982, 
10013, 10043, 10074, 10105, 10135, 10166, 10196, 10227, 10258, 
10286, 10317, 10347, 10378, 10408, 10439), class = "Date"), Year = c(1997, 
1997, 1997, 1997, 1997, 1997, 1997, 1997, 1997, 1997, 1997, 1997, 
1998, 1998, 1998, 1998, 1998, 1998, 1998, 1998), VCI = c(48.7488482440362, 
51.8662335250972, 54.4212125411374, 61.7338808190779, 63.9870065731148, 
61.3375357670741, 62.6197335631611, 63.0950799644754, 61.6276947895731, 
61.1298324406371, 64.4422427513358, 60.3823204404222, 60.5883337239537, 
61.8918834440238, 59.1304135098709, 62.1668350554954, 61.9352586665065, 
55.75795384214, 50.3371363875305, 52.5748728440737), TCI = c(53.7071192239754, 
53.6178820221828, 57.7831310561669, 57.3996088686986, 49.8613200877384, 
54.9673093834738, 42.4962626269047, 33.542249807155, 36.9526033996693, 
46.0464770178552, 49.5240246297537, 49.6298842520857, 47.9889200846868, 
40.3862301499032, 36.8605803231892, 38.8799158911488, 39.0120455451407, 
45.9071510330717, 55.8730250709158, 60.4339176493461), SPEI = c(0.385767341805337, 
-0.240467091114443, 0.218601001011986, 0.392296211626228, -0.0041472667529566, 
0.36089672045203, -0.415596363086708, -0.694577131096395, -0.53422184521265, 
0.372791671097943, 0.0714646484375678, 0.100567550879492, 0.484279813014397, 
-0.478876226785371, -0.591222448288627, -0.473201395390211, -0.347352514594038, 
-0.432571106796894, -0.259775061906046, 0.114961224539346)), row.names = c(NA, 
20L), class = "data.frame")

Here is the code:

## Plot first set of data and draw its axis
par(mar = c(5, 5, 4, 4))
#VCI index
plot(variables$Time, variables$VCI, pch=20, cex=.9, axes=FALSE, ylim=c(0,100), xlab="", ylab="", 
     type="l",col="Aquamarine3", main="Temporal trend - drought indices, growing season")
axis(2, ylim=c(0,100),col="black",las=1)  ## las=1 makes horizontal labels
mtext("VCI and TCI",side=2,line=2.5)
box()
abline(h = 40, col = "black", lty = "dotted", lwd= 2)


## Allow a second plot on the same graph
par(new=TRUE)

#TCI index
plot(variables$Time, variables$TCI, pch=21, cex = 1.2, axes=FALSE, ylim=c(0,100), xlab="", ylab="", 
     type="l",col="Chocolate2", main="Temporal trend - drought indices, growing season")
axis(2, ylim=c(0,100),col="grey",las=1)  ## las=1 makes horizontal labels
mtext("VCI and TCI",side=2,line=2.5)
box()


## Allow a second plot on the same graph
par(new=TRUE)

## Plot the second plot and put axis scale on right
plot(variables$Time, variables$SPEI, pch=15, cex=.4, xlab="", ylab="", ylim=c(-2.5,2.5), 
     axes=FALSE, type="l", col="darkorchid3")
abline(h = -1.5, col = "black", lty = "dashed", lwd= 2)

## a little farther out (line=4) to make room for labels
mtext("SPEI",side=4,line=2.5)
axis(4, ylim=c(-2.5, 2.5), col="black",col.axis="black",las=1)

par(new=TRUE)
## Draw the time axis
axis(1, spei.df$Time, format(spei.df$Time, "%Y"), 20)
mtext("Time",side=1,col="black",line=2.5)  

## Add Legend
legend("topleft",legend=c("VCI","TCI", "SPEI"),
       text.col=c("black","black", "black"), lty=1, lwd=2, col=c("Aquamarine3","Chocolate2", "darkorchid3"))

CodePudding user response：

I'm not sure how you do this in base R, but you can align a second axis with a little math and trial and error with ggplot2

library(tidyverse)

dat |>
  mutate(SPEI = (SPEI (25/15))*15)|>
  pivot_longer(VCI:SPEI) |>
  mutate(name = factor(name, levels = c("VCI", "TCI", "SPEI"))) |>
  ggplot(aes(Time, value, color = name)) 
  geom_line() 
  geom_hline(linetype = "dotted", yintercept = 40)  
  geom_hline(linetype = "dashed", yintercept = 2.5) 
  scale_color_manual(values = c("Aquamarine3", "Chocolate2", "darkorchid3")) 
  scale_y_continuous(name = "VCI and TCI", 
                     breaks = seq(0, 100, by = 20),
                     limits = c(0, 100),
                     sec.axis = sec_axis(trans = ~(./15)-(25/15), name = "SPEI", 
                                         breaks = seq(-2,5, by = 1))) 
  ggtitle("Temporal trend - drought indices, growing season") 
  labs(color = "") 
  theme(panel.background = element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.key = element_blank(),
        legend.box.background = element_rect(fill = NA, color = "black"),
        legend.position = c(0.07,0.87),
        legend.title = element_blank())

UPDATE edited to make 40 = -1

library(tidyverse)

dat |>
  mutate(SPEI = ((SPEI (55/15))*15))|>
  pivot_longer(VCI:SPEI) |>
  mutate(name = factor(name, levels = c("VCI", "TCI", "SPEI"))) |>
  ggplot(aes(Time, value, color = name)) 
  geom_line() 
  geom_hline(linetype = "dotted", yintercept = 40)  
  geom_hline(linetype = "dashed", yintercept = 32.5) 
  scale_color_manual(values = c("Aquamarine3", "Chocolate2", "darkorchid3")) 
  scale_y_continuous(name = "VCI and TCI", 
                     breaks = seq(0, 100, by = 20),
                     limits = c(0, 100),
                     sec.axis = sec_axis(trans = ~(./15)-(55/15), name = "SPEI", 
                                         breaks = seq(-4,3, by = 1))) 
  ggtitle("Temporal trend - drought indices, growing season") 
  labs(color = "") 
  theme(panel.background = element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.key = element_blank(),
        legend.box.background = element_rect(fill = NA, color = "black"),
        legend.position = c(0.07,0.87),
        legend.title = element_blank())

CodePudding user response：

The question can be solved by calculating the maximum value of one of the axes from the other using the linear interpolation formula, see Wikipedia for an explanation.

Here, I prefer to give a more compact and generic answer, that can be re-used in similar cases. The solution below uses random data and omits problem-specific labeling.

It works as follows:

Create some random data with one series for x and two series ya and yb for y.
Set the y axis limits for ya, the minimum for yb and the matching values for both axes, but omit the maximum value for the second axis.
Calculate the maximum for yb.
The rest can then be done as usual.

## create a few random test data
set.seed(123)
x <- 1:10
ya <- runif(x, min = 20, max = 70)
yb <- runif(x, min = -1.5, max= 1)

## set limits for the two y axes, but omit yb_max
ya_min <- 0; ya_max <- 100; yb_min <- -3

## set value, where the two axes should match
ya_match <- 40; yb_match <- -1

## calculate yb_max using linear interpolation (here: extrapolation)
yb_max <- yb_min   (ya_max - ya_min) * (yb_match - yb_min) / (ya_match - ya_min)

## set graphical parameter to increase margin for right axis
par(mar=c(5, 4, 5, 5)   .1)

## plot data set "a" 
plot(x, ya, type = "l", col = "black", ylim = c(ya_min, ya_max))
abline(h = ya_match, col = "grey", lty = "dashed")

## add a second plot for data set "b", omitting axes and axis labels
par(new = TRUE)
plot(x, yb, type = "l", col = "red", ylim = c(yb_min, yb_max),
     axes = FALSE, xlab = "", ylab = "")

## add right axis and axis label
axis(4)
mtext("yb", side = 4, line = 3, col = "red")

plot with two matching y axes