Copied data folder from persefone-model

README.md

@@ -37,12 +37,12 @@ releases.
 
 ![Persefone.jl Desktop screenshot](screenshot.png)
 
-The main window component is the map view. This displays a land cover map of the
+The main window component is the **map view**. This displays a land cover map of the
 simulated region: dark green are forests, light green grassland, yellow fields, red
 built-up areas and blue water. On it, little circles show the position of individual
 animals, with different species denoted by different colours.
 
-Below the map is the control bar, with the following elements (from left to right):
+Below the map is the **control bar**, with the following elements (from left to right):
 
 - **Back button:** Rewind the simulation by one day.
 
@@ -58,7 +58,7 @@ dates of the simulation.
 
 - **Date:** Shows the simulation date currently displayed on the map.
 
-Above the map is the menu bar, with the following options:
+Above the map is the **menu bar**, with the following options:
 
 **Simulation:**
 

data/crops/almass/convert_almass_data.py

+#!/usr/bin/python3
+###
+### ALMaSS uses a pretty convoluted file format to specify its crop growth parameters.
+### This script converts the data to a more easily understandable long-table CSV format.
+###
+### Daniel Vedder, 02/08/2023
+###
+
+input_file = "almass_crop_growth_curves.pre"
+output_file = "almass_crop_growth_curves.csv"
+
+"""
+Read and parse an ALMaSS crop parameter file, returning it as a long-form array.
+
+Original data format:
+
+- first line: number of curves defined in the file
+
+- curve definition:
+  - first line: curve ID
+	- last figure: crop number
+	- if <100 -> normal growth, if >100 -> low nutrient growth
+  - each curve includes five growth phases
+  - each growth phase consists of a series of 15 inflection points and slopes
+  - four variable columns:
+	- col 1: growing degree-days at inflection point or start marker
+	- col 2: leaf area total slope
+	- col 3: leaf area green slope
+	- col 4: height slope
+  - if col 1 == -1: set variables to values in this row (start of new growth phase)
+  - if col 1 == 99999: no further growth in this growth phase
+
+See `PlantGrowthData::PlantGrowthData()` in `Landscape/plants.cpp` in the ALMaSS code.
+"""
+def parse_almass_file(filename):
+    # read in all lines and split them into their component numbers
+    cf = open(filename, 'r')
+    raw_data = cf.readlines()[1:] # the first line is irrelevant
+    cf.close()
+    for i in range(len(raw_data)):
+        raw_data[i] = list(map(float, raw_data[i].strip().split(sep='\t')))
+    # reformat the data
+    data = []
+    curve_id = 0
+    crop_name = ""
+    nutrient_status = ""
+    phases = ["janfirst", "sow", "marchfirst", "harvest1", "harvest2"]
+    growth_phase = -1
+    inflection_point = 0
+    points_per_phase = 15
+    for l in raw_data:
+        # the beginning of a new curve definition is denoted by a single crop number
+        if len(l) == 1:
+            curve_id = int(l[0])
+            inflection_point = 0
+            growth_phase = -1
+            crop_name, nutrient_status = interpret_crop_number(curve_id)
+            #print("Adding crop "+crop_name+" ("+nutrient_status+" nutrients)")
+            continue
+        # there are 15 inflection points per growth phase
+        if inflection_point % points_per_phase == 0: growth_phase += 1
+        inflection_point += 1
+        if (inflection_point-1) % points_per_phase > 0 and (l == [0,0,0,0] or l == [99999,0,0,0]):
+            continue # we can ignore empty inflection points
+        # save each relevant inflection or start point with their associated values
+        if l[0] == -1: point_type = "start"
+        else: point_type = "slope"
+        #if l[0] == -1 or l[0] == 99999: l[0] = 0 #XXX should we keep sign-post values?
+        data.append([curve_id, crop_name, nutrient_status,
+                     phases[growth_phase], point_type]+l)
+    return data    
+
+"""
+Each growth curve is associated with a crop number, which specifies the crop type
+and nutrient growing conditions.
+
+See `PlantGrowthData::VegTypeToCurveNum()` in `Landscape/plants.cpp` in the ALMaSS code.
+"""
+def interpret_crop_number(i):
+    if i in [1, 101]: crop_type = "spring barley"
+    elif i in [2, 102]: crop_type = "winter barley"
+    elif i in [3, 103]: crop_type = "spring wheat"
+    elif i in [4, 104]: crop_type = "winter wheat"
+    elif i in [5, 105]: crop_type = "winter rye"
+    elif i in [6, 106]: crop_type = "oats"
+    elif i in [7, 107]: crop_type = "triticale"
+    elif i in [8, 108]: crop_type = "maize"
+    elif i in [13, 113]: crop_type = "undersown spring barley"
+    elif i in [21, 121]: crop_type = "spring rape"
+    elif i in [22, 122]: crop_type = "winter rape"
+    elif i == 25: crop_type = "permanent grassland (low yield)"
+    elif i == 26: crop_type = "permanent grassland (grazed)"
+    elif i == 27: crop_type = "permanent grassland (seeded)"
+    #elif i == 28: crop_type = "NA" #XXX unknown number that appears in the data
+    elif i == 29: crop_type = "fodder/clover"
+    elif i == 30: crop_type = "peas/beans"
+    elif i in [41, 141]: crop_type = "carrots"
+    elif i in [50, 150]: crop_type = "potatoes"
+    elif i in [60, 160]: crop_type = "beet"
+    elif i in [65, 165]: crop_type = "silage clover/grass"
+    elif i == 70: crop_type = "lucerne"
+    elif i == 80: crop_type = "tulips"
+    elif i == 90: crop_type = "natural grass"
+    elif i == 91: crop_type = "no growth"
+    elif i == 92: crop_type = "permanent set-aside"
+    elif i == 94: crop_type = "lawn"
+    #elif i == 99: crop_type = "NA" #XXX unknown number that appears in the data
+    elif i == 112: crop_type = "heath"
+    else:
+        crop_type = "NA"
+        print("WARNING: unknown crop type "+str(i))
+    if i > 100: nutrient_status = "low"
+    else: nutrient_status = "high"
+    return crop_type, nutrient_status
+
+"""
+Reformat the crop parameters and output them as a long-form CSV table.
+"""
+def print_csv_file(crop_data, filename):
+    cf = open(filename, 'w')
+    cf.write("curve_id,crop_name,nutrient_status,growth_phase,point_type,GDD,LAI_total,LAI_green,height\n")
+    for entry in crop_data:
+        line = ""
+        for e in entry:
+            line += str(e)+","
+        cf.write(line[:-1]+'\n')
+    cf.close()
+
+
+if __name__ == "__main__":
+    data = parse_almass_file(input_file)
+    print_csv_file(data, output_file)

data/crops/almass/crop_data_general.csv

+name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp
+"winter barley","15 September","30 September",NA,NA,0
+"spring barley","1 March","10 April",NA,NA,0
+"peas/beans",NA,NA,NA,NA,5
+"spring rape",NA,NA,NA,NA,NA
+"winter rape",NA,NA,NA,NA,NA
+"winter rye",NA,NA,NA,NA,NA
+"winter wheat","15 October","31 October",NA,NA,0
+"beet",NA,NA,NA,NA,NA
+"maize",NA,NA,NA,NA,8
+"permanent grassland (grazed)",NA,NA,NA,NA,NA
+"permanent grassland (seeded)",NA,NA,NA,NA,NA
+"fodder/clover",NA,NA,NA,NA,NA
+"natural grass",NA,NA,NA,NA,NA
+"potatoes",NA,NA,NA,NA,4
+"undersown spring barley",NA,NA,NA,NA,0
+"carrots",NA,NA,NA,NA,NA
+"oats",NA,NA,NA,NA,NA
+"permanent set-aside",NA,NA,NA,NA,NA
+"lucerne",NA,NA,NA,NA,NA
+"triticale","25 September","10 October",NA,NA,NA

data/regions/extract_weather_data.R

+#!/usr/bin/Rscript
+###
+### Extract the relevant data from DWD weather files. (See the HTML documentation
+### for instructions on how to obtain the data files.)
+###
+### Daniel Vedder, 27/07/2023
+### 
+
+library(tidyverse)
+
+## replace this with the correct file name
+weatherfile = "daily_weather_jena_dwd/produkt_klima_tag_18210101_20221231_02444.txt"
+
+data = read.table(weatherfile, sep=";", header=T)
+
+weather = data %>%
+    ## drop values before 2000 to save space
+	filter(MESS_DATUM>=20000101) %>%
+    ## select relevant variables and convert place-holder values to NA
+    select(MESS_DATUM, FM, RSK, SDK, VPM, TMK, TXK, TNK) %>%
+    mutate(date=MESS_DATUM, MESS_DATUM=NULL,
+           mean_windspeed=na_if(FM, -999), FM=NULL,
+           precipitation=na_if(RSK, -999), RSK=NULL,
+           sunshine_hours=na_if(SDK, -999), SDK=NULL,
+           mean_vapour_pressure=na_if(VPM, -999), VPM=NULL,
+           mean_temperature=na_if(TMK, -999), TMK=NULL,
+           max_temperature=na_if(TXK, -999), TXK=NULL,
+           min_temperature=na_if(TNK, -999), TNK=NULL)
+
+## replace with the desired file name
+write.csv(weather, file="weather_jena.csv", row.names=FALSE)