Select Git revision
cats_analysis.R
cats_analysis.R 6.17 KiB
library(proj4)
library(terra)
read_stats <- function(stats_file) {
if (file.exists(stats_file)) {
return(utils::read.csv(stats_file, sep=",", comment.char="#"))
}
return (NULL)
}
get_carrying_capacity_maximum <- function(catsconfig, run.name)
{
cc <- NA
config <- catsconfig
if (class(catsconfig) == "character") {
if (file.exists(catsconfig)) {
config <- load_cats_configuration(catsconfig)
} else {
warning(sprintf("Could not load CATS configuration file %s", catsconfig))
}
}
if (class(config) == "data.frame") {
config <- config[which(config$general.run_name == run.name), ]
if (nrow(config) == 1) {
cc <- as.numeric(config$species.carrying_capacity_maximum)
}
}
if (is.na(cc))
{
warning("Could not load carrying capacity maximum from configuration")
}
return(cc)
}
#' Plot the results of a CATS simulation
#'
#' Plots the start and end population distribution of a CATS simulation,
#' the range gains and losses, and the development of the number of populated cells
#' (fit and unfit) over time.
#'
#' @param results the data frame returned by [run_cats()]
#' @param catsconfig the configuration passed to [run_cats()], either a configuration data frame or file name
#' @param run.name the name of the CATS simulation to be plotted. Default NULL. If not specified, the run name will be taken
#' from the configuration file, or form the first row of the configuration data frame.
#' @param replicate the replicate number of the CATS simulation to be plotted. Default 0.
#'
#' @return NULL
#' @export
#'
plot_cats_run <- function(results, catsconfig, run.name = NULL, replicate = 0) {
if (is.null(run.name) && class(catsconfig) == "character") {
run.name <- get_run_name(catsconfig)
} else if (is.null(run.name) && class(catsconfig) == "data.frame") {
run.name = catsconfig$general.run_name[1]
}
if (is.null(run.name)) {
stop("Run name not specified and not available from catsconfig")
}
cc <- get_carrying_capacity_maximum(catsconfig, run.name)
stats <- NA
stats_rows <- results[which(results$file.type == "grid-stats" & results$run.name == run.name), ]
stats_rows <- stats_rows[which(stats_rows$replicate == replicate), ]
if (nrow(stats_rows) == 1) {
stats_fn <- stats_rows$file.name
if (is.null(stats_fn)) {
stop(sprintf("Could not load stats file for run %s and replicate %d", run.name, replicate))
}
stats <- utils::read.csv(stats_fn, comment.char = "#", sep=",")
} else {
stop(sprintf("Could not load stats file for run %s and replicate %d", run.name, replicate))
}
to_plot <- results[which(results$run.name == run.name & results$replicate == replicate & results$file.type=="population"), ]
first_plot <- utils::head(to_plot, 1)
last_plot <- utils::tail(to_plot, 1)
to_plot <- rbind(first_plot, last_plot)
old_par <- graphics::par("mfrow", "mar", "mgp")
stats_file <- results[which(results$run.name == run.name & results$replicate == replicate & results$file.type=="grid-stats"), ]$file.name
graphics::par(mfrow = c(2, 2))
graphics::par(mgp = c(2,1,0))
first <- NULL
last <- NULL
first_year <- NULL
last_year <- NULL
#tryCatch(
# {
for (fn in to_plot$file.name) {
this <- to_plot[which(to_plot$file.name == fn), ]
r <- terra::rast(fn)
title <- sprintf("%s - year %d", run.name, this$year)
if (is.null(first)) {
first <- r
first_year <- this$year
}
last_year <- this$year
if (is.na(cc)) {
terra::plot(r, main = title)
} else {
terra::plot(r, range = c(0, cc), main = title)
}
last <- r
}
if (is.null(first) == F && is.null(last) == F) {
p_colors <- data.frame(value = c(-1, 0, 1), color = c("red", "gray", "green"))
first <- (first > 0)*1
last <- (last > 0)*1
title <- sprintf("%s\nrange gains (+1, green) and losses (-1, red) from %d to %d", run.name, first_year, last_year)
x <- last - first
no_pop <- (last == 0 & first == 0)*1
no_pop <- terra::mask(x * 0, mask=no_pop, maskvalues=0) * 0
losses <-(x == -1)*1
losses[losses != 1] <- NA
gains <- (x == 1)*1
gains[gains != 1] <- NA
same <- (x == 0)* 1
same[same != 1] <- NA
terra::plot(same, col = "gray", main=title, legend=F)
terra::plot(losses , col = "red", add=T, legend=F, axes=F)
terra::plot(gains , col = "green", add=T, legend=F, axes=F)
terra::plot(no_pop, col="lightgray", add=T, legend=F, axes=F)
}
stats <- read_stats(stats_file)
years <- c(0)
populated <- c(0)
populated_fit <- c(0)
populated_unfit <- c(0)
max_cells <- 0
if (! is.null(stats)) {
stats_data <- stats[which(stats$phase == "S"), ]
populated <- stats_data$populated
populated_fit <- stats_data$populated_fit
populated_unfit <- stats_data$populated_unfit
max_cells <- max(populated)
years <- stats_data$year
}
plot(populated ~ years, type="l", ylim=c(0, max_cells), main="Range size (cells)", xlab="Year", ylab="Cells")
graphics::lines(populated_fit ~ years, col="blue")
graphics::lines(populated_unfit ~ years, col="red")
graphics::legend(x = min(years), y=max(populated/2), legend=c("Populated cells (all)", "Populated cells (fit)", "Populated cells (unfit)"), col=c("black", "blue", "red"), lty=1)
#},
#error = function(cond) {
# print(cond)
#},
#finally = function(cond){
# graphics::par(old_par)
#}
#)
return(invisible(NULL))
}