#!/usr/bin/Rscript
### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file visualises the output of the nature model.
###

##TODO replace this with Julia code using Makie (issue #47)

library(tidyverse)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(terra) #https://rspatial.github.io/terra/index.html

autorun = TRUE ## automatically run all analyses?

datadir = "results" ## "results" by default, can be specified via commandline
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.csv(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.csv(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) {
    ## the last commandline argument is interpreted as the output directory
    ##TODO make map files configurable
    arg = commandArgs()[length(commandArgs())]
    if (dir.exists(arg)) datadir = arg
    populationTrends()
    visualiseMap()
}