How to make 2 axes cross at 0 in RStudio all the time (even when resizing the plot device)

Question

I want to be able to plot 2 axes on a graph so the they cross at 0. The question was asked about it here and solved for the crossing axes part. I can't find a way to make the two axes resize in an aspect ratio in RStudio. The graph made by biplot.prcomp is actually doing this, but I'm not able to find the code that enables the resize-axes-in-a-certain-ratio thing.

new_lim <- function(a, type = 1) {
  newdata_ratio <-  NULL
  i <- type * 2 - 1
  old_lim <- par("usr")[i:(i+1)] + c(diff(par("usr")[i:(i+1)]) * 0.04 / 1.08, 
                                     diff(par("usr")[i:(i+1)]) * -0.04 / 1.08)
  old_ratio <- old_lim[1] / old_lim[2]
  newdata_ratio <- if (max(a) <= 0) -1.0e+6 else min(a) / max(a)
  if (old_ratio >= newdata_ratio ) {
    new_min <- min(a)
    new_max <- min(a) / old_ratio
  } else {
    new_min <- max(a) * old_ratio
    new_max <- max(a)
  }
  c(new_min, new_max)
}

x2 <- -10:20
y2 <- seq(40, 10, length.out = length(x2))
library(vegan)                               
### Fixed elements are prepared.
par(mfrow=c(2,3))
for (i in 1:6) {


### pettern 1 / plot(pca); par(new=T); plot(~, new_lim())
s1= rnorm(50,mean = 12); s2= rnorm(50, mean = 17); s3= rnorm(50, mean = 20)
pca=rda(cbind(s1,s2,s3))
pca.scoop=scores(pca, scaling = 2)
plot(pca, xlab = "x1", ylab = "y1",
     type = c("p"),
     main= "main",
     scaling = 2,
     choices = c(1,2),
     xlim =c(min(pca.scoop$sites[,1]),max(pca.scoop$sites[,1])),
     ylim = c(min(pca.scoop$sites[,2]),max(pca.scoop$sites[,2])),
     bty = "o",#"l"
     pch=4)
abline(v = 0, lty = 2); abline(h = 0, lty = 2)
par(new =TRUE)
plot(x2, y2, 
     xlim = new_lim(x2), 
     ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red")
axis(4, col = "red", col.axis = "red")
mtext("x2", side = 3, line = 2.5, col = "red")
mtext("y2", side = 4, line = 2.5, col = "red")

## pattern 2 / biplot(pca); par(new=T); plot(~, new_lim())
s1= rnorm(50,mean = 12)
s2= rnorm(50, mean = 17)
s3= rnorm(50, mean = 20)
pca=rda(cbind(s1,s2,s3))
pca.scoop=scores(pca, scaling = 2)
par(new =TRUE)
plot(x2, y2, 
     xlim = new_lim(x2), 
     ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red")
axis(4, col = "red", col.axis = "red")
mtext("x2", side = 3, line = 2.5, col = "red")
mtext("y2", side = 4, line = 2.5, col = "red")
}

Here is the exact same graph, but when I resize the graph manually, the top red x axis is shifted.

Here is an example with prcomp:

for (i in 1:6) {
s1= rnorm(50,mean = 12); s2= rnorm(50, mean = 17); s3= rnorm(50, mean = 20)
pr=prcomp(cbind(s1,s2,s3))
biplot(pr)
abline(h=0,v=0, lty =3)
par(new =TRUE)
plot(x2, y2, axes = F, ann = F)
}

The rescaled image looks the same!:

Edit

This code, when I want to base the second graph on the PC vectors, it's not aligning. Is there a way to align everything and conserve asp = 1?

par(mfrow=c(2,3))
for (i in 1:6) {

  ### pettern 1 / plot(pca); par(new=T); plot(~, new_lim())
  s1= rnorm(50,mean = 12); s2= rnorm(50, mean = 17); s3= rnorm(50, mean = 20)
  pca=rda(cbind(s1,s2,s3))
  pca.scoop=scores(pca, scaling = 2)
  myccaplot(    ### (4)
    pca, xlab = "x1", ylab = "y1",
    type = c("p"),
    main= "main",
    scaling = 2,
    choices = c(1,2),
    xlim =c(min(pca.scoop$sites[,1]),max(pca.scoop$sites[,1])),
    ylim = c(min(pca.scoop$sites[,2]),max(pca.scoop$sites[,2])),
    bty = "o",pch=4,
    axes = F) ### (1)
  axis(1)  ### (2)
  axis(2)  ### (2)
  box()    ### (2)
  x2 = c(min(pca.scoop$species[,1]),max(pca.scoop$species[,1]))
  y2 = c(min(pca.scoop$species[,2]),max(pca.scoop$species[,2]))
  par(new =TRUE)
  plot(x2, y2, 
       xlim = new_lim(x2), 
       ylim = new_lim(y2, 2), axes = F, ann = F)
  axis(3, col = "red", col.axis = "red")
  axis(4, col = "red", col.axis = "red")
  abline(v = 0, lty = 2)      ### (3)
  abline(h = 0, lty = 2)      ### (3)
  mtext("x2", side = 3, line = 2.5, col = "red")
  mtext("y2", side = 4, line = 2.5, col = "red")
}

Answer 1

The solution is a little bit tricky:

1. First, don't draw the axis with plot.cca (the first call of plot of object pca).
2. You can plot the axes yourself afterwards like you do it in the second call of plot.
3. Draw the ablines after the second plot.
4. When testing you will see that there are two types of ablines now, one from plot.cca and one we plotted in 3. When looking into the function plot.cca (which is a S3 method for plot), e.g. with getAnywhere(plot.cca) we see that abline is called twice and there is no option to prevent that. So, a little bit tricky: we define our own plot function and remove the two lines with abline within plot.cca. We can call that function e.g. myccaplot. Now it should work like expected.

In code snippets:

for (i in 1:6) {

### pettern 1 / plot(pca); par(new=T); plot(~, new_lim())
s1= rnorm(50,mean = 12); s2= rnorm(50, mean = 17); s3= rnorm(50, mean = 20)
pca=rda(cbind(s1,s2,s3))
pca.scoop=scores(pca, scaling = 2)
myccaplot(    ### (4)
     pca, xlab = "x1", ylab = "y1",
     type = c("p"),main= "main",scaling = 2,choices = c(1,2),
     xlim =c(min(pca.scoop$sites[,1]),max(pca.scoop$sites[,1])),
     ylim = c(min(pca.scoop$sites[,2]),max(pca.scoop$sites[,2])),
     bty = "o",pch=4,
     axes = F) ### (1)
axis(1)  ### (2)
axis(2)  ### (2)
box()    ### (2)

par(new =TRUE)
plot(x2, y2, 
   xlim = new_lim(x2), 
   ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red")
axis(4, col = "red", col.axis = "red")
abline(v = 0, lty = 2)      ### (3)
abline(h = 0, lty = 2)      ### (3)
mtext("x2", side = 3, line = 2.5, col = "red")
mtext("y2", side = 4, line = 2.5, col = "red")
}

### (4) ## copy modified content of getAnywhere(plot.cca) here
    myccaplot <- function (x, choices = c(1, 2), display = c("sp", "wa", "cn"), 
          scaling = "species", type, xlim, ylim, const, correlation = FALSE, 
          hill = FALSE, ...) 
{
  TYPES <- c("text", "points", "none")
  g <- scores(x, choices, display, scaling, const, correlation = correlation, 
              hill = hill)
  if (length(g) == 0 || all(is.na(g))) 
    stop("nothing to plot: requested scores do not exist")
  if (!is.list(g)) 
    g <- list(default = g)
  for (i in seq_along(g)) {
    if (length(dim(g[[i]])) > 1) 
      rownames(g[[i]]) <- rownames(g[[i]], do.NULL = FALSE, 
                                   prefix = substr(names(g)[i], 1, 3))
  }
  if (!is.null(g$centroids)) {
    if (is.null(g$biplot)) 
      g$biplot <- scores(x, choices, "bp", scaling)
    if (!is.na(g$centroids)[1]) {
      bipnam <- rownames(g$biplot)
      cntnam <- rownames(g$centroids)
      g$biplot <- g$biplot[!(bipnam %in% cntnam), , drop = FALSE]
      if (nrow(g$biplot) == 0) 
        g$biplot <- NULL
    }
  }
  if (missing(type)) {
    nitlimit <- 80
    nit <- max(nrow(g$spe), nrow(g$sit), nrow(g$con), nrow(g$def))
    if (nit > nitlimit) 
      type <- "points"
    else type <- "text"
  }
  else type <- match.arg(type, TYPES)
  if (length(choices) == 1) {
    if (length(g) == 1) 
      pl <- linestack(g[[1]], ...)
    else {
      hasSpec <- names(g)[1] == "species"
      ylim <- range(c(g[[1]], g[[2]]), na.rm = TRUE)
      pl <- linestack(g[[1]], ylim = ylim, side = ifelse(hasSpec, 
                                                         "left", "right"), ...)
      linestack(g[[2]], ylim = ylim, side = ifelse(hasSpec, 
                                                   "right", "left"), add = TRUE, ...)
    }
    return(invisible(pl))
  }
  if (missing(xlim)) {
    xlim <- range(g$species[, 1], g$sites[, 1], g$constraints[, 
                                                              1], g$biplot[, 1], if (length(g$centroids) > 0 && 
                                                                                     is.na(g$centroids)) NA else g$centroids[, 1], g$default[, 
                                                                                                                                             1], na.rm = TRUE)
  }
  if (!any(is.finite(xlim))) 
    stop("no finite scores to plot")
  if (missing(ylim)) {
    ylim <- range(g$species[, 2], g$sites[, 2], g$constraints[, 
                                                              2], g$biplot[, 2], if (length(g$centroids) > 0 && 
                                                                                     is.na(g$centroids)) NA else g$centroids[, 2], g$default[, 
                                                                                                                                             2], na.rm = TRUE)
  }
  plot(g[[1]], xlim = xlim, ylim = ylim, type = "n", asp = 1, 
       ...)
  # abline(h = 0, lty = 3) # removed
  # abline(v = 0, lty = 3) # removed
  if (!is.null(g$species)) {
    if (type == "text") 
      text(g$species, rownames(g$species), col = "red", 
           cex = 0.7)
    else if (type == "points") 
      points(g$species, pch = "+", col = "red", cex = 0.7)
  }
  if (!is.null(g$sites)) {
    if (type == "text") 
      text(g$sites, rownames(g$sites), cex = 0.7)
    else if (type == "points") 
      points(g$sites, pch = 1, cex = 0.7)
  }
  if (!is.null(g$constraints)) {
    if (type == "text") 
      text(g$constraints, rownames(g$constraints), cex = 0.7, 
           col = "darkgreen")
    else if (type == "points") 
      points(g$constraints, pch = 2, cex = 0.7, col = "darkgreen")
  }
  if (!is.null(g$biplot) && nrow(g$biplot) > 0 && type != "none") {
    if (length(display) > 1) {
      mul <- ordiArrowMul(g$biplot)
    }
    else mul <- 1
    attr(g$biplot, "arrow.mul") <- mul
    arrows(0, 0, mul * g$biplot[, 1], mul * g$biplot[, 2], 
           length = 0.05, col = "blue")
    biplabs <- ordiArrowTextXY(mul * g$biplot, rownames(g$biplot))
    text(biplabs, rownames(g$biplot), col = "blue")
    axis(3, at = c(-mul, 0, mul), labels = rep("", 3), col = "blue")
    axis(4, at = c(-mul, 0, mul), labels = c(-1, 0, 1), col = "blue")
  }
  if (!is.null(g$centroids) && !is.na(g$centroids) && type != 
      "none") {
    if (type == "text") 
      text(g$centroids, rownames(g$centroids), col = "blue")
    else if (type == "points") 
      points(g$centroids, pch = "x", col = "blue")
  }
  if (!is.null(g$default) && type != "none") {
    if (type == "text") 
      text(g$default, rownames(g$default), cex = 0.7)
    else if (type == "points") 
      points(g$default, pch = 1, cex = 0.7)
  }
  class(g) <- "ordiplot"
  invisible(g)
}