Skip to content
Snippets Groups Projects
Select Git revision
  • e0de57c783b26a9ff8bc098da55c2ceb35bba56a
  • main default protected
2 results

cats_analysis.R

Blame
  • 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))
    }