Merge branch 'marco/weatherdata' into 'master'

Change weather data representation

See merge request persefone/persefone-model!4

CHANGELOG.md

@@ -24,6 +24,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - added soil map section to GIS docs
 
+- the `Weather` type now stores all weather information for the whole
+  simulation as a "struct-of-arrays" and with a function interface,
+  e.g. `sunshine(weather, date)` (previously, `Weather` was a struct
+  that stored the weather data for one day only, and
+  `AgricultureModel` had a field `weather::Dict{Date,Weather}`, making
+  it a "dict-of-struct" data layout)
+
+- `AgricultureModel` and `Weather` are now defined with `@kwdef`,
+  allowing for keyword arguments in their constructors
+
+- when reading weather data, we now throw an error when there are
+  missing values for the fields `min_temperature`, `max_temperature`,
+  `mean_temperature`, `precipitation`, and
+  `potential_evapotranspiration` (in the future missing values could
+  also be imputed)
+
 ### Deprecated
 
 ### Removed

src/core/simulation.jl

@@ -12,20 +12,20 @@ This is the heart of the model - a struct that holds all data and state
 for one simulation run. It is created by [`initialise`](@ref) and passed
 as input to most model functions.
 """
-mutable struct AgricultureModel{Tcroptype,Tcropstate} <: SimulationModel
-    settings::Dict{String,Any}
+@kwdef mutable struct AgricultureModel{Tcroptype, Tcropstate} <: SimulationModel
+    settings = Dict{String, Any}()
     rng::AbstractRNG
     logger::AbstractLogger
-    dataoutputs::Dict{String,DataOutput}
+    dataoutputs = Dict{String, DataOutput}()
     date::Date
-    landscape::Matrix{Pixel}
-    weather::Dict{Date,Weather}
-    crops::Dict{String,Tcroptype}
-    farmers::Vector{Farmer}
-    farmplots::Vector{FarmPlot{Tcropstate}}
-    animals::Vector{Union{Animal,Nothing}}
-    migrants::Vector{Pair{Animal,AnnualDate}}
-    events::Vector{FarmEvent}
+    landscape = Matrix{Pixel}()
+    weather::Weather
+    crops = Dict{String, Tcroptype}()
+    farmers = Vector{Farmer}()
+    farmplots = Vector{FarmPlot{Tcropstate}}()
+    animals = Vector{Union{Animal, Nothing}}()
+    migrants = Vector{Pair{Animal, AnnualDate}}()
+    events = Vector{FarmEvent}()
 end
 
 """
@@ -36,7 +36,31 @@ Return the total number of agents in a model object.
 function nagents(model::AgricultureModel)
     length(model.animals)+length(model.farmers)+length(model.farmplots)
 end
-    
+
+windspeed(model::AgricultureModel, date::Date) = windspeed(model.weather, date)
+windspeed(model::AgricultureModel) = windspeed(model, model.date)
+
+precipitation(model::AgricultureModel, date::Date) = precipitation(model.weather, date)
+precipitation(model::AgricultureModel) = precipitation(model, model.date)
+
+sunshine(model::AgricultureModel, date::Date) = sunshine(model.weather, date)
+sunshine(model::AgricultureModel) = sunshine(model, model.date)
+
+humidity(model::AgricultureModel, date::Date) = humidity(model.weather, date)
+humidity(model::AgricultureModel) = humidity(model, model.date)
+
+meantemp(model::AgricultureModel, date::Date) = meantemp(model.weather, date)
+meantemp(model::AgricultureModel) = meantemp(model, model.date)
+
+maxtemp(model::AgricultureModel, date::Date) = maxtemp(model.weather, date)
+maxtemp(model::AgricultureModel) = maxtemp(model, model.date)
+
+mintemp(model::AgricultureModel, date::Date) = mintemp(model.weather, date)
+mintemp(model::AgricultureModel) = mintemp(model, model.date)
+
+evapotranspiration(model::AgricultureModel, date::Date) = evapotranspiration(model.weather, date)
+evapotranspiration(model::AgricultureModel) = evapotranspiration(model, model.date)
+
 """
     stepagent!(agent, model)
 
@@ -117,22 +141,14 @@ function initmodel(settings::Dict{String, Any})
                               settings["core.enddate"])
         crops, Tcroptype, Tcropstate = initcropmodel(settings["crop.cropmodel"],
                                                      settings["crop.cropdirectory"])
-        farmers = Vector{Farmer}()
-        farmplots = Vector{FarmPlot{Tcropstate}}()
-        model = AgricultureModel{Tcroptype,Tcropstate}(
+        model = AgricultureModel{Tcroptype, Tcropstate}(;
             settings,
-            StableRNG(settings["core.seed"]),
+            rng = StableRNG(settings["core.seed"]),
             logger,
-            Dict{String,DataOutput}(),
-            settings["core.startdate"],
+            date = settings["core.startdate"],
             landscape,
             weather,
-            crops,
-            farmers,
-            farmplots,
-            Vector{Union{Animal,Nothing}}(),
-            Vector{Pair{Animal, Date}}(),
-            Vector{FarmEvent}()
+            crops
         )
         saveinputfiles(model)
 

src/world/weather.jl

@@ -7,18 +7,81 @@
 """
     Weather
 
-A single weather datum, combining the observations from one day.
-"""
-struct Weather
-    windspeed::Union{Missing,Float64}
-    precipitation::Union{Missing,Float64}
-    sunshine::Union{Missing,Float64}
-    cloudcover::Union{Missing,Float64}
-    humidity::Union{Missing,Float64}
-    meantemp::Union{Missing,Float64}
-    maxtemp::Union{Missing,Float64}
-    mintemp::Union{Missing,Float64}
-    evapotranspiration::Union{Missing,Float64}
+Holds the weather information for the whole simulation period.
+"""
+@kwdef struct Weather
+    firstdate::Date
+    lastdate::Date
+    # TODO: types for some of these fields can be just `Vector{Float64}` after
+    #       fixing missing measurement values (replacing with values
+    #       approximated from previous/later measurements)
+    windspeed::Vector{Union{Missing,Float64}}
+    precipitation::Vector{Float64}
+    sunshine::Vector{Union{Missing,Float64}}
+    cloudcover::Vector{Union{Missing,Float64}}
+    humidity::Vector{Union{Missing,Float64}}
+    meantemp::Vector{Float64}
+    maxtemp::Vector{Float64}
+    mintemp::Vector{Float64}
+    evapotranspiration::Vector{Float64}
+end
+
+function Base.length(weather::Weather)
+    return Dates.value(weather.lastdate - weather.firstdate) + 1
+end
+
+"""
+    daynumber(weather, date)
+
+Returns the number of days, counting `weather.firstdate` as day `1`.
+"""
+function daynumber(weather::Weather, date::Date)
+    @assert weather.firstdate <= date <= weather.lastdate
+    return Dates.value(date - weather.firstdate) + 1
+end
+
+"""
+    findspans(predicate_fn, array) -> Vector{Tuple{Int, Int}}
+
+Returns spans of indices in a 1-d `array` where a `predicate_fn`
+returns `true`.  The spans are returned as a `Vector{Tuple{Int,
+Int}}`, where each tuple is of the form `(start_index, end_index)`.
+"""
+function findspans(predicate_fn, arr::AbstractVector)
+    spans = Tuple{Int, Int}[]
+    startidx = nothing
+    for i = firstindex(arr):lastindex(arr)
+        if predicate_fn(arr[i])
+            if startidx === nothing
+                startidx = i
+            end
+        elseif startidx !== nothing
+            push!(spans, (startidx, i - 1))
+            startidx = nothing
+        end
+    end
+    if startidx !== nothing
+        push!(spans, (startidx, lastindex(arr)))
+    end
+    return spans
+end
+
+"""
+    check_missing_weatherdata(dataframe)
+
+Check the weather input data for missing values in columns where input
+values are required.
+"""
+function check_missing_weatherdata(df::DataFrame)
+    colnames = ["min_temperature", "max_temperature", "mean_temperature",
+                "precipitation", "potential_evapotranspiration"]
+    for colname in colnames
+        spans = findspans(ismissing, getproperty(df, colname))
+        # TODO: fix missing values and warn here instead of error
+        if ! isempty(spans)
+            error("Column $colname has missing values: $spans")
+        end
+    end
 end
 
 """
@@ -33,97 +96,103 @@ mapped to dates.
 """
 function initweather(weatherfile::String, startdate::Date, enddate::Date)
     @debug "Initialising weather"
-    data = CSV.File(weatherfile, missingstring="NA", dateformat="yyyy-mm-dd",
-                    types=[Date, Float64, Float64, Float64, Float64, Float64,
-                           Float64, Float64, Float64, Float64])
-    weather = Dict{Date,Weather}()
-    for row in data
-        if row.date >= startdate && row.date <= enddate
-            weather[row.date] = Weather(row.mean_windspeed, row.precipitation,
-                                        row.sunshine_hours, row.mean_cloud_cover,
-                                        row.mean_humidity, row.mean_temperature,
-                                        row.max_temperature, row.min_temperature,
-                                        row.potential_evapotranspiration)
-        end
+    if startdate > enddate
+        error("Startdate is after enddate: $startdate > $enddate.")
     end
-    if length(weather) <= Dates.value(enddate-startdate)
-        @warn ("There are missing days in the weather input file:"
-               * " expected $(Dates.value(enddate-startdate) + 1), got $(length(weather)).")
+    Tfloat = Union{Missing, Float64}
+    df = CSV.read(weatherfile, DataFrame;
+                  missingstring="NA", dateformat="yyyy-mm-dd",
+                  types=[Date, Tfloat, Tfloat, Tfloat, Tfloat, Tfloat,
+                         Tfloat, Tfloat, Tfloat, Tfloat])
+    needed_cols = [
+        "date", "mean_windspeed", "precipitation", "sunshine_hours",
+        "mean_cloud_cover", "mean_humidity", "mean_temperature",
+        "max_temperature", "min_temperature", "potential_evapotranspiration"
+    ]
+    if ! issubset(needed_cols, names(df))
+        error("Missing columns in weather data file: $(setdiff(needed_cols, names(df)))")
     end
-    weather
+    filter!(r -> r.date >= startdate && r.date <= enddate, df)
+    sort!(df, :date)
+    if startdate:enddate != df.date
+        error("There are missing days (rows) in the weather data file."
+              * " expected $(Dates.value(enddate - startdate) + 1), got $(nrow(df)))."
+              * " Missing dates are $(setdiff(startdate:enddate, df.date)).")
+    end
+    check_missing_weatherdata(df)
+    return Weather(; firstdate = minimum(df.date),
+                   lastdate = maximum(df.date),
+                   windspeed = df.mean_windspeed,
+                   precipitation = df.precipitation,
+                   sunshine = df.sunshine_hours,
+                   cloudcover = df.mean_cloud_cover,
+                   humidity = df.mean_humidity,
+                   meantemp = df.mean_temperature,
+                   maxtemp = df.max_temperature,
+                   mintemp = df.min_temperature,
+                   evapotranspiration = df.potential_evapotranspiration)
 end
 
 """
-    windspeed(model)
+    windspeed(weather, date)
 
-Return today's average windspeed in m/s.
+Return the average windspeed in m/s on `date`.
 """
-function windspeed(model::SimulationModel)
-    model.weather[model.date].windspeed
-end
+windspeed(weather::Weather, date::Date) =
+    weather.windspeed[daynumber(weather, date)]
 
 """
-    precipitation(model)
+    precipitation(weather, date)
 
-Return today's total precipitation in mm.
+Return the total precipitation in mm on `date`.
 """
-function precipitation(model::SimulationModel)
-    model.weather[model.date].precipitation
-end
+precipitation(weather::Weather, date::Date) =
+    weather.precipitation[daynumber(weather, date)]
 
 """
-    sunshine(model)
+    sunshine(weather, date)
 
-Return today's sunshine duration in hours.
+Return the sunshine duration in hours on `date`.
 """
-function sunshine(model::SimulationModel)
-    model.weather[model.date].sunshine
-end
+sunshine(weather::Weather, date::Date) =
+    weather.sunshine[daynumber(weather, date)]
 
 """
-    humidity(model)
+    humidity(weather, date)
 
 Return today's average vapour pressure in %.
 """
-function humidity(model::SimulationModel)
-    model.weather[model.date].humidity
-end
+humidity(weather::Weather, date::Date) =
+    weather.humidity[daynumber(weather, date)]
 
 """
-    meantemp(model)
+    meantemp(weather, date)
 
-Return today's mean temperature in °C.
+Return the mean temperature in °C on `date`.
 """
-function meantemp(model::SimulationModel)
-    model.weather[model.date].meantemp
-end
+meantemp(weather::Weather, date::Date) =
+    weather.meantemp[daynumber(weather, date)]
 
 """
-    maxtemp(model)
+    maxtemp(weather, date)
 
-Return today's maximum temperature in °C.
+Return the maximum temperature in °C on `date`.
 """
-function maxtemp(model::SimulationModel)
-    model.weather[model.date].maxtemp
-end
+maxtemp(weather::Weather, date::Date) =
+    weather.maxtemp[daynumber(weather, date)]
 
 """
-    mintemp(model)
+    mintemp(weather, date)
 
-Return today's minimum temperature in °C.
+Return the minimum temperature in °C on `date`.
 """
-function mintemp(model::SimulationModel)
-    model.weather[model.date].mintemp
-end
+mintemp(weather::Weather, date::Date) =
+    weather.mintemp[daynumber(weather, date)]
 
 """
-    evapotranspiration(model)
+    evapotranspiration(weather, date)
 
-Return today's potential evapotranspiration (ETo), based on
-the
+Return today's potential evapotranspiration (ETo) on `date`.
 """
-function evapotranspiration(model::SimulationModel)
-    model.weather[model.date].evapotranspiration
-end
-
-#TODO add functions for evapotranspiration
+evapotranspiration(weather::Weather, date::Date) =
+    weather.evapotranspiration[daynumber(weather, date)]

test/io_tests.jl

@@ -9,8 +9,8 @@
     @test @param(core.startdate) == Date(2022, 2, 1)
     @test @param(core.loglevel) == "warn"
     @test @param(nature.targetspecies) == ["Wolpertinger", "Wyvern"]
-    @param(core.enddate) = Date(2022,1,3)
-    @test @param(core.enddate) == Date(2022,1,3)
+    @param(core.enddate) = Date(2022,2,3)
+    @test @param(core.enddate) == Date(2022,2,3)
 end
 
 @testset "Output functions" begin

test/landscape_tests.jl

@@ -64,10 +64,10 @@ end
     @test Ps.distancetoedge((6,6), model) == 20m
 end
 
-@testset "Weather initialisation" begin
+@testset "Weather interface" begin
     # these tests are specific to the Jena weather file
     model = inittestmodel()
-    @test length(keys(model.weather)) == 59
+    @test length(model.weather) == 59
     @test ismissing(Ps.windspeed(model))
     @test Ps.precipitation(model) == 1.3
     @test ismissing(Ps.sunshine(model))

test/runtests.jl

@@ -45,22 +45,14 @@ function inittestmodel(smallmap=true)
                              TESTSETTINGS["core.enddate"])
     # TODO: support other crop models besides ALMaSS
     crops = Ps.ALMaSS.readcropparameters(TESTSETTINGS["crop.cropdirectory"])
-    farmers = Vector{Farmer}()
-    farmplots = Vector{Ps.FarmPlot{Ps.ALMaSS.CropState}}()
-    model = AgricultureModel{Ps.ALMaSS.CropType,Ps.ALMaSS.CropState}(
-        TESTSETTINGS,
-        StableRNG(TESTSETTINGS["core.seed"]),
-        global_logger(),
-        Dict{String,DataOutput}(),
-        TESTSETTINGS["core.startdate"],
+    model = AgricultureModel{Ps.ALMaSS.CropType,Ps.ALMaSS.CropState}(;
+        settings = copy(TESTSETTINGS),
+        rng = StableRNG(TESTSETTINGS["core.seed"]),
+        logger = global_logger(),
+        date = TESTSETTINGS["core.startdate"],
         landscape,
         weather,
-        crops,
-        farmers,
-        farmplots,
-        Vector{Union{Animal,Nothing}}(),
-        Vector{Pair{Animal, Date}}(),
-        Vector{FarmEvent}()
+        crops
     )
     model
 end    
@@ -100,6 +92,9 @@ end
         include("landscape_tests.jl")
         include("simulation_tests.jl")
     end
+    @testset "Weather model" begin
+        include("weather_tests.jl")
+    end
     @testset "Nature model" begin
         include("nature_tests.jl")
     end

test/weather_tests.jl

+
+@testset "Constructor and interface" begin
+    float3 = [0.0, 1.0, 2.0]
+    missing3 = Union{Missing,Float64}[missing, missing, missing]
+    weather = Weather(; firstdate = Date("2000-01-01"),
+                      lastdate = Date("2000-01-03"),
+                      windspeed = missing3,
+                      precipitation = float3,
+                      sunshine = missing3,
+                      cloudcover = missing3,
+                      humidity = missing3,
+                      meantemp = float3,
+                      maxtemp = float3,
+                      mintemp = float3,
+                      evapotranspiration = float3)
+    date = Date("2000-01-02")
+    @test length(weather) == 3
+    @test Ps.daynumber(weather, date) == 2
+    @test Ps.windspeed(weather, date) |> ismissing
+    @test Ps.precipitation(weather, date) == 1.0
+    @test Ps.sunshine(weather, date) |> ismissing
+    @test Ps.humidity(weather, date) |> ismissing
+    @test Ps.meantemp(weather, date) == 1.0
+    @test Ps.maxtemp(weather, date) == 1.0
+    @test Ps.mintemp(weather, date) == 1.0
+    @test Ps.evapotranspiration(weather, date) == 1.0
+end