Commit 72a82dbb authored by Sebastian Ostberg's avatar Sebastian Ostberg
Browse files

Initial commit including scripts

parent 90848c4a
# WM_LPJmL_indicators
# World Modelers LPJmL crop yield indicators
This project contains code to aggregate raw outputs from LPJmL crop yield forecasts and/or historical yield simulations into indicators that are more readily usable.
\ No newline at end of file
This project contains code to aggregate raw outputs from LPJmL crop yield forecasts and/or historical yield simulations into indicators that are more readily usable.
## Functionality
- `average_yield.R`: function `average_yield()` returns an array of crop data averaged across time. Spatial and crop dimensions are kept
- `download_yield_data.R`: function `download_yield_data()` attempts to download the raw output file for one LPJmL simulation (either historical time series or forecast) from the remote server and returns the local file if downloaded successfully.
- `get_forecast_options.R`: function `get_forecast_options` returns a list of character vectors with elements 'irrigation_options', 'fertilizer_options', 'sowing_options' giving all available options on the remote server
- `get_weather_options.R`: function `get_weather_options` returns either a character vector of weather options available on the remote server for a specific combination of irrigation, fertilizer and sowing option, or a list containing a character vector of weather options and a character vector with corresponding remote file names.
- `rsync_setup.R`: sets up basic information to retrieve data from remote server
## Indicators
### risk_yield_loss.R
This indicator describes the risk that forecasted yields will fall below historical yield levels by more than a certain percentage threshold. The risk is calculated based on the fraction of weather options that show the respective yield loss.
User-specified parameters:
- `hist_ref_startyear` and `hist_ref_endyear`: start and end of the historical reference period against which forecasts are compared (default: 2016-2020)
- `yield_loss`: yield loss threshold in percent compared to historical reference to be counted towards yield loss risk (default: 10)
- `loss_probability`: risk threshold in %, i.e. the minimum fraction of all weather options that need to show a yield loss greater than the `yield_loss` threshold (default 66.7)
- `forecast_year`: the forecast year to be compared against the historical reference yield (default: 2021)
- `irrigation_scenario`: forecast intervention option (only "reference" available for August 2021 experiment)
- `fertilizer_scenario`: forecast intervention option (default: reference)
- `sowing_scenario`: forecast intervention option (default: reference)
- `experiment`: allows future extension to LPJmL results from different experiments (only "August2021_experiment" currently available)
The script first downloads LPJmL outputs for the historical time period from the remote server and calculates for each crop and each grid cell the average yield during the reference period specified by `hist_ref_startyear` and `hist_ref_endyear`.
In a second step, outputs from LPJmL crop yield forecast simulations are downloaded for all weather options available on the remote server for the specified combination of `irrigation_scenario`, `fertilizer_scenario` and `sowing_scenario`. Forecasted yield during the `forecast_year` is calculated for each crop and and grid cell and compared against the historical reference yield. This risk of yield loss is calculated for each crop and grid cell as the fraction of weather options that show a yield loss (compared to reference yield) of at least `yield_loss` percent.
To do: Grid-cell level risks are aggregated to the country level by summing grid cells where the risk of yield loss is at least `loss_probability` percent. The area at risk at the country level will be specified as absolute area (in hectare) as well as in percent of total crop area (e.g. 25% of all rainfed wheat areas in the country).
To do: Confirm format of returned risk indicator with CauseMos.
#' @title Calculate average yield
#'
#' @description Calculates the average across time for a crop data time series from LPJmL.
#' The function always calculates the mean across different years. Multiple time steps within
#' the same year may be either aggregated as mean or as sum.
#'
#' @param filename File name of the NetCDF file containing the time series data
#' @param avg_variable Name of the NetCDF variable containing the crop time series data. This should
#' normally be crop yield but could also be other crop variables such as crop production or
#' harvested area
#' @param avg_startyear Beginning of period used for averaging
#' @param avg_endyear End of period used for averaging
#' @param multiple_per_year Function to use to aggregate multiple time steps within the same year
#' (either "mean" or "sum", default: "mean")
#' @param multiple_weight_variable Optional name of NetCDF variable used as weights if
#' multiple_per_year == "mean" (default: NULL)
#'
#' @return Array of crop data averaged across time. Spatial and crop dimensions are kept
#'
#' @export
average_yield <- function(filename, avg_variable, avg_startyear, avg_endyear, multiple_per_year = "mean", multiple_weight_variable = NULL) {
## check for valid setting for multiple_per_year
if(!multiple_per_year %in% c("mean", "sum")) {
stop(paste("multiple_per_year must be either 'mean' or 'sum'. You provided:", sQuote(multiple_per_year)))
}
## check that file exists
if(!file.exists(filename)) {
stop(paste(sQuote(filename), "does not exist."))
}
## open file
nc <- nc_open(filename)
## check if avg_variable is available in file
if(!avg_variable %in% names(nc$var)) {
stop(paste("Variable", sQuote(avg_variable), "not available in file", filename))
}
if(!is.null(multiple_weight_variable) && !multiple_weight_variable %in% names(nc$var)) {
stop(paste("Variable", sQuote(multiple_weight_variable), "not available in file", filename))
}
## get time
# unit string in file
timeunit <- nc$dim$time$units
# Is time unit relative to an origin?
if(grepl("since", timeunit)) {
# detect origin date
origin <- as.POSIXct(sub("years|days|seconds since ", "", timeunit), tz="UTC")
# detect unit of offset compared to origin
origin_unit <- regmatches(timeunit, regexec("^([a-z]+) since", timeunit))[[1]][-1]
time_is_relative <- TRUE
} else {
time_is_relative <- FALSE
if(!grepl("year", timeunit)) {
# currently only support "year" as absolute time unit
stop(paste("Time unit", sQuote(timeunit), "not supported in file", filename))
}
}
# load values of time axis
timevals <- ncvar_get(nc, "time")
## find timesteps within period avg_startyear to avg_endyear
if(!time_is_relative) {
# timevals equal year values
avg_steps <- which(timevals >= avg_startyear & timevals <= avg_endyear)
if(length(avg_steps) < 1) {
stop(paste0("avg_startyear ", sQuote(avg_startyear), " or avg_endyear ", sQuote(avg_endyear), " outside of time period covered covered by file ", filename, " (", paste(range(timevals), collapse=("-")), ")"))
}
} else {
# need to convert timevals into year values using origin and origin_unit
year_timeval <- switch(origin_unit,
years = as.double(format(origin, "%Y"))+timevals,
days = as.double(format(as.Date(origin)+timevals, "%Y")),
seconds = as.double(format(origin+timevals, "%Y")))
avg_steps <- which(year_timeval >= avg_startyear & year_timeval <= avg_endyear)
if(length(avg_steps) < 1) {
stop(paste0("Start year ", sQuote(avg_startyear), " or end year ", sQuote(avg_endyear), " outside of time period covered covered by file ", filename, " (", paste(range(year_timeval), collapse=("-")), ")"))
}
# set timevals to year values
timevals <- year_timeval
}
if(avg_startyear < min(timevals)) {
stop(paste("avg_startyear", avg_startyear, "is earlier than period covered by of file", filename))
}
if(avg_endyear > max(timevals)) {
stop(paste("avg_endyear", avg_endyear, "is later than period covered by file", filename))
}
## names of variable dimensions
var_dims <- names(nc$dim)[nc$var[[avg_variable]]$dimids+1]
## index of time dimension in variable dimensions
time_dim <- which(var_dims == "time")
var_dims_has_names <- rep(FALSE, length(var_dims))
names(var_dims_has_names) <- names(var_dims)
## find dimensions for which name attributes are saved as character variables
for(namevar in names(which(apply(sapply(var_dims, grepl, x = grep("Name", names(nc$var), value = TRUE, ignore.case = TRUE), ignore.case = TRUE), 2, any)))) {
nc_var <- intersect(grep("Name", names(nc$var), value = TRUE, ignore.case = TRUE), grep(namevar, names(nc$var), value = TRUE, ignore.case = TRUE))
name_values <- ncvar_get(nc, nc_var)
assign(paste0(namevar, "_names"), name_values)
var_dims_has_names[namevar] <- TRUE
}
result_array <- array(dim=nc$var[[avg_variable]]$varsize[-time_dim])
for(namevar in names(which(var_dims_has_names))) {
dimnames(result_array)[[which(var_dims[-time_dim] == namevar)]] <- get(paste0(namevar, "_names"))
}
nvalid <- array(0, dim=nc$var[[avg_variable]]$varsize[-time_dim])
for(y in avg_startyear:avg_endyear) {
year_result <- array(dim=nc$var[[avg_variable]]$varsize[-time_dim])
year_nvalid <- array(0, dim=nc$var[[avg_variable]]$varsize[-time_dim])
for(yearstep in which(timevals==y)) {
filedata <- ncvar_get(nc, avg_variable, start=ifelse(var_dims == "time", yearstep, 1), count=ifelse(var_dims == "time", 1, -1))
valid_index <- which(!is.na(filedata))
update_index <- intersect(which(!is.na(year_result)), valid_index)
set_index <- intersect(which(is.na(year_result)), valid_index)
if(multiple_per_year == "mean" && !is.null(multiple_weight_variable)) {
## Normally harvested area for several timesteps within the same year
## should be identical for LPJmL results.However, if specified
## multiple_weight_variable is used as weight. This variable must have
## the same dimensions as avg_variable
weightdata <- ncvar_get(nc, multiple_weight_variable, start=ifelse(var_dims == "time", yearstep, 1), count=ifelse(var_dims == "time", 1, -1))
if(anyNA(weightdata[valid_index])) {
stop(paste("Inconsistency between", sQuote(avg_variable), "and", sQuote(multiple_weight_variable), "data in year", y))
}
weightdata[-valid_index] <- NA
year_nvalid[valid_index] <- year_nvalid[valid_index] + weightdata[valid_index]
year_result[set_index] <- filedata[set_index] * weightdata[set_index]
year_result[update_index] <- year_result[update_index] + filedata[update_index] * weightdata[update_index]
} else {
year_nvalid[valid_index] <- year_nvalid[valid_index]+1
year_result[set_index] <- filedata[set_index]
year_result[update_index] <- year_result[update_index] + filedata[update_index]
}
}
if(multiple_per_year == "mean") {
year_result <- year_result/year_nvalid
year_result[which(year_nvalid == 0)] <- NA
}
valid_index <- which(!is.na(year_result))
update_index <- intersect(which(!is.na(result_array)), valid_index)
set_index <- intersect(which(is.na(result_array)), valid_index)
nvalid[valid_index] <- nvalid[valid_index]+1
result_array[set_index] <- year_result[set_index]
result_array[update_index] <- result_array[update_index] + year_result[update_index]
}
result_array <- result_array/nvalid
result_array[which(nvalid == 0)] <- NA
nc_close(nc)
rm(nvalid, year_result, year_nvalid)
return(result_array)
}
#' @title Download output crop output data for LPJmL from PIK rsync server
#'
#' @description This function can be used to download historical time
#' series data as well as forecast data for LPJmL
#'
#' @param runtype Type of LPJmL crop simulation. Either 'historical' or 'forecast'
#' @param localname Local file name to be used for downloaded data
#' @param experiment The name of the set of crop simulations
#' @param irrigation_scenario Irrigation scenario for which to download
#' data (if runtype == 'forecast')
#' @param fertilizer_scenario Fertilizer scenario for which to download
#' data (if runtype == 'forecast')
#' @param sowing_scenario Sowing date scenario for which to download
#' data (if runtype == 'forecast')
#'
#' @return The name of the local file if downloaded successfully. This
#' function will throw an error if any of the provided parameters are invalid.
#'
#' @export
download_yield_data <- function(runtype, localname, experiment = "August2021_experiment", irrigation_scenario = NULL, fertilizer_scenario = NULL, sowing_scenario = NULL, weather_year = NULL) {
if(!runtype %in% c("historical", "forecast")) {
stop(paste("Invalid runtype", sQuote(runtype), "provided. Must be 'historical' or 'forecast'"))
}
## check for valid experiment
if((runtype == "historical" && !experiment %in% names(historical_remote_dir)) || (runtype == "forecast" && !experiment %in% names(forecast_remote_dir))) {
stop(paste("Invalid experiment", sQuote(experiment)))
}
if(runtype == "historical") {
rsync_url <- paste0(rsync_baseurl, experiment, "/lpjml-results/", historical_remote_dir[experiment], "/")
## get contents of remote directory
remote_contents <- system2(rsync_command, args=c(rsync_args, rsync_url), stdout=TRUE)
source_file <- grep(".nc$", unname(sapply(sapply(remote_contents, strsplit, split = " "), function(indata) return(indata[length(indata)]))), value = TRUE)
if(length(source_file) != 1) {
stop(paste("Cannot determine source file from contents of source directory:", toString(sQuote(unname(sapply(sapply(remote_contents, strsplit, split = " "), function(indata) return(indata[length(indata)])))))))
}
} else if(runtype == "forecast") {
for(opt in c("irrigation_scenario", "fertilizer_scenario", "sowing_scenario", "weather_year")) {
if(!exists(opt) || is.null(get(opt))) {
stop(paste(sQuote(opt), "not set"))
}
}
weather_options <- get_weather_options(experiment = experiment,
remote_dir = forecast_remote_dir[experiment],
irrigation_scenario = irrigation_scenario,
fertilizer_scenario = fertilizer_scenario,
sowing_scenario = sowing_scenario,
include_filenames = TRUE)
if(!weather_year %in% weather_options$weather_options) {
stop(paste("Invalid weather_year", sQuote(weather_year)))
}
source_file <- grep(weather_year, weather_options$filenames, value = TRUE)
if(length(source_file) != 1) {
stop(paste0("Cannot determine matching filename for irrigation_scenario ",
sQuote(irrigation_scenario),
", fertilizer_scenario ",
sQuote(fertilizer_scenario),
", sowing_scenario ",
sQuote(sowing_scenario),
", weather_year ", sQuote(weather_year)))
}
rsync_url <- paste0(rsync_baseurl, experiment, "/lpjml-results/", forecast_remote_dir[experiment], "/irrigation_target_", irrigation_scenario, "_n_rate_shift_", fertilizer_scenario, "_sowing_delta_", sowing_scenario, "/")
} else {
stop(paste("runtype", sQuote(runtype), "not defined"))
}
rsync_url <- paste0(rsync_url, source_file)
## Download file
system2(rsync_command, args=c(rsync_args, rsync_url, localname))
if(file.exists(localname)) {
return(localname)
} else {
return(FALSE)
}
}
#' @title Get available intervention scenario options for a set of
#' forecast simulations
#'
#' @description Finds all available irrigation, fertilizer and sowing
#' date scenario options available for a set of forecast simulations
#'
#' @param experiment The name of the set of forecast simulations
#' @param remote_dir Remote directory on the rsync server containing
#' the forecast files
#'
#' @return A list of character vectors with elements 'irrigation_options',
#' 'fertilizer_options', 'sowing_options'
#'
#' @export
get_forecast_options <- function(experiment, remote_dir) {
if(!experiment %in% names(remote_dir)) {
stop(paste("Invalid experiment", sQuote(experiment)))
}
## parts of rsync command
rsync_url <- paste0(rsync_baseurl, experiment, "/lpjml-results/", remote_dir[experiment], "/")
## get contents of remote directory
remote_contents <- system2(rsync_command, args=c(rsync_args, rsync_url), stdout=TRUE)
forecast_versions <- unname(sapply(sapply(remote_contents, strsplit, split = " "), function(indata) return(indata[length(indata)])))
forecast_versions <- forecast_versions[grep("irrigation", forecast_versions)]
scenario_options <- regmatches(forecast_versions, regexec("^irrigation_target_(.+)_n_rate_shift_(.+)_sowing_delta_(.+)", forecast_versions))
if(any(sapply(scenario_options, length) != 4)) {
stop(paste("Invalid forecast scenario(s):", toString(sQuote(sapply(scenario_options[which(sapply(scenario_options, length) != 4)], function(indata) indata[1])))))
}
irrigation_options <- unique(sapply(scenario_options, function(indata) indata[2]))
fertilizer_options <- unique(sapply(scenario_options, function(indata) indata[3]))
sowing_options <- unique(sapply(scenario_options, function(indata) indata[4]))
return(list(irrigation_options = irrigation_options, fertilizer_options = fertilizer_options, sowing_options = sowing_options))
}
#' @title Get available weather options for a specific forecast set
#'
#' @description Finds all available weather options corresponding to a
#' specific set of LPJmL forecasts.
#'
#' @param experiment The name of the set of forecast simulations
#' @param remote_dir Remote directory on the rsync server containing
#' the forecast files
#' @param irrigation_scenario Irrigation scenario for which to query
#' weather options
#' @param fertilizer_scenario Fertilizer scenario for which to query
#' weather options
#' @param sowing_scenario Sowing date scenario for which to query
#' weather options
#' @param include_filenames If set to TRUE, include file names
#' corresponding to weather options
#'
#' @return A character vector of weather options if include_filenames is FALSE.
#' Otherwise a list containing a character vector of weather options and
#' a character vector with corresponding file names. Function fails if
#' any of irrigation_scenario, fertilizer_scenario or sowing_scenario
#' does not exist
#'
#' @export
get_weather_options <- function(experiment, remote_dir, irrigation_scenario, fertilizer_scenario, sowing_scenario, include_filenames = FALSE) {
scenario_options <- get_forecast_options(experiment, remote_dir)
irrigation_options <- scenario_options$irrigation_options
fertilizer_options <- scenario_options$fertilizer_options
sowing_options <- scenario_options$sowing_options
if(!irrigation_scenario %in% irrigation_options) {
stop(paste("Invalid irrigation_scenario", sQuote(irrigation_scenario)))
}
if(!fertilizer_scenario %in% fertilizer_options) {
stop(paste("Invalid fertilizer_scenario", sQuote(fertilizer_scenario)))
}
if(!sowing_scenario %in% sowing_options) {
stop(paste("Invalid sowing_scenario", sQuote(sowing_scenario)))
}
## get contents of forecast version directory
rsync_url <- paste0(rsync_baseurl, experiment, "/lpjml-results/", remote_dir, "/irrigation_target_", irrigation_scenario, "_n_rate_shift_", fertilizer_scenario, "_sowing_delta_", sowing_scenario, "/")
remote_contents <- system2(rsync_command, args=c(rsync_args, rsync_url), stdout=TRUE)
forecast_versions <- unname(sapply(sapply(remote_contents, strsplit, split = " "), function(indata) return(indata[length(indata)])))
forecast_versions <- forecast_versions[grep("forecast", forecast_versions)]
weather_options <- regmatches(forecast_versions, regexec("_wy([0-9]{4}-[0-9]{4}).nc", forecast_versions))
if(any(sapply(weather_options, length) != 2)) {
stop(paste("Cannot determine weather years from forecast version(s):", toString(sQuote(forecast_versions[which(sapply(weather_options, length) != 2)]))))
}
weather_options <- sapply(weather_options, function(indata) indata[2])
if(include_filenames) {
return(list(weather_options = weather_options, filenames = forecast_versions))
} else {
return(weather_options)
}
}
rm(list=ls())
## default values for command line parameters
# experiment name (must be defined in rsync_setup.R)
experiment <- "August2021_experiment"
# reference period
hist_ref_startyear <- 2016
hist_ref_endyear <- 2020
# threshold for yield loss compared to average yield during reference period
yield_loss <- 10 # in percent
# fraction (in %) of weather options that must show yield loss of at least yield_loss %
loss_probability <- (2/3)*100
# forecast year (2021 or 2022) which is compared to historical yield
forecast_year <- 2021
# intervention scenario options
irrigation_scenario <- "reference"
fertilizer_scenario <- "reference"
sowing_scenario <- "reference"
## load R packages
library(ncdf4)
## read command line parameters (if any were provided)
if(length(commandArgs(trailingOnly=TRUE)) > 0) {
if(any(grepl("hist_ref_startyear", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("hist_ref_startyear", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
startyear <- as.integer(grep("hist_ref_startyear", unlist(param), invert=T, value=TRUE))
if(is.finite(startyear)) {
print(paste("Parameter hist_ref_startyear provided as command line argument:", startyear))
hist_ref_startyear <- startyear
} else {
stop(paste("Invalid parameter hist_ref_startyear", sQuote(startyear)))
}
}
if(any(grepl("hist_ref_endyear", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("hist_ref_endyear", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
endyear <- as.integer(grep("hist_ref_endyear", unlist(param), invert=T, value=TRUE))
if(is.finite(startyear)) {
print(paste("Parameter hist_ref_endyear provided as command line argument:", endyear))
hist_ref_endyear <- endyear
} else {
stop(paste("Invalid parameter hist_ref_endyear", sQuote(endyear)))
}
}
if(any(grepl("yield_loss", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("yield_loss", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
loss <- as.double(grep("yield_loss", unlist(param), invert=T, value=TRUE))
if(is.finite(loss)) {
print(paste("Parameter yield_loss provided as command line argument:", loss))
yield_loss <- loss
} else {
stop(paste("Invalid parameter yield_loss", sQuote(loss)))
}
}
if(any(grepl("loss_probability", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("loss_probability", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
loss <- as.double(grep("loss_probability", unlist(param), invert=T, value=TRUE))
print(paste("Parameter loss_probability provided as command line argument:", loss))
loss_probability <- loss
}
if(any(grepl("forecast_year", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("forecast_year", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
year <- as.integer(grep("forecast_year", unlist(param), invert=T, value=TRUE))
if(is.finite(year)) {
print(paste("Parameter forecast_year provided as command line argument:", year))
forecast_year <- year
} else {
stop(paste("Invalid parameter forecast_year", sQuote(year)))
}
}
if(any(grepl("irrigation_scenario", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("irrigation_scenario", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
scen <- grep("irrigation_scenario", unlist(param), invert=T, value=TRUE)
if(nchar(scen)>0) {
print(paste("Parameter irrigation_scenario provided as command line argument:", scen))
irrigation_scenario <- scen
} else {
stop(paste("Invalid parameter irrigation_scenario", sQuote(scen)))
}
}
if(any(grepl("fertilizer_scenario", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("fertilizer_scenario", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
scen <- grep("fertilizer_scenario", unlist(param), invert=T, value=TRUE)
if(nchar(scen)>0) {
print(paste("Parameter fertilizer_scenario provided as command line argument:", scen))
fertilizer_scenario <- scen
} else {
stop(paste("Invalid parameter fertilizer_scenario", sQuote(scen)))
}
}
if(any(grepl("sowing_scenario", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("sowing_scenario", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
scen <- grep("sowing_scenario", unlist(param), invert=T, value=TRUE)
if(nchar(scen)>0) {
print(paste("Parameter sowing_scenario provided as command line argument:", scen))
sowing_scenario <- scen
} else {
stop(paste("Invalid parameter sowing_scenario", sQuote(scen)))
}
}
if(any(grepl("experiment", commandArgs(trailingOnly=TRUE)))) {
param <- strsplit(grep("experiment", commandArgs(trailingOnly=TRUE), value=TRUE), split="=")
exp <- grep("experiment", unlist(param), invert=T, value=TRUE)
if(nchar(exp)>0) {
print(paste("Parameter experiment provided as command line argument:", exp))
experiment <- exp
} else {
stop(paste("Invalid parameter experiment", sQuote(exp)))
}
}
}
## rsync functionality
source("rsync_setup.R")
source("get_forecast_options.R")
source("get_weather_options.R")
source("download_yield_data.R")
source("average_yield.R")
if(!experiment %in% names(historical_remote_dir) || !experiment %in% names(forecast_remote_dir)) {
stop(paste("Invalid experiment", sQuote(experiment)))
}
## download historical time series data
if(!file.exists("output")) {
dir.create("output")
}
hist_timeseries_filename <- download_yield_data(runtype = "historical",
localname = "output/hist_timeseries.nc",
experiment = experiment,
irrigation_scenario = NULL,
fertilizer_scenario = NULL,
sowing_scenario = NULL,
weather_year = NULL)
if(!is.character(hist_timeseries_filename) || !file.exists(hist_timeseries_filename)) {
stop("Error downloading historical time series data to calculate reference yield.")
} else {
print(paste("Historical yield time series downloaded to temporary file", hist_timeseries_filename))
}
## calculate average yield for reference period
print(paste0("Calculate average yield for reference period ", hist_ref_startyear, "-", hist_ref_endyear))
hist_ref_yield <- average_yield(filename = hist_timeseries_filename,
avg_variable = "crop_yield",
avg_startyear = hist_ref_startyear,
avg_endyear = hist_ref_endyear,
multiple_per_year = "mean",
multiple_weight_variable = "harvested_area")
## calculate threshold for yield loss
if(yield_loss < 0 || yield_loss > 100) {
stop(paste("yield_loss threshold must be between 0 and 100%. Provided value:", yield_loss))
}
yield_loss_threshold <- hist_ref_yield*(100-yield_loss)/100
## delete dowloaded file to conserve space
file.remove(hist_timeseries_filename)
## get all available weather options for selected forecast scenario options
forecast_weather_options <- get_weather_options(experiment = experiment,
remote_dir = forecast_remote_dir[experiment],
irrigation_scenario = irrigation_scenario,
fertilizer_scenario = fertilizer_scenario,
sowing_scenario = sowing_scenario,
include_filenames = FALSE)
## check how many weather options fall below yield_loss_threshold
forecast_below_threshold_count <- array(dim=dim(yield_loss_threshold), dimnames=dimnames(yield_loss_threshold))
harvested_area_forecast_below_threshold <- array(dim=dim(yield_loss_threshold), dimnames=dimnames(yield_loss_threshold))
harvested_area_all <- array(dim=dim(yield_loss_threshold), dimnames=dimnames(yield_loss_threshold))
harvested_area_all_count <- array(0, dim=dim(yield_loss_threshold), dimnames=dimnames(yield_loss_threshold))
for(wy in forecast_weather_options) {
forecast_filename <- download_yield_data(runtype = "forecast",
localname = "output/forecast.nc",
experiment = experiment,
irrigation_scenario = irrigation_scenario,
fertilizer_scenario = fertilizer_scenario,
sowing_scenario = sowing_scenario,
weather_year = wy)
if(!is.character(forecast_filename) || !file.exists(forecast_filename)) {
stop(paste("Error downloading forecast data for weather option", sQuote(wy)))
} else {
print(paste("Forecast data for weather option", sQuote(wy), "downloaded to temporary file", forecast_filename))
}
print(paste("Calculate yield change in", forecast_year, "compared to historical reference yield"))
forecast_yield <- average_yield(filename = forecast_filename,
avg_variable = "crop_yield",
avg_startyear = forecast_year,
avg_endyear = forecast_year,
multiple_per_year = "sum",
multiple_weight_variable = NULL)
harvested_area <- average_yield(filename = forecast_filename,
avg_variable = "harvested_area",
avg_startyear = forecast_year,
avg_endyear = forecast_year,
multiple_per_year = "sum",
multiple_weight_variable = NULL)
## delete downloaded file to ensure next weather option is downloaded correctly
file.remove(forecast_filename)
if(any(dim(forecast_yield) != dim(yield_loss_threshold))) {
stop("Forecast and reference yield data have different dimensions")
}
if(any(dim(harvested_area) != dim(yield_loss_threshold))) {
stop("Forecast harvested area and reference yield data have different dimensions")
}
## harmonize sequence of crops and cropping systems (rainfed/irrigated)
if(!is.null(dimnames(forecast_yield)[[3]]) && !is.null(dimnames(yield_loss_threshold)[[3]])) {
if(length(intersect(dimnames(forecast_yield)[[3]], dimnames(yield_loss_threshold)[[3]])) != dim(yield_loss_threshold)[3]) {
stop("Crop mismatch between forecast and reference yield data")
}
bakdim <- dim(forecast_yield)
bakdimnames <- dimnames(forecast_yield)
bakdimnames[[3]] <- dimnames(yield_loss_threshold)[[3]]
forecast_yield <- forecast_yield[,,match(dimnames(yield_loss_threshold)[[3]], dimnames(forecast_yield)[[3]]),]
dim(forecast_yield) <- bakdim
dimnames(forecast_yield) <- bakdimnames
}
if(!is.null(dimnames(forecast_yield)[[4]]) && !is.null(dimnames(yield_loss_threshold)[[4]])) {
if(length(intersect(dimnames(forecast_yield)[[4]], dimnames(yield_loss_threshold)[[4]])) != dim(yield_loss_threshold)[4]) {
stop("Cropping system mismatch between forecast and reference yield data")
}
bakdim <- dim(forecast_yield)
bakdimnames <- dimnames(forecast_yield)
bakdimnames[[4]] <- dimnames(yield_loss_threshold)[[4]]
forecast_yield <- forecast_yield[,,,match(dimnames(yield_loss_threshold)[[4]], dimnames(forecast_yield)[[4]])]
dim(forecast_yield) <- bakdim
dimnames(forecast_yield) <- bakdimnames
}
if(!is.null(dimnames(harvested_area)[[3]]) && !is.null(dimnames(yield_loss_threshold)[[3]])) {
if(length(intersect(dimnames(harvested_area)[[3]], dimnames(yield_loss_threshold)[[3]])) !=