R - update boxplot axis range after adding points

I have a boxplot which summarizes ~60000 turbidity data points into quartiles, median, whiskers and sometimes outliers. Often a few outliers are so high up that the whole plot is compressed at the bottom, and I therefor choose to omit the outliers. However, I also have added averages to the plots as points, and I want these to be plotted always. The problem is that the y-axis of the boxplot does not adjust to the added average points, so when averages are far above the box they are simply plotted outside the chart window (see X-point for 2020, but none for 2021 or 2022). Normally with this parameter, the average will be between the whisker end and the most extreme outliers. This is normal, and expected in the data.

I have tried to capture the boxplot y-axis range to compare with the average, and then setting the ylim if needed, but I just don't know how to retrieve these axis ranges.

My code is just
boxplot(...)
points(...)
and works as far as plotting the points. Just not adjusting the y-axis.

Question 1: is it not possible to get the boxplot to redraw with the new points data? I thought this was standard in R plots.
Question 2: if not, how can I dynamically adjust the y-axis range?

Let's try to show a concrete example of the problem with some simulated data:

set.seed(1)

df <- data.frame(y = c(rexp(99), 150), x = rep(c("A", "B"), each = 50))

Here, group "B" has a single outlier at 150, even though most values are a couple of orders of magnitude lower. That means that if we try to draw a boxplot, the boxes get squished at the bottom of the plot:

boxplot(y ~ x, data = df, col = "lightblue")

If we remove outliers, the boxes plot nicely:

boxplot(y ~ x, data = df, col = "lightblue", outline = FALSE)

The problem comes when we want to add a point indicating the mean value for each boxplot, since the mean of "B" lies outside the plot limits. Let's calculate and plot the means:

mean_vals <- sapply(split(df$y, df$x), mean)
mean_vals
#>         A         B 
#> 0.9840417 4.0703334

boxplot(y ~ x, data = df, col = "lightblue", outline = FALSE)
points(1:2, mean_vals, cex = 2, pch = 16, col = "red")

The mean for "B" is missing because it lies above the upper range of the plot.

The secret here is to use boxplot.stats to get the limits of the whiskers. By concatenating our vector of means to this vector of stats and getting its range, we can set our plot limits exactly where they need to be:

y_limits <- range(c(boxplot.stats(df$y)$stats, mean_vals))

Now we apply these limits to a new boxplot and draw it with the points:

boxplot(y ~ x, data = df, outline = FALSE, ylim = y_limits, col = "lightblue")
points(1:2, mean_vals, cex = 2, pch = 16, col = "red")

For comparison, you could do the whole thing in ggplot like this:

library(ggplot2)

ggplot(df, aes(x, y)) +
  geom_boxplot(fill = "lightblue", outlier.shape = NA) +
  geom_point(size = 3, color = "red", stat = "summary", fun = mean) +
  coord_cartesian(ylim = range(c(range(c(boxplot.stats(df$y)$stats, 
                                         mean_vals))))) +
  theme_classic(base_size = 16)

^{Created on 2023-02-05 with reprex v2.0.2}