R how to draw ellipses and their axis (major and minor) over a chart with ggplot2

Question

I am trying to reproduce an existing chart with R and ggplot2. This chart has 3 ellipses for tolerance level (50% - 75% - 95%) with a specific angle. What I am able to do is to design the chart (x and y axis), but I am not able to reproduce the ellipses. Here there is what I achieved:

And here there is the code used to make this chart (i have a background image that helps me to know if I am close to the original chart).

ggplot() +  
  annotation_custom(imgonchart, -Inf, Inf, -Inf, Inf) +
  geom_point(data=df, aes(rz,xc)) +
  geom_point(data=biadf, aes(mean(Rzm), mean(Xcm))) +
  geom_ellipse(aes(x0 = 298.6, y0 = 30.8, a = 43, b = 8.9, angle = pi / 39.30 ), color="green") +
  coord_fixed(ratio=10)+
  scale_x_continuous(expand = c(0, 0), limits = c(0,600), breaks = seq(0, 600, by= 50)) +  # set limit, steps, and how much can expand
  scale_y_continuous(expand = c(0, 0), limits = c(0,60), breaks = seq(0, 60, 5)) +
  labs(x = expression(paste("Rz/H (", Omega, "/m)")), 
       y = expression(paste("Xc/H (", Omega, "/m)")),
       ) + # LEGENDA
  theme(plot.title = element_text(hjust = 0.5),
        panel.grid.minor.y =  element_blank(),
        panel.grid.minor.x = element_blank()) # to adjust the title in the middle

As you can see, the green ellipse (the more vivid one) is the one made by me. But when I try to rotate that ellipse, it goes out of the chart. There is a scientific papers that describe how to design the above graph, here the text:

Using the scale of Figure 1 for R/H and Xc/H, draw ellipses with major and minor axes slopes of 69.30° and -20.70° respectively. The semiaxes lenghts of the male panel ellipses are 89 and 43 Ohm/m for the 50% tolerance, 127 and 61 Ohm for the 75% tolerance and 187 and 89 Ohm/m for the 95% tolerance ellipses.

I was also trying to do it using the a "fake" sample like described in the paper with the rnorm() and the stat_ellipse(). But I was not able to reach the desired result.

Any help will be appreciated. Marco

Answer 1

I think the easiest way to do this would be to stick to the units used in the paper, and fake the axis instead by putting labels at 1/10th the value.

Starting with some dummy data for the df in your code, and converting the angle stated in the paper to radians, we have:

library(ggplot2)
library(ggforce)

df        <- data.frame(rz = c(250, 300, 360), xc = c(35.1, 30.8, 34.8))
angle     <- -20.7
rad_angle <- angle / 180 * pi

And the plot would be:

ggplot() +  
  coord_equal() +
  geom_point(data = df, aes(rz, xc * 10)) +
  geom_ellipse(aes(x0 = 298.6, y0 = 308, a = 5,  b = 10,  angle = rad_angle)) +
  geom_ellipse(aes(x0 = 298.6, y0 = 308, a = 43, b = 89,  angle = rad_angle),
               color = 'green4') +
  geom_ellipse(aes(x0 = 298.6, y0 = 308, a = 61, b = 127, angle = rad_angle),
               color = 'orange') +
  geom_ellipse(aes(x0 = 298.6, y0 = 308, a = 89, b = 187, angle = rad_angle),
               color = 'red4') +
  geom_segment(aes(x = 298.6 - 89 * cos(rad_angle), 
                   xend = 298.6 + 89 * cos(rad_angle),
                   y = 308 - 89 * sin(rad_angle),
                   yend = 308 + 89 * sin(rad_angle)), linetype = 2, 
               color = 'gray50') +
  geom_segment(aes(x = 298.6 + 187 * sin(rad_angle), 
                   xend = 298.6 - 187 * sin(rad_angle),
                   y = 308 - 187 * cos(rad_angle),
                   yend = 308 + 187 * cos(rad_angle)), linetype = 2, 
               color = 'gray50') +
  scale_x_continuous(expand = c(0, 0), limits = c(0,600), 
                     breaks = seq(0, 600, by = 50)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0,600), 
                     breaks = seq(0, 600, 50), labels = ~.x/10) +
  labs(x = expression(paste("Rz/H (", Omega, "/m)")), 
       y = expression(paste("Xc/H (", Omega, "/m)"))) +
  theme(plot.title = element_text(hjust = 0.5),
        panel.grid.minor.y =  element_blank(),
        panel.grid.minor.x = element_blank())