Wrote analysis script to plot population trends and maps

src/analysis/analyse_nature.R

+#!/usr/bin/Rscript
+### Persephone - a socio-economic-ecological model of European agricultural landscapes.
+###
+### This file visualises the output of the nature model.
+###
+
+library(tidyverse)
+library(ggplot2)
+library(ggsci)
+library(RColorBrewer)
+library(terra)
+
+autorun = TRUE ## automatically run all analyses?
+
+datadir = "results" ## "results" by default
+popfile = "populations.csv"
+indfile = "individuals.csv"
+mapfile = "landcover_jena.tif"
+
+population_output_file = "population_trends"
+map_output_file = "landscape_map"
+
+populationTrends = function() {
+    print("Plotting population trends over time.")
+    popdata = read.csv2(paste(datadir, popfile, sep="/")) %>%
+        mutate(Date = as.POSIXct(strptime(Date,format="%Y-%m-%d")))
+    ggplot(data=popdata, aes(x=Date, y=Abundance, color=Species)) +
+        geom_point() +
+        geom_smooth() +
+        scale_color_viridis_d() +
+        theme_bw()
+    ggsave(paste0(datadir, "/", population_output_file, ".pdf"), width=8, height=4)
+}
+
+visualiseMap = function() {
+    print("Visualising individuals on the landscape map.")
+    landcover = rast(paste(datadir, mapfile, sep="/"))
+    inddata = read.csv2(paste(datadir, indfile, sep="/")) %>% select(Date,Species,X,Y) %>%
+        mutate(Date = as.POSIXct(strptime(Date,format="%Y-%m-%d")))
+    for (d in unique(inddata$Date)) {
+        ## somehow, d is changed into a number by the for loop, so we have to convert back
+        d = format(as.POSIXct(d,origin="1970-01-01"))
+        ##colors = brewer.pal(n=length(unique(inddata$species)), name="Dark2")
+        colors = c("blue", "red") #XXX replace with line above after testing
+        pdf(paste0(datadir, "/", map_output_file, "_", format(d, format="%Y-%m-%d"), ".pdf"))
+        plot(landcover, axes=FALSE)
+        for (s in unique(inddata$Species)) {
+            species_points = inddata %>% filter(Species==s) %>% select(X,Y) %>%
+                mutate(X = xFromCol(landcover, X), Y = yFromRow(landcover, Y)) %>%
+                as.data.frame() %>% vect(geom=c("X","Y"), keepgeom=TRUE)
+            points(species_points, col=colors[1])
+            colors = colors[-1]
+        }
+        dev.off()
+    }
+    ##TODO change plotting to ggplot (https://dieghernan.github.io/tidyterra/)    
+}
+
+## If autorun is set, run the experiment specified via commandline argument
+if (autorun) {
+    populationTrends()
+    visualiseMap()
+    ##TODO set variables from commandline arguments
+    ## arg = commandArgs()[length(commandArgs())]
+    ## if (arg %in% c("tolerance", "habitat", "mutation", "linkage","tolerance_long")) {
+    ##     experiment = arg
+    ##     results = loadData(experiment)
+    ##     plotAll(results)
+    ## }
+    ## else { print(paste("Unknown experiment" , arg)) }
+}

src/nature/wolpertinger.jl

@@ -35,7 +35,7 @@ function updatewolpertinger!(w::Animal, model::AgentBasedModel)
     end
     w.energy -= speed
     # reproduce every once in a blue moon
-    if rand() < 0.05
+    if rand() < trait(w, "fecundity")
         @debug "Wolpertinger $(w.id) has reproduced."
         add_agent!(w.pos, Animal, model, getspecies("Wolpertinger"),
                    hermaphrodite, 0, trait(w, "birthenergy"))
@@ -47,4 +47,5 @@ newspecies("Wolpertinger",
            updatewolpertinger!,
            Dict("popdensity"=>1/10000,
                 "birthenergy"=>400,
+                "fecundity"=>0.02,
                 "maxspeed"=>5))

src/nature/wyvern.jl

@@ -26,6 +26,11 @@ Wyverns are ferocious hunters, scouring the landscape for their favourite
 prey: wolpertingers...
 """
 function updatewyvern!(w::Animal, model::AgentBasedModel)
+    # check if the wyvern has reached its end-of-life
+    if w.age == trait(w, "maxage")
+        kill_agent!(w, model)
+        return
+    end
     # check if a wolpertinger is in pouncing distance
     for a in nearby_agents(w, model, trait(w, "speed"))
         (a.species.name != "Wolpertinger") && continue
@@ -52,7 +57,7 @@ function updatewyvern!(w::Animal, model::AgentBasedModel)
     end
     # reproduce every once in a blue moon
     @label reproduce
-    if rand() < 0.01
+    if rand() < trait(w, "fecundity")
         @debug "Wyvern $(w.id) has reproduced."
         add_agent!(w.pos, Animal, model, getspecies("Wyvern"),
                    hermaphrodite, 0, trait(w, "birthenergy"))
@@ -64,6 +69,8 @@ newspecies("Wyvern",
            updatewyvern!,
            Dict("initpop"=>100,
                 "birthenergy"=>1000,
+                "fecundity"=>0.01,
+                "maxage"=>365,
                 "speed"=>10,
                 "vision"=>50,
                 "pounceenergy"=>20,