# define output files and see if they have already been generated
out_file_mm <- file.path(samp_dir, "mm", "pca_full.tif")
out_file_nk <- file.path(samp_dir, "nk", "pca_full.tif")
if (!all(file.exists(out_file_mm), file.exists(out_file_nk))){
# fix names
names(samp_rast_mm) <- sub("\\.", "_", names(samp_rast_mm))
names(samp_rast_mm) <- sub("-", "_", names(samp_rast_mm))
names(samp_rast_mm) <- tolower(names(samp_rast_mm))
names(samp_rast_nk) <- sub("\\.", "_", names(samp_rast_nk))
names(samp_rast_nk) <- sub("-", "_", names(samp_rast_nk))
names(samp_rast_nk) <- tolower(names(samp_rast_nk))
# run pca
pca_mm <- prcomp(samp_rast_mm, center = T, scale. = T, maxcell = 1e6)
pca_nk <- prcomp(samp_rast_nk, center = T, scale. = T, maxcell = 1e6)
# get percent variance explained
pve_mm <- pca_mm$sdev^2 / sum(pca_mm$sdev^2)
pve_nk <- pca_nk$sdev^2 / sum(pca_nk$sdev^2)
# get cumulative variance explained
cve_mm <- cumsum(pve_mm)
cve_nk <- cumsum(pve_nk)
# get point at which more than 99% of variance is explained
cutoff_mm <- which(cve_mm < 0.99)
cutoff_nk <- which(cve_nk < 0.99)
# get pca prediction
pca_rast_mm <- predict(samp_rast_mm, pca_mm, index = cutoff_mm)
pca_rast_nk <- predict(samp_rast_nk, pca_nk, index = cutoff_nk)
# write to files
writeRaster(pca_rast_mm, out_file_mm, overwrite = T)
writeRaster(pca_rast_nk, out_file_nk, overwrite = T)
# save pca info
saveRDS(pca_mm, sub("\\.tif", ".rds", out_file_mm))
saveRDS(pca_nk, sub("\\.tif", ".rds", out_file_nk))
} else {
# read the existing files in
pca_rast_mm <- rast(out_file_mm)
pca_rast_nk <- rast(out_file_nk)
pca_mm <- readRDS(sub("\\.tif", ".rds", out_file_mm))
pca_nk <- readRDS(sub("\\.tif", ".rds", out_file_nk))
}
# get percent variance explained
pve_mm <- pca_mm$sdev^2 / sum(pca_mm$sdev^2)
pve_nk <- pca_nk$sdev^2 / sum(pca_nk$sdev^2)
# get cumulative variance explained
cve_mm <- cumsum(pve_mm)
cve_nk <- cumsum(pve_nk)
# get point at which more than 99% of variance is explained
cutoff_mm <- which(cve_mm < 0.99) |> max()
cutoff_nk <- which(cve_nk < 0.99) |> max()
# plot out pca variance explained results -- mm
par(mfrow = c(1,2), las = 1)
plot(cve_mm, pch = 16, xlab = "PC", ylab = "Variance Explained",
main = "Monroe Mountain", col = "gray")
lines(
x = c(par("usr")[1], cutoff_mm),
y = rep(cve_mm[cutoff_mm],2),
col = "red"
)
lines(
x = rep(cutoff_mm,2),
y = c(par("usr")[3], cve_mm[cutoff_mm]),
col = "red"
)
points(
x = cutoff_mm,
y = cve_mm[cutoff_mm],
pch = 16,
col = "red"
)
text(
x = cutoff_mm,
y = mean(par("usr")[3:4]),
pos = 4,
labels = paste0("99% Var Expl\nPC", cutoff_mm),
col = "red"
)
# plot out pca variance explained results -- nk
plot(cve_nk, pch = 16, xlab = "PC", ylab = "Variance Explained",
main = "North Kaibab", col = "gray")
lines(
x = c(par("usr")[1], cutoff_nk),
y = rep(cve_nk[cutoff_nk],2),
col = "red"
)
lines(
x = rep(cutoff_nk,2),
y = c(par("usr")[3], cve_nk[cutoff_nk]),
col = "red"
)
points(
x = cutoff_nk,
y = cve_nk[cutoff_nk],
pch = 16,
col = "red"
)
text(
x = cutoff_nk,
y = mean(par("usr")[3:4]),
pos = 4,
labels = paste0("99% Var Expl\nPC", cutoff_nk),
col = "red"
)
