Principal components analysis of South-east Australian environment

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf"))
library(eucs)

vars <- c(
  "Relief" = "relief_1000m_radius", "TWI" = "topographic_wetness_index_3sec",
  "Moisture" = "mean_moisture_index_of_low_qtr", "Nitrogen" = "total_nitrogen",
  "Radiation" = "annual_mean_radiation",
  "Precipitation" = "annual_precipitation",
  "Temperature" = "annual_mean_temperature",
  "Phosphorus" = "total_phosphorus", "Sand" = "sand",
  "Soil Depth" = "depth_of_soil"
)
rgn_names <- sort(c("Greater Grampians", unique(southeast$ibra_subregion)))
pca_data <- na.omit(rbind(covariates_grampians, covariates_southeast))
pca_data <- subset(pca_data, ibra_subregion %in% rgn_names)
pca <- prcomp(pca_data[, vars], center = TRUE, scale. = TRUE)

rgns <- factor(pca_data$ibra_subregion)
n <- NROW(pca$x)
lam <- pca$sdev[1:2] * sqrt(n)
y <- t(t(pca$rotation[, 1:2]) * lam)
x <- t(t(pca$x[, 1:2]) / lam)
n <- nrow(x)
p <- nrow(y)
unsigned.range <-
  function(x) c(-abs(min(x, na.rm = TRUE)), abs(max(x, na.rm = TRUE)))
rangx1 <- unsigned.range(x[, 1])
rangx2 <- unsigned.range(x[, 2])
rangy1 <- unsigned.range(y[, 1])
rangy2 <- unsigned.range(y[, 2])
xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2)
ratio <- max(rangy1 / rangx1, rangy2 / rangx2) / .9

colors <- rainbow(nlevels(rgns), alpha = .5)
names(colors) <- levels(rgns)
colors[
  c(
    "Bateman", "Bondo", "East Gippsland Lowlands", "Ettrema", "Gippsland Plain",
    "Highlands-Northern Fall", "Highlands-Southern Fall", "Kybeyan-Gourock",
    "Monaro", "Murrumbateman", "South East Coastal Ranges", "Strzelecki Ranges",
    "Victorian Alps"
  )
] <- grey.colors(1, alpha = .1)
par(mar = c(5, 4, 8, 3))
plot(
  x, xlim = xlim, ylim = ylim, cex = .25,
  col = colors[rgns], pch = 19,
  panel.first = grid()
)
par(xpd = TRUE)
legend(
  "topleft", inset = c(0, -.21), legend = levels(rgns),
  fill = colors, ncol = 4, bty = "n", x.intersp = .8, cex = .7
)
par(new = TRUE)
dev.hold()

plot(
  y, axes = FALSE, type = "n", xlim = xlim * ratio, ylim = ylim * ratio,
  xlab = "", ylab = ""
)
box()
text(y, labels = names(vars), cex = .6)
arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, length = 0.1)

dev.flush()