Plots and analyses produced for ‘The Weather, land and crops in the Indus Village model’ (Quaternary, 2021)
Andreas Angourakis
07/01/2022
This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/Quaternary_Angourakis-et-al-2022
Load packages and define helper functions:
Load packages:
require(pander)## Loading required package: pander## Warning: package 'pander' was built under R version 3.6.3require(reshape2)## Loading required package: reshape2## Warning: package 'reshape2' was built under R version 3.6.3require(ggplot2)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.6.3require(patchwork)## Loading required package: patchwork## Warning: package 'patchwork' was built under R version 3.6.3require(colorspace)## Loading required package: colorspace## Warning: package 'colorspace' was built under R version 3.6.3require(grid)## Loading required package: gridrequire(gridExtra)## Loading required package: gridExtra## Warning: package 'gridExtra' was built under R version 3.6.3require(png)## Loading required package: pngDefine lookup functions to access data between tables using identifiers or ‘keys’:
getRelatedVariableWithKey <- function(keyInCallingTable, keyInCalledTable, variableInCalledTable)
{
return(
sapply(keyInCallingTable, function(x) variableInCalledTable[match(x, keyInCalledTable)])
)
}
getRelatedVariableWithTwoKeys <- function(firstKeyInCallingTable, secondKeyInCallingTable,
firstKeyInCalledTable, secondKeyInCalledTable, variableInCalledTable)
{
temp <- 1:length(firstKeyInCallingTable)
for (i in 1:length(temp))
{
temp[i] <- variableInCalledTable[firstKeyInCalledTable == firstKeyInCallingTable[i] &
secondKeyInCalledTable == secondKeyInCallingTable[i]]
}
return(temp)
}Define function to build table with Pearson’s correlations and linear regression models for a pair of dependent and independent variables per each crop in dataset:
linearModels.table <- function(dataset, dependentVariableName, independentVariableName)
{
tableTemp <- data.frame(
levels(dataset$crop),
numeric(nlevels(dataset$crop)),
numeric(nlevels(dataset$crop)),
numeric(nlevels(dataset$crop)),
numeric(nlevels(dataset$crop)),
numeric(nlevels(dataset$crop)),
numeric(nlevels(dataset$crop))
)
names(tableTemp) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(dataset$crop))
{
tableTemp$correlation_Pearson[tableTemp$crop == aCrop] <-
cor(x = dataset[dataset$crop == aCrop, dependentVariableName],
y = dataset[dataset$crop == aCrop, independentVariableName],
use = "pairwise.complete.obs")
linearModel <- lm(as.formula(paste0(dependentVariableName, "~", independentVariableName)), data = dataset[yieldData$crop == aCrop,])
tableTemp$intercept[tableTemp$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
tableTemp$intercept_p[tableTemp$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
tableTemp$speed[tableTemp$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
tableTemp$speed_p[tableTemp$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
tableTemp$adjrsquared[tableTemp$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop, linearModel)
return(tableTemp)
}Define auxiliar function for calculating cumulative curve (extracted from weatherModel.R; see documentation folder in Weather model):
getCumulativePrecipitationOfYear <- function(dailyPrecipitationYear)
{
return(cumsum(dailyPrecipitationYear) / sum(dailyPrecipitationYear))
}Global settings
Define the format of image files to be generated (only “png” and “eps”):
listOfImageFormats <- c("eps", "png")
# eps-specific arguments (cairo_ps)
pointSize = 12
fallbackResolution = 600Crop table
Load crop table:
cropTable <- read.csv("models/cropsTable.csv", skip = 4)Clean empty columns (generated with automatic names starting with X):
cropTable <- cropTable[,names(cropTable)[grep("^(?!X).*$", names(cropTable), perl = TRUE)]]Reorder crops for a more logical display:
WARNING: doing this manually makes the script sensitive to the value of crop-selection in NetLogo. Thus, this chunk must be updated if crop-selection is changed before running experiments.
cropTable$crop <- factor(cropTable$crop, levels = c("proso millet", "pearl millet", "rice", "barley", "wheat 1", "wheat 2"))Get vector of colours to represent crops:
cropColours <- qualitative_hcl(nlevels(cropTable$crop))
# alternative using only base R graphics
#cropColours <- rainbow(nlevels(cropTable$crop), s = 0.8, v = 0.8, end = 0.85)knitr::kable(cropTable)| crop | species | cultivar | T_sum | HI | I_50A | I_50B | T_base | T_opt | RUE | I_50maxH | I_50maxW | T_heat | T_ext | S_CO2 | S_water | sugSowingDay | sugHarvestingDay | root.zone.depth..mm. |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| wheat 1 | Triticum aestivum | Yecora Rojo | 2200 | 0.36 | 480.0 | 200.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| wheat 2 | Triticum aestivum | Batten | 2150 | 0.34 | 280.0 | 50.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| rice | Oryza sativa | IR72 | 2300 | 0.47 | 850.0 | 200.00 | 9 | 26 | 1.24 | 100 | 10 | 34 | 50 | 0.08 | 1.00 | 170 | 320 | 400 |
| barley | Hordeum vulgare | hypothetical; based on DSSAT v4.7 | 1762 | 0.42 | 350.0 | 170.00 | 0 | 15 | 1.24 | 100 | 20 | 34 | 45 | 0.08 | 0.40 | 330 | 100 | 1000 |
| pearl millet | Pennisetum glaucum | hypothetical; based on DSSAT v4.7 - MLCER047.CUL - Middle Variety | 1220 | 0.25 | 245.0 | 120.00 | 10 | 33 | 1.90 | 100 | 5 | 35 | 47 | 0.08 | 0.05 | 200 | 320 | 1000 |
| proso millet | Panicum miliaceum L. | hypothetical | 1328 | 0.29 | 157.3 | 96.75 | 0 | 18 | 3.00 | 100 | 5 | 34 | 45 | 0.08 | 0.05 | 107 | 220 | 1000 |
Experiment 0
Experiment 0
Experiment 0 simulates the growth and yield of crops included in ‘cropTable’ using the NetLogo implementation of SIMPLE crop model (Zhao et al. 2019), adapted as a module of the Indus Village model, and using daily weather variables imported from external file containing real-world data.
As the reference real-world data, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ between 01/01/1984 and 31/12/2007. The exact location used as a reference is:
- Rakhigarhi, Haryana, India (Latitude: 29.1687, Longitude: 76.0687)
We selected the ICASA Format’s parameters:
- Precipitation (PRECTOT)
- Wind speed at 2m (WS2M)
- Relative Humidity at 2m (RH2M)
- Dew/frost point at 2m (T2MDEW)
- Maximum temperature at 2m (T2M_MAX)
- Minimum temperature at 2m (T2M_MIN)
- All sky insolation incident on a horizontal surface (ALLSKY_SFC_SW_DWN)
- Temperature at 2m (T2M)
Load Soil Water (or Soil Water Balance) model code to calculate ARID given real-world data:
source("models/soilWaterModel/estimateETr.R")
source("models/soilWaterModel/watbal.model.R")Load and format the real-world dataset so that they can be used in the Soil Water model:
weatherData_rakhigarhi <- watbal.weather.file("models/weather_input/POWER_SinglePoint_Daily_19840101_20071231_029d17N_076d70E_5b401917.csv",
year = 1984:2007)Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):
weatherData_rakhigarhi$ETr <- estimateETr(weatherData_rakhigarhi$I,
weatherData_rakhigarhi$T2M,
weatherData_rakhigarhi$Tmax,
weatherData_rakhigarhi$Tmin,
weatherData_rakhigarhi$T2MDEW,
weatherData_rakhigarhi$WS2M)Load simulation data:
yieldData_exp0 <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=weather-input-data_experiment-name=exp0_initRandomSeed=0.csv")Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData_exp0$yield[is.na(yieldData_exp0$meanARID_grow)] <- NAForce the same order of crops used by cropTable:
yieldData_exp0$crop <- factor(yieldData_exp0$crop, levels = levels(cropTable$crop))Calculate the range of ARID and productivity (yield):
minARID_exp0 = round(min(c(yieldData_exp0$meanARID, yieldData_exp0$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID_exp0 = round(max(c(yieldData_exp0$meanARID, yieldData_exp0$meanARID_grow), na.rm = TRUE), digits = 2)
minYield_exp0 = round(min(yieldData_exp0$yield, na.rm = TRUE), digits = -1)
maxYield_exp0 = round(max(yieldData_exp0$yield, na.rm = TRUE), digits = -1)Replace precipitation_yearlyMean and precipitation_dailyCum_plateauValue_yearlyMean variables with values calculated from weatherData_rakhigarhi to compare experiment 0 with experiment 2.
This is a function to estimate the “plateau value” of a cumulative precipitation curve (one year), assuming it approximates a double logistic curve:
estimatePlateauValue <- function(dailyPrecipitation,
# cumulative curve slope with less than this is considered a plateau:
plateauMaxRate = 1/365, # daily precipitation with up to ~0.27% annualSum, which corresponds to a straight slope from 0 to 1 with no plateaus
# this is the percentage of annualSum
# which dismiss plateaus at the extremes of the cumulative curve:
toleranceOfProximityToExtremes = 0.005 # 0.5% of annualSum
)
{
annualSum = sum(dailyPrecipitation)
cumulativeCurve <- getCumulativePrecipitationOfYear(dailyPrecipitation)
# filter out values not in plateaus
daysInPlateaus <- dailyPrecipitation <= plateauMaxRate * annualSum
# filter out values too near the first and last plateaus
daysMoreLikelyToBeMiddlePlateau <-
cumulativeCurve > toleranceOfProximityToExtremes &
cumulativeCurve < 1 - toleranceOfProximityToExtremes
return(
getMode(
cumulativeCurve[
daysInPlateaus & daysMoreLikelyToBeMiddlePlateau
]
)
)
}
getMode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}yieldData_exp0$precipitation_yearlyMean <- rep(mean(yieldData_exp0$precipitation_yearTotal), nrow(yieldData_exp0))
yieldData_exp0$precipitation_dailyCum_plateauValue_yearlyMean <- mean(aggregate(
weatherData_rakhigarhi$RAIN,
by = list(weatherData_rakhigarhi$year),
FUN = estimatePlateauValue)[,2]
)Preliminary plotting of yield and ARID per crop:
ggplot(data = yieldData_exp0, aes(x = meanARID_grow, y = yield, color = crop, fill = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean ARID") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * 1) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
Preliminary plotting of yield and T_max per crop:
ggplot(data = yieldData_exp0, aes(x = meanT_max_grow, y = yield, color = crop, fill = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean T_max") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * 1) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
Experiment 0 results are plotted together with experiment 1, serving as points of comparison of the Crop model behaviour with and without the input of the Weather model.
Experiment 1
Experiment 1 is analogous to experiment 0, except that weather variables are simulated using the Weather model under the default parameter configuration, which is aimed to approximate the conditions in Rakhigarhi, Haryana, NW India, as reflected in the dataset used in experiment 0.
Loading and preparation
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NAForce the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))Calculate the range of ARID and productivity (yield):
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)Load daily simulated weather and ARID data:
weatherData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_weather_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")Calculate ARID in real-world data:
Use Soil Water model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.
The first step is to set all soil parameters, here assumed to be constant (the same values used in experiment 1 simulations):
watbalParameters <- watbal.define.param()Run the model inputing parameters and dataset:
watbalOutput_rakhigarhi <- watbal.model(watbalParameters,
weatherData_rakhigarhi,
typeOfParValue = "nominal")watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.
Show parameter values used here in Soil Water model (both in watbal.model and in the SIMPLE crop model). In experiemnt 1, they should be the same:
knitr::kable(
cbind(
parameters = c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF"),
SIMPLE = as.vector(unlist(yieldData[1, c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF")])),
watbal = c(200, # these values of elevation and albedo are the default in estimateETr
0.23,
as.vector(watbalParameters[1, c("CN", "DC", "z", "FC", "WHC", "WP", "MUF")]))
)
)| parameters | SIMPLE | watbal |
|---|---|---|
| elevation | 200 | 200 |
| albedo | 0.23 | 0.23 |
| CN | 65 | 65 |
| DC | 0.55 | 0.55 |
| z | 400 | 400 |
| FC | 0.21 | 0.21 |
| WHC | 0.15 | 0.15 |
| WP | 0.06 | 0.06 |
| MUF | 0.096 | 0.096 |
Summaries per day of year
Create data frames containing the daily summary statistics of simulations and real-world data:
- Create names vector (used for both datasets)
weatherSummaryNames <- c("dayOfYear",
"solarRadiation.mean",
"solarRadiation.sd",
"solarRadiation.max",
"solarRadiation.min",
"solarRadiation.error",
"temperature.mean",
"temperature.sd",
"temperature.max",
"temperature.min",
"temperature.error",
"maxTemperature.mean",
"maxTemperature.max",
"maxTemperature.min",
"maxTemperature.error",
"minTemperature.mean",
"minTemperature.max",
"minTemperature.min",
"minTemperature.error",
"temperature.lowerDeviation",
"temperature.lowerDeviation.error",
"temperature.upperDeviation",
"temperature.upperDeviation.error",
"precipitation.mean",
"precipitation.max",
"precipitation.min",
"precipitation.error",
"ARID.mean",
"ARID.max",
"ARID.min",
"ARID.error")- Simulation data
weatherSummary <- vector("list", length(weatherSummaryNames))
names(weatherSummary) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:365)
{
weatherSummary$dayOfYear <- c(weatherSummary$dayOfYear, day)
tempData_thisDay <- weatherData[weatherData$currentDayOfYear == day,]
# solar radiation
weatherSummary$solarRadiation.mean <- c(
weatherSummary$solarRadiation.mean,
mean(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.sd <- c(
weatherSummary$solarRadiation.sd,
sd(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.max <- c(
weatherSummary$solarRadiation.max,
max(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.min <- c(
weatherSummary$solarRadiation.min,
min(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.error <- c(
weatherSummary$solarRadiation.error,
qt(0.975,
length(tempData_thisDay$solarRadiation) - 1) *
sd(tempData_thisDay$solarRadiation, na.rm = TRUE) /
sqrt(length(tempData_thisDay$solarRadiation)))
# temperature
## daily mean
weatherSummary$temperature.mean <- c(
weatherSummary$temperature.mean,
mean(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.sd <- c(
weatherSummary$temperature.sd,
sd(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.max <- c(
weatherSummary$temperature.max,
max(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.min <- c(
weatherSummary$temperature.min,
min(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.error <- c(
weatherSummary$temperature.error,
qt(0.975,
length(tempData_thisDay$temperature) - 1) *
sd(tempData_thisDay$temperature, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature)))
## daily max
weatherSummary$maxTemperature.mean <- c(
weatherSummary$maxTemperature.mean,
mean(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.max <- c(
weatherSummary$maxTemperature.max,
max(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.min <- c(
weatherSummary$maxTemperature.min,
min(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.error <- c(
weatherSummary$maxTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
## daily min
weatherSummary$minTemperature.mean <- c(
weatherSummary$minTemperature.mean,
mean(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.max <- c(
weatherSummary$minTemperature.max,
max(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.min <- c(
weatherSummary$minTemperature.min,
min(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.error <- c(
weatherSummary$minTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature_min, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
## daily upper and lower deviation
weatherSummary$temperature.lowerDeviation <- c(
weatherSummary$temperature.lowerDeviation,
mean(tempData_thisDay$temperature -
tempData_thisDay$temperature_min))
weatherSummary$temperature.lowerDeviation.error <- c(
weatherSummary$temperature.lowerDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature -
tempData_thisDay$temperature_min,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
weatherSummary$temperature.upperDeviation <- c(
weatherSummary$temperature.upperDeviation,
mean(tempData_thisDay$temperature_max -
tempData_thisDay$temperature))
weatherSummary$temperature.upperDeviation.error <- c(
weatherSummary$temperature.upperDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max -
tempData_thisDay$temperature,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
# precipitation
weatherSummary$precipitation.mean <- c(
weatherSummary$precipitation.mean,
mean(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.max <- c(
weatherSummary$precipitation.max,
max(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.min <- c(
weatherSummary$precipitation.min,
min(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.error <- c(
weatherSummary$precipitation.error,
qt(0.975,
length(tempData_thisDay$RAIN) - 1) *
sd(tempData_thisDay$RAIN, na.rm = TRUE) /
sqrt(length(tempData_thisDay$RAIN)))
# ARID
weatherSummary$ARID.mean <- c(
weatherSummary$ARID.mean,
mean(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.max <- c(
weatherSummary$ARID.max,
max(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.min <- c(
weatherSummary$ARID.min,
min(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.error <- c(
weatherSummary$ARID.error,
qt(0.975,
length(tempData_thisDay$ARID) - 1) *
sd(tempData_thisDay$ARID, na.rm = TRUE) /
sqrt(length(tempData_thisDay$ARID)))
}
rm(day, tempData_thisDay)
weatherSummary <- data.frame(weatherSummary)- Real-world data
weatherSummary_rakhigarhi <- vector("list", length(weatherSummaryNames))
names(weatherSummary_rakhigarhi) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:366)
{
weatherSummary_rakhigarhi$dayOfYear <- c(weatherSummary_rakhigarhi$dayOfYear, day)
tempData <- weatherData_rakhigarhi[weatherData_rakhigarhi$day == day,]
tempWatbalData <- watbalOutput_rakhigarhi[watbalOutput_rakhigarhi$day == day,]
# solar radiation
weatherSummary_rakhigarhi$solarRadiation.mean <- c(
weatherSummary_rakhigarhi$solarRadiation.mean,
mean(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.sd <- c(
weatherSummary_rakhigarhi$solarRadiation.sd,
sd(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.max <- c(
weatherSummary_rakhigarhi$solarRadiation.max,
max(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.min <- c(
weatherSummary_rakhigarhi$solarRadiation.min,
min(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.error <- c(
weatherSummary_rakhigarhi$solarRadiation.error,
qt(0.975, length(tempData$I) - 1) *
sd(tempData$I, na.rm = TRUE) /
sqrt(length(tempData$I)))
# temperature
## daily mean
weatherSummary_rakhigarhi$temperature.mean <- c(
weatherSummary_rakhigarhi$temperature.mean,
mean(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.sd <- c(
weatherSummary_rakhigarhi$temperature.sd,
sd(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.max <- c(
weatherSummary_rakhigarhi$temperature.max,
max(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.min <- c(
weatherSummary_rakhigarhi$temperature.min,
min(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.error <- c(
weatherSummary_rakhigarhi$temperature.error,
qt(0.975, length(tempData$T2M) - 1) *
sd(tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$T2M)))
## daily max
weatherSummary_rakhigarhi$maxTemperature.mean <- c(
weatherSummary_rakhigarhi$maxTemperature.mean,
mean(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.max <- c(
weatherSummary_rakhigarhi$maxTemperature.max,
max(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.min <- c(
weatherSummary_rakhigarhi$maxTemperature.min,
min(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.error <- c(
weatherSummary_rakhigarhi$maxTemperature.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
## daily min
weatherSummary_rakhigarhi$minTemperature.mean <- c(
weatherSummary_rakhigarhi$minTemperature.mean,
mean(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.max <- c(
weatherSummary_rakhigarhi$minTemperature.max,
max(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.min <- c(
weatherSummary_rakhigarhi$minTemperature.min,
min(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.error <- c(
weatherSummary_rakhigarhi$minTemperature.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
## daily lower and upper deviation
weatherSummary_rakhigarhi$temperature.lowerDeviation <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation,
mean(tempData$T2M - tempData$Tmin))
weatherSummary_rakhigarhi$temperature.lowerDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$T2M - tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
weatherSummary_rakhigarhi$temperature.upperDeviation <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation,
mean(tempData$Tmax - tempData$T2M))
weatherSummary_rakhigarhi$temperature.upperDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax - tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
# precipitation
weatherSummary_rakhigarhi$precipitation.mean <- c(
weatherSummary_rakhigarhi$precipitation.mean,
mean(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.max <- c(
weatherSummary_rakhigarhi$precipitation.max,
max(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.min <- c(
weatherSummary_rakhigarhi$precipitation.min,
min(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.error <- c(
weatherSummary_rakhigarhi$precipitation.error,
qt(0.975, length(tempData$RAIN) - 1) *
sd(tempData$RAIN, na.rm = TRUE) /
sqrt(length(tempData$RAIN)))
# ARID
weatherSummary_rakhigarhi$ARID.mean <- c(
weatherSummary_rakhigarhi$ARID.mean,
mean(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.max <- c(
weatherSummary_rakhigarhi$ARID.max,
max(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.min <- c(
weatherSummary_rakhigarhi$ARID.min,
min(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.error <- c(
weatherSummary_rakhigarhi$ARID.error,
qt(0.975,
length(tempWatbalData$ARID) - 1) *
sd(tempWatbalData$ARID, na.rm = TRUE) /
sqrt(length(tempWatbalData$ARID)))
}
rm(day, tempData, tempWatbalData)
weatherSummary_rakhigarhi <- data.frame(weatherSummary_rakhigarhi)Get ranges of weather variables:
roundToMultiple <- function(i, baseOfMultiple, roundFunction = round)
{
return(match.fun(roundFunction)(i/baseOfMultiple) * baseOfMultiple)
}
rangeSolar = c(
roundToMultiple(min(c(weatherSummary$solarRadiation.min, weatherSummary_rakhigarhi$solarRadiation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$solarRadiation.max, weatherSummary_rakhigarhi$solarRadiation.max)), 5, ceiling))
rangeTemp = c(
roundToMultiple(min(c(weatherSummary$minTemperature.min, weatherSummary_rakhigarhi$minTemperature.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$maxTemperature.max, weatherSummary_rakhigarhi$maxTemperature.max)), 5, ceiling))
rangePrecip = c(
roundToMultiple(min(c(weatherSummary$precipitation.min, weatherSummary_rakhigarhi$precipitation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$precipitation.max, weatherSummary_rakhigarhi$precipitation.max)), 5, ceiling))Set colours for maximum and minimum temperature:
maxTemperatureColour = hsv(7.3/360, 74.6/100, 70/100)
minTemperatureColour = hsv(232/360, 64.6/100, 73/100)Fig - Weather variables
Create figure:
yearLengthInDays_sim = nlevels(factor(weatherSummary$dayOfYear))
yearLengthInDays_real = nlevels(factor(weatherSummary_rakhigarhi$dayOfYear))
yLabs <- c(expression(paste("solar radiation (", MJ/m^-2, ")")),
"temperature (C)", "precipitation (mm)", "ARID", "crop seasons")
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_weatherAndARID.png"
grScale = 2
fontRescale = 0
axisTextRescale = 0
marginTextRescale = 0
leftPlotMargin = 2
rightPlotMargin = 4
cumPrecipMargins <- c(8, 1.5)
textMarginWeight = 2
cropRangesHeadLength = 0.1
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 14
cropRangesLabelsSize = 0.95
png(plotName, width = grScale * 800, height = grScale * 1000)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_weatherAndARID.eps"
grScale = 1.2
fontRescale = 0.1
axisTextRescale = -0.1
marginTextRescale = -0.5
leftPlotMargin = 1.5
rightPlotMargin = 3.8
cumPrecipMargins <- c(6, 1.5)
textMarginWeight = 1
cropRangesHeadLength = 0.08
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 10
cropRangesLabelsSize = 0.9
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = grScale * 8,
height = grScale * 10,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
nColumns = 3
nRowsExceptTopAndBottom = 5
layout(rbind(1:nColumns, # top row (headers)
# main block with plots
matrix((nColumns + 1):(nColumns + (nColumns * nRowsExceptTopAndBottom)),
nrow = nRowsExceptTopAndBottom,
ncol = nColumns,
byrow = FALSE),
# bottom row (x axis title)
c((nColumns + 1) + (nColumns * nRowsExceptTopAndBottom),
rep(((nColumns + 1) +(nColumns * nRowsExceptTopAndBottom) + 1), nColumns - 1))),
widths = c(1, 10, 12),
heights = c(textMarginWeight,
rep(10, nRowsExceptTopAndBottom) + c(1, rep(0, nRowsExceptTopAndBottom - 2), -2),
textMarginWeight)
)
# top row: empty and two headers
par(mar = c(0, 0, 0, 0))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
par(mar = c(0, leftPlotMargin, 0, 0))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (2 + fontRescale),
labels = "Rakhigarhi (1984-2007)")
par(mar = c(0, 0, 0, rightPlotMargin))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.45, y = 0.5, font = 4,
cex = grScale * (2 + fontRescale),
labels = "Experiment 1 (750 years)")
# First column: y axis titles
par(cex = grScale,
cex.axis = grScale * (0.8 + axisTextRescale))
par(mar = c(0, 0, 0, 0.4))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[1])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.4, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[3])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[4])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[5])
#==============================================================================
# Second column: real-world data
# 1. solar radiation
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean + weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean - weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.max,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.min,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeSolar[1], rangeSolar[2], 5))
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean + weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean - weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.max,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.min,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean + weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean - weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.max,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.min,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean + weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean - weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.max,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.min,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeTemp[1], rangeTemp[2], 5))
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], leftPlotMargin, 0.1, 0.1))
plot(c(1, yearLengthInDays_real), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (year in unique(weatherData_rakhigarhi$year))
{
lines(1:nrow(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year,]),
getCumulativePrecipitationOfYear(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
rm(year)
# daily values
par(new = T,
mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean + weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean - weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.max,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.min,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangePrecip[1], rangePrecip[2], 10))
# 4. ARID
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean + weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean - weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.max,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.min,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(0, 1, 0.1))
#5. crop growing seasons
par(mar = c(3, leftPlotMargin, 0.1, 0.1))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = 365, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
#==============================================================================
# Third column: simulated data
# 1. solar radiation
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean + weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean - weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.max,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.min,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_sim,
weatherSummary$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean + weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean - weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.max,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.min,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_sim,
weatherSummary$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean + weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean - weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.max,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.min,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_sim,
weatherSummary$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean + weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean - weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.max,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.min,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], 0.1, 0.1, rightPlotMargin))
plot(c(1, yearLengthInDays_sim), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (randomSeed in unique(weatherData$randomSeed))
{
for (year in unique(weatherData$currentYear))
{
lines(1:nrow(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year,]),
getCumulativePrecipitationOfYear(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
}
rm(year, randomSeed)
axis(4, at = seq(0, 1, 0.25))
mtext("cumulative annual sum", 4, line = 2,
cex = grScale * (1.5 + marginTextRescale))
# daily values
par(new = T,
mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean + weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean - weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.max,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.min,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 4. ARID
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean + weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean - weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.max,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.min,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
#5. crop growing seasons
par(mar = c(3, 0.1, 0.1, rightPlotMargin))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = yearLengthInDays_sim, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
mtext(aCrop,
side = 4, col = cropColours[aCropIndex],
padj = ((1 - aCropIndex/nrow(cropTable)) - 0.4) * cropRangesLabelsVerticalSpread,
las = 2, cex = grScale * cropRangesLabelsSize)
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
# bottom row: empty and "day of year" or x axis title
par(mar = c(0, 0, 0, 0))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.45, y = 0.5, font = 4,
cex = grScale * (0.8 + fontRescale),
labels = "day of year")
dev.off()
}
rm(imageFormat, grScale, fontRescale, yLabs, nColumns, nRowsExceptTopAndBottom,
axisTextRescale, marginTextRescale, leftPlotMargin, rightPlotMargin, cumPrecipMargins, textMarginWeight,
tempData, aCrop, aCropIndex, cropRangesHeadLength, cropRangesLineWidth)knitr::include_graphics(plotName)
Display weather parameters:
knitr::kable(t(yieldData[1, 2:26]), col.names = "default (NW India echo)")| default (NW India echo) | |
|---|---|
| temperature_annualMaxAt2m | 37.00 |
| temperature_annualMinAt2m | 12.80 |
| temperature_meanDailyFluctuation | 2.20 |
| temperature_dailyLowerDeviation | 6.80 |
| temperature_dailyUpperDeviation | 7.90 |
| solar_annualMax | 24.20 |
| solar_annualMin | 9.20 |
| solar_meanDailyFluctuation | 3.30 |
| CO2_annualMin | 245.00 |
| CO2_annualMax | 255.00 |
| CO2_meanDailyFluctuation | 1.00 |
| precipitation_yearlyMean | 489.00 |
| precipitation_yearlySd | 142.20 |
| precipitation_dailyCum_nSamples | 200.00 |
| precipitation_dailyCum_maxSampleSize | 10.00 |
| precipitation_dailyCum_plateauValue_yearlyMean | 0.25 |
| precipitation_dailyCum_plateauValue_yearlySd | 0.10 |
| precipitation_dailyCum_inflection1_yearlyMean | 40.00 |
| precipitation_dailyCum_inflection1_yearlySd | 5.00 |
| precipitation_dailyCum_rate1_yearlyMean | 0.07 |
| precipitation_dailyCum_rate1_yearlySd | 0.02 |
| precipitation_dailyCum_inflection2_yearlyMean | 240.00 |
| precipitation_dailyCum_inflection2_yearlySd | 20.00 |
| precipitation_dailyCum_rate2_yearlyMean | 0.08 |
| precipitation_dailyCum_rate2_yearlySd | 0.02 |
Fig - productivity (yield) per crop type (with default parameter setting)
Summary statistics per crop:
yieldData_summary <- tapply(as.numeric(as.character(yieldData$yield)), yieldData$crop, summary)
yieldData_summary <- data.frame(Reduce(rbind, yieldData_summary), row.names = names(yieldData_summary))## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)knitr::kable(yieldData_summary)| Min. | X1st.Qu. | Median | Mean | X3rd.Qu. | Max. | |
|---|---|---|---|---|---|---|
| proso millet | 886.0298 | 1059.6563 | 1110.0890 | 1108.8338 | 1157.3545 | 1373.5597 |
| pearl millet | 300.5779 | 332.6118 | 342.6917 | 342.9984 | 352.8890 | 392.0383 |
| rice | 0.0000 | 230.7686 | 265.1131 | 261.6839 | 299.6243 | 397.6978 |
| barley | 337.8926 | 389.9014 | 412.2040 | 414.7265 | 437.7929 | 520.7963 |
| wheat 1 | 280.4535 | 327.1036 | 345.8900 | 347.8794 | 367.3557 | 437.3639 |
| wheat 2 | 296.1396 | 338.7688 | 357.1001 | 359.2555 | 378.5212 | 447.0028 |
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldPerCrop.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldPerCrop.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
par(mar = c(6,6,1,1), cex.lab = grScale, cex.axis = 0.8 * grScale)
boxplot(yield ~ factor(crop), data = yieldData,
ylab = expression(paste("yield (", g/m^2, ")")), xlab = "",
col = cropColours)
# boxplot(yield ~ factor(crop), data = yieldData_exp0,
# ylab = expression(paste("yield (", g/m^2, ")")), xlab = "",
# col = "red")
stripchart(yield ~ factor(crop), data = yieldData_exp0,
method = "jitter",
col = "red", pch = 4,
vertical = TRUE, add = TRUE)
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldPerCrop_ggplot.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldPerCrop_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
ggplot(yieldData, aes(y = yield, x = crop, fill = crop)) +
geom_violin(position = "dodge") +
geom_jitter(data = yieldData_exp0, aes(y = yield, x = crop, fill = crop), color = "red", shape = 4, width = 0.1) +
#coord_flip() +
#geom_boxplot() +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
ylab(expression(paste("yield (", g/m^2, ")")))
})
dev.off()
}## Warning: Removed 75 rows containing non-finite values (stat_ydensity).## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 75 rows containing non-finite values (stat_ydensity).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Fig - mean daily ARID during each crop growing season versus mean ARID of the harvest calendar year
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_ARIDVsARID_grow.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_ARIDVsARID_grow.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
layout(matrix(c(1, 2), nrow = 1, ncol = 2),
widths = c(10, 3))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
boxplot(meanARID_grow ~ factor(crop), data = yieldData,#[yieldData$yield > 0,], # show only non-zero yield
ylab = "mean ARID", xlab = "growing seasons",
#las = 2,
col = cropColours)
stripchart(meanARID_grow ~ factor(crop), data = yieldData_exp0,
method = "jitter",
col = "red", pch = 4,
vertical = TRUE, add = TRUE)
text(x = 0.5, y = minARID * 1.005, labels = "a", font = 2, cex = grScale)
par(mar = c(5,1,1,1))
xhist <- hist(yieldData$meanARID, breaks = 20, plot = FALSE)
barplot(xhist$counts, axes = FALSE, space = 0, horiz = TRUE,
xlab = "current year")
text(x = max(xhist$counts) * 0.9, y = minARID * 1.005, labels = "b", font = 2, cex = grScale)
xhist <- hist(yieldData_exp0$meanARID, breaks = 20, plot = FALSE)
barplot(xhist$counts, axes = FALSE, space = 0, horiz = TRUE, add = TRUE, col = "red")
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_ARIDVsARID_grow_ggplot.png"
grScale = 2
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_ARIDVsARID_grow_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 6 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
# build temporary data frames where year meanARID is coded as an extra crop type
tempData <- yieldData[, c("crop", "meanARID_grow")]
tempData$crop <- sapply(tempData$crop, as.character)
for (year in unique(yieldData$currentYear))
{
for (randomSeed in unique(yieldData$randomSeed))
{
tempData <- rbind(tempData,
cbind(crop = "year",
meanARID_grow = yieldData[yieldData$randomSeed == randomSeed & yieldData$currentYear == year,]$meanARID[1]))
}
}
tempData$meanARID_grow <- as.numeric(tempData$meanARID_grow)
tempData$crop <- factor(tempData$crop, levels = c(levels(cropTable$crop), "year"))
tempData_exp0 <- yieldData_exp0[, cbind("crop", "meanARID_grow")]
tempData_exp0$crop <- sapply(tempData_exp0$crop, as.character)
for (year in unique(yieldData_exp0$currentYear))
{
tempData_exp0 <- rbind(tempData_exp0, c(crop = "year",
meanARID_grow = yieldData_exp0[yieldData_exp0$currentYear == year,]$meanARID[1]))
}
tempData_exp0$meanARID_grow <- as.numeric(tempData_exp0$meanARID_grow)
tempData_exp0$crop <- factor(tempData_exp0$crop, levels = c(levels(cropTable$crop), "year"))
print({
(ggplot(tempData, aes(fill = crop, y = meanARID_grow, x = crop)) +
geom_violin(position = "dodge") +
geom_jitter(data = tempData_exp0, aes(y = meanARID_grow, x = crop, fill = crop), color = "red", shape = 4, width = 0.1) +
#geom_boxplot() +
geom_hline(yintercept = mean(tempData$meanARID_grow[tempData$crop == "year"]), size = grScale) +
geom_hline(yintercept = mean(tempData_exp0$meanARID_grow[tempData_exp0$crop == "year"]), size = 0.5 * grScale, color = "red") +
scale_fill_manual(values = c(cropColours, "grey"), name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
theme(axis.title.x = element_blank()) +
ylab("mean ARID")) +
ylim(minARID, maxARID)
})
dev.off()
}## Warning: Removed 76 rows containing non-finite values (stat_ydensity).## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 76 rows containing non-finite values (stat_ydensity).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Use One-Way ANOVA to express the differences between crop and year-round means (simulation results):
NOTE: using the tempDatacreated in the last code chunk.
anovaTable <- data.frame(
cbind(
Count = lapply(levels(tempData$crop), function(x) { nrow(tempData[tempData$crop == x, ]) }),
Mean = lapply(levels(tempData$crop), function(x) { mean(tempData$meanARID_grow[tempData$crop == x], na.rm = TRUE) }),
Variance = lapply(levels(tempData$crop), function(x) { var(tempData$meanARID_grow[tempData$crop == x], na.rm = TRUE) })
),
row.names = levels(tempData$crop)
)
knitr::kable(anovaTable)| Count | Mean | Variance | |
|---|---|---|---|
| proso millet | 750 | 0.9144407 | 0.003808305 |
| pearl millet | 750 | 0.6223765 | 0.01008866 |
| rice | 750 | 0.6776866 | 0.005708732 |
| barley | 750 | 0.7854098 | 0.01043134 |
| wheat 1 | 750 | 0.7909198 | 0.009883204 |
| wheat 2 | 750 | 0.7909198 | 0.009883204 |
| year | 750 | 0.7762129 | 0.003579471 |
anova_results <- aov(meanARID_grow ~ crop, data = tempData)
pander::pander(summary(anova_results),
caption = "One-Way ANOVA (dep: ARID_grow, indep: crop)")| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| crop | 6 | 39.09 | 6.516 | 858.4 | 0 |
| Residuals | 5168 | 39.23 | 0.007591 | NA | NA |
Post hoc test (Holm’s method):
pairwise_ttest_results <- pairwise.t.test(tempData$meanARID_grow, tempData$crop, p.adjust.method = "holm")
pander::pander(pairwise_ttest_results$p.value,
caption = "Post hoc pairwise t-test (Holm's method)")| proso millet | pearl millet | rice | barley | |
|---|---|---|---|---|
| pearl millet | 0 | NA | NA | NA |
| rice | 0 | 2.058e-33 | NA | NA |
| barley | 3.168e-164 | 2.563e-251 | 2.239e-117 | NA |
| wheat 1 | 1.766e-151 | 1.362e-266 | 7.951e-129 | 0.6858 |
| wheat 2 | 1.766e-151 | 1.362e-266 | 7.951e-129 | 0.6858 |
| year | 2.33e-189 | 5.611e-230 | 7.902e-101 | 0.1709 |
| wheat 1 | wheat 2 | |
|---|---|---|
| pearl millet | NA | NA |
| rice | NA | NA |
| barley | NA | NA |
| wheat 1 | NA | NA |
| wheat 2 | 1 | NA |
| year | 0.007191 | 0.007191 |
Tukey Honest Significant Differences:
pander::pander(TukeyHSD(anova_results, "crop")$crop,
caption = "Tukey Honest Significant Differences")| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| pearl millet-proso millet | -0.2921 | -0.3053 | -0.2788 | 0 |
| rice-proso millet | -0.2368 | -0.25 | -0.2235 | 0 |
| barley-proso millet | -0.129 | -0.1424 | -0.1156 | 0 |
| wheat 1-proso millet | -0.1235 | -0.1369 | -0.1101 | 0 |
| wheat 2-proso millet | -0.1235 | -0.1369 | -0.1101 | 0 |
| year-proso millet | -0.1382 | -0.1515 | -0.125 | 0 |
| rice-pearl millet | 0.05531 | 0.04204 | 0.06858 | 0 |
| barley-pearl millet | 0.163 | 0.1496 | 0.1764 | 0 |
| wheat 1-pearl millet | 0.1685 | 0.1552 | 0.1819 | 0 |
| wheat 2-pearl millet | 0.1685 | 0.1552 | 0.1819 | 0 |
| year-pearl millet | 0.1538 | 0.1406 | 0.1671 | 0 |
| barley-rice | 0.1077 | 0.09434 | 0.1211 | 0 |
| wheat 1-rice | 0.1132 | 0.09985 | 0.1266 | 0 |
| wheat 2-rice | 0.1132 | 0.09985 | 0.1266 | 0 |
| year-rice | 0.09853 | 0.08526 | 0.1118 | 0 |
| wheat 1-barley | 0.00551 | -0.007987 | 0.01901 | 0.8928 |
| wheat 2-barley | 0.00551 | -0.007987 | 0.01901 | 0.8928 |
| year-barley | -0.009197 | -0.02258 | 0.004187 | 0.3977 |
| wheat 2-wheat 1 | -2.887e-15 | -0.0135 | 0.0135 | 1 |
| year-wheat 1 | -0.01471 | -0.02809 | -0.001323 | 0.0205 |
| year-wheat 2 | -0.01471 | -0.02809 | -0.001323 | 0.0205 |
rm(tempData, tempData_exp0, anova_results, pairwise_ttest_results, anovaTable)Fig - mean ARID vs productivity (yield)
Plotting absolute values of yield vs ARID:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldVsARID.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldVsARID.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (harvest year)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "a", font = 2, cex = grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (grow season)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID_grow[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$S_water[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "S_water",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale)knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldVsARID_ggplot.png"
grScale = 2
png(plotName, width = grScale * 350, height = grScale * 500)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldVsARID_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 3.5 * grScale,
height = 5 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
ggplot(data = yieldData, aes(x = meanARID_grow, y = yield, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', se = FALSE) +
geom_point(data = yieldData_exp0, aes(x = meanARID_grow, y = yield), size = 2, shape = 4, color = 'red') +
geom_smooth(data = yieldData_exp0, aes(x = meanARID_grow, y = yield), method = 'lm', color = 'red', se = FALSE, linetype = "dashed") +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean ARID") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
knitr::kable(linearModels.table(yieldData, "yield", "meanARID_grow"), caption = "Pearson's correlation and linear model between ARID (x) and Yield (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.0202109 | 1087.0279 | 0 | 23.84613 | 0.5805157 | -0.0009279 |
| pearl millet | -0.0482022 | 347.6334 | 0 | -7.44716 | 0.1872903 | 0.0009897 |
| rice | -0.9094099 | 693.8981 | 0 | -637.77881 | 0.0000000 | 0.8267950 |
| barley | -0.8629164 | 629.6293 | 0 | -273.61862 | 0.0000000 | 0.7442715 |
| wheat 1 | -0.8470997 | 531.3303 | 0 | -231.94635 | 0.0000000 | 0.7171872 |
| wheat 2 | -0.8582940 | 540.0265 | 0 | -228.55794 | 0.0000000 | 0.7363044 |
Fig - mean T_max vs productivity (yield)
Plotting absolute values of yield vs ARID:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldVsT_max_ggplot.png"
grScale = 2
png(plotName, width = grScale * 350, height = grScale * 500)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldVsT_max_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 3.5 * grScale,
height = 5 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
ggplot(data = yieldData, aes(x = meanT_max_grow, y = yield, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', se = FALSE) +
geom_point(data = yieldData_exp0, aes(x = meanT_max_grow, y = yield), size = 2, shape = 4, color = 'red') +
geom_smooth(data = yieldData_exp0, aes(x = meanT_max_grow, y = yield), method = 'lm', color = 'red', se = FALSE, linetype = "dashed") +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean T_max") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness (simulated data):
knitr::kable(linearModels.table(yieldData, "yield", "meanT_max_grow"), caption = "Pearson's correlation and linear model between T_max (x) and Yield (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | -0.4051720 | 7323.5024 | 0.0000000 | -144.657750 | 0.0000000 | 0.1630469 |
| pearl millet | 0.2197011 | -222.1159 | 0.0157229 | 16.935428 | 0.0000000 | 0.0469962 |
| rice | -0.0183107 | 455.9481 | 0.2401432 | -5.458571 | 0.6166068 | -0.0010012 |
| barley | -0.0921359 | 843.4203 | 0.0000012 | -16.622115 | 0.0130708 | 0.0071176 |
| wheat 1 | -0.0742136 | 648.4933 | 0.0000181 | -11.473209 | 0.0457633 | 0.0041322 |
| wheat 2 | -0.0774328 | 664.2962 | 0.0000064 | -11.642162 | 0.0371170 | 0.0046210 |
Integrated plot with yield vs mean ARID during growing season and daily maximum temperature:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp1_yieldvsARIDAndTmax.png"
grScale = 2
fontRescale = 1
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp1_yieldvsARIDAndTmax.eps"
grScale = 2
fontRescale = 0.8
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 3.5 * grScale,
height = 5 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
(ggplot(data = yieldData, aes(y = yield, x = meanARID_grow, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', se = FALSE) +
geom_point(data = yieldData_exp0, aes(x = meanARID_grow, y = yield), size = 2, shape = 4, color = 'red') +
geom_smooth(data = yieldData_exp0, aes(x = meanARID_grow, y = yield), method = 'lm', color = 'red', se = FALSE, linetype = "dashed") +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab(expression(paste("yield (", g/m^2, ")"))) +
xlab("mean ARID") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = yield, x = meanT_max_grow, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', se = FALSE) +
geom_point(data = yieldData_exp0, aes(x = meanT_max_grow, y = yield), size = 2, shape = 4, color = 'red') +
geom_smooth(data = yieldData_exp0, aes(x = meanT_max_grow, y = yield), method = 'lm', color = 'red', se = FALSE, linetype = "dashed") +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab("") +
xlab("mean T_max") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 1, widths = c(1, 1))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Apply random forest to identify ranking of parameters and indicators that best explain variation in yield.
Name parameters to consider:
# those that make sense when considering all crops (includes crop parameters)
pars_overall <- c(
"T_sum", "HI", "I_50A", "I_50B", "T_base", "T_opt", "RUE",
"I_50maxH", "I_50maxW", "T_heat", "T_extreme", "S_water",
"sowingDay", "harvestDay", "meanARID_grow", "meanT_max_grow"
)
# those that make sense to analise variance within each crop (weather-related indicators)
pars_perCrop <- c(
"meanARID_grow", "meanT_max_grow"
)Overall, across crops:
if(!file.exists('randomForestOutputData/randomForest_yield.RData'))
{
randomForest_yield <-
randomForest::randomForest(
as.formula(paste0("yield ~ ", paste(pars_overall, collapse = '+'))),
data = yieldData[complete.cases(yieldData),],
mtry = 5,
ntree = 1000,
importance = T
)
# save
save(randomForest_yield, file = 'randomForestOutputData/randomForest_yield.RData')
}
# load
load('randomForestOutputData/randomForest_yield.RData')Variable importance plot:
randomForest::varImpPlot(randomForest_yield)
NOTE: sowingDayis so high in ranking only because proso millet, with its distinct season, returns substantially higher yields than all other crops.
Per each crop:
fileName = paste0('randomForestOutputData/randomForest_yield_perCrop.RData')
if(!file.exists(fileName))
{
for (aCrop in cropTable$crop)
{
randomForest_yield[[aCrop]] <-
randomForest::randomForest(
as.formula(paste0("yield ~ ", paste(pars_perCrop, collapse = '+'))),
data = yieldData[yieldData$crop == aCrop & complete.cases(yieldData),],
mtry = 4,
ntree = 1000,
importance = T
)
}
# save
save(randomForest_yield, file = fileName)
}
# load
load(fileName)Variable importance plot:
randomForest::varImpPlot(randomForest_yield[['proso millet']])
randomForest::varImpPlot(randomForest_yield[['pearl millet']])
randomForest::varImpPlot(randomForest_yield[['rice']])
randomForest::varImpPlot(randomForest_yield[['barley']])
randomForest::varImpPlot(randomForest_yield[['wheat 1']])
randomForest::varImpPlot(randomForest_yield[['wheat 2']])
Clear temporary objects:
rm(watbalOutput_rakhigarhi, watbalParameters,
weatherData, weatherData_rakhigarhi, weatherSummary, weatherSummary_rakhigarhi, weatherSummaryNames,
xhist, yieldData, yieldData_summary,
yearLengthInDays_real, yearLengthInDays_sim,
maxARID, minARID, maxYield, minYield, rangePrecip, rangeSolar, rangeTemp,
maxTemperatureColour, minTemperatureColour,
plotName)Experiment 2
Experiment 2 perform the same runs of experiment 1, except by varying stochastically two weather model parameters concerning precipitation: precipitation_yearlyMean (average annual sum of precipitation) and precipitation_dailyCum_plateauValue_yearlyMean (average plateau value of the cumulative curve of daily precipitation).
Loading and preparation
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=precipitation-variation_experiment-name=exp2_initRandomSeed=0.csv")Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NAForce the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))Calculate the range of ARID, yield, precipitation annual sum (yearly mean) and the “plateau value” (yearly mean) of the cumulative curve of daily precipitation:
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)
minPrecipitation_yearlyMean = round(min(yieldData$precipitation_yearlyMean), digits = -2)
maxPrecipitation_yearlyMean = round(max(yieldData$precipitation_yearlyMean), digits = -2)
minPrecipitation_dailyCum_plateauValue_yearlyMean = round(min(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
maxPrecipitation_dailyCum_plateauValue_yearlyMean = round(max(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)Display explored values of precipitation_yearlyMean and precipitation_dailyCum_plateauValue_yearlyMean:
layout(matrix(1:2, ncol = 2))
hist(yieldData$precipitation_yearlyMean, main = "precipitation_yearlyMean")
hist(yieldData$precipitation_dailyCum_plateauValue_yearlyMean, main = "precipitation_dailyCum_plateauValue_yearlyMean")
Fig - Effect of precipitation parametres on ARID
Plotting ARID vs precipitation yearly mean and seasonal balance (plateau value):
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp2_ARIDvsPrecipitation.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp2_ARIDvsPrecipitation.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minPrecipitation_yearlyMean, maxPrecipitation_yearlyMean),
c(minARID, maxARID),
xlab = "precipitation annual sum (mm) (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_yearlyMean * 0.95, y = minARID * 1.2, labels = "a", font = 2, cex = grScale)
plot(c(minPrecipitation_dailyCum_plateauValue_yearlyMean, maxPrecipitation_dailyCum_plateauValue_yearlyMean),
c(minARID, maxARID),
xlab = "cumulative daily precipitation plateau value (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_dailyCum_plateauValue_yearlyMean * 0.95, y = minARID * 1.2, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$sugSowingDay[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "sowing day of year",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale, cex.points)knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp2_ARIDvsPrecipitation_ggplot.png"
grScale = 2
fontRescale = 1
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp2_ARIDvsPrecipitation_ggplot.eps"
grScale = 2
fontRescale = 0.8
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_point(data = yieldData_exp0, aes(y = meanARID_grow, x = precipitation_yearlyMean), color = "red", shape = 4) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("mean ARID") +
xlab("precipitation annual sum (mm) (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_point(data = yieldData_exp0, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean), color = "red", shape = 4) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("") +
xlab("cumulative daily precipitation plateau value (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 1, widths = c(1, 1))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness (simulated data):
knitr::kable(linearModels.table(yieldData, "meanARID_grow", "precipitation_yearlyMean"), caption = "Pearson's correlation and linear model between precipitation annual sum yearly mean (x) and ARID (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | -0.5760686 | 0.9960945 | 0 | -0.0001314 | 0 | 0.3309618 |
| pearl millet | -0.7070320 | 0.9428775 | 0 | -0.0003964 | 0 | 0.4992257 |
| rice | -0.7308760 | 0.9545028 | 0 | -0.0003454 | 0 | 0.5335569 |
| barley | -0.7133386 | 0.8883206 | 0 | -0.0004509 | 0 | 0.5081726 |
| wheat 1 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
| wheat 2 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
knitr::kable(linearModels.table(yieldData, "meanARID_grow", "precipitation_dailyCum_plateauValue_yearlyMean"), caption = "Pearson's correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and ARID (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.1633846 | 0.8971123 | 0 | 0.0628231 | 6.9e-06 | 0.0253933 |
| pearl millet | 0.4151101 | 0.5367639 | 0 | 0.3922293 | 0.0e+00 | 0.1712098 |
| rice | 0.4268577 | 0.6015239 | 0 | 0.3400367 | 0.0e+00 | 0.1811142 |
| barley | -0.3736394 | 0.8743004 | 0 | -0.3980747 | 0.0e+00 | 0.1384163 |
| wheat 1 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
| wheat 2 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
Fig - Effect of precipitation parametres on yield
Plotting yield vs precipitation yearly mean and seasonal balance (plateau value):
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp2_yieldvsPrecipitation_ggplot.png"
grScale = 2
fontRescale = 1
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/exp2_yieldvsPrecipitation_ggplot.eps"
grScale = 2
fontRescale = 0.8
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 5 * grScale,
height = 3 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
(ggplot(data = yieldData, aes(y = yield, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_point(data = yieldData_exp0, aes(y = yield, x = precipitation_yearlyMean), color = "red", shape = 4) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab(expression(paste("yield (", g/m^2, ")"))) +
xlab("precipitation annual sum (mm) (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = yield, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_point(data = yieldData_exp0, aes(y = yield, x = precipitation_dailyCum_plateauValue_yearlyMean), color = "red", shape = 4) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab("") +
xlab("cumulative daily precipitation plateau value (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 1, widths = c(1, 1))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## Warning: Removed 3 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
knitr::kable(linearModels.table(yieldData, "yield", "precipitation_yearlyMean"), caption = "Pearson's correlation and linear model between precipitation annual sum yearly mean (x) and yield (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | -0.0622356 | 1117.31608 | 0 | -0.0158849 | 0.0885303 | 0.0025416 |
| pearl millet | -0.0461655 | 343.26635 | 0 | -0.0024760 | 0.2066404 | 0.0007972 |
| rice | 0.7255094 | 56.17851 | 0 | 0.2541253 | 0.0000000 | 0.5257308 |
| barley | 0.6671483 | 389.71758 | 0 | 0.1223045 | 0.0000000 | 0.4443194 |
| wheat 1 | 0.6650579 | 327.07130 | 0 | 0.1017196 | 0.0000000 | 0.4415307 |
| wheat 2 | 0.6704850 | 338.27122 | 0 | 0.1007498 | 0.0000000 | 0.4487888 |
knitr::kable(linearModels.table(yieldData, "yield", "precipitation_dailyCum_plateauValue_yearlyMean"), caption = "Pearson's correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and yield (y) per crop")| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.0330008 | 1101.8486 | 0 | 14.196494 | 0.3667893 | -0.0002464 |
| pearl millet | 0.0188182 | 341.1269 | 0 | 1.701073 | 0.6068719 | -0.0009823 |
| rice | -0.4003169 | 308.5154 | 0 | -236.329820 | 0.0000000 | 0.1591310 |
| barley | 0.3789046 | 388.6906 | 0 | 117.073782 | 0.0000000 | 0.1423841 |
| wheat 1 | 0.3783210 | 326.1346 | 0 | 97.524925 | 0.0000000 | 0.1419416 |
| wheat 2 | 0.3755728 | 338.1276 | 0 | 95.117190 | 0.0000000 | 0.1398669 |
Fig - Effects of precipitation parametres on yield and ARID
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/exp2_yieldAndARIDvsPrecipitation_ggplot.png"
grScale = 2
fontRescale = 1
png(plotName, width = grScale * 600, height = grScale * 600)
}
if (imageFormat == "eps")
{
plotName = "figures/exp2_yieldAndARIDvsPrecipitation_ggplot.eps"
grScale = 2
fontRescale = 0.8
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
width = 6 * grScale,
height = 6 * grScale,
pointsize = pointSize,
fallback_resolution = fallbackResolution
)
}
print({
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("mean ARID") +
xlab("") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("") +
xlab("") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = yield, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab(expression(paste("yield (", g/m^2, ")"))) +
xlab("precipitation annual sum (mm) (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = yield, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1, alpha = 0.3) +
geom_smooth(method = 'lm', color = 'black', fill = 'black', se = FALSE) +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
ylab("") +
xlab("cumulative daily precipitation plateau value (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15 * fontRescale),
axis.title = element_text(size = 15 * fontRescale),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 2, ncol = 2, widths = c(1, 1), heights = c(1, 1))
})
dev.off()
}## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).rm(imageFormat, grScale)knitr::include_graphics(plotName)
Clear temporary objects:
rm(yieldData,
maxARID, minARID, maxYield, minYield,
maxPrecipitation_yearlyMean, minPrecipitation_yearlyMean,
maxPrecipitation_dailyCum_plateauValue_yearlyMean, minPrecipitation_dailyCum_plateauValue_yearlyMean,
plotName)Fig - terrains generated with the Land model
Loading and preparation
Defining terrain random number generator seed (must be available inside “terrains” folder:
listOfTerrainSeeds <- c("0", "35", "56", "72", "92")Setting dimensions:
figScale = 2
terrainDim = 450 * figScale
topHeight = 50 * figScale
bottomHeight = 150 * figScale
leftWidth = 150 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 5 + leftWidth + 4 * internalMargin
figHeight = terrainDim * length(listOfTerrainSeeds) + topHeight + bottomHeight + (length(listOfTerrainSeeds) - 1) * internalMarginCreate separate component images
Create margin header images:
headers <- c("seed",
"Elevation",
"Soil texture",
"Soil texture type",
"Ecological composition",
"Cover type")
namesInFile <- c("seed", "elev", "soilText", "soilTextType", "ecolComp", "coverType")
for (i in 1:length(headers))
{
#if (i == 1) { width = leftWidth } else { width = terrainDim }
width = ifelse(i == 1, leftWidth, terrainDim)
png(paste0("FigTerrains_tempFiles/FigTerrains_", namesInFile[i], ".PNG"),
width = width,
height = topHeight)
par(mar = rep(0, 4))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.7, cex = figScale * 3.7, adj = c(0.5, 0.5),
labels = headers[i])
dev.off()
}
rm(i)
for (seed in listOfTerrainSeeds)
{
png(paste0("FigTerrains_tempFiles/FigTerrains_seed=", seed, ".PNG"),
width = leftWidth,
height = terrainDim)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.5, cex = figScale * 5, adj = .5,
labels = seed)
dev.off()
}
rm(seed)Create legend images:
# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 10)
png("FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = c(2,2,1,2), cex = figScale * 2)
barplot(rep(100, 10),
col = rgb(elev - (100/255), elev, 0),
border = rgb(elev - (100/255), elev, 0),
axes = FALSE, space = 0)
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)
dev.off()## png
## 2#======================================================================
soilTexture <- c("100% sand", "100% silt", "100% clay")
soilTexture_colors <- c("red", "green", "blue")
xy <- list(c(0.2, 1.05),
c(0.2, 0.75),
c(0.2, 0.45))
png("FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTexture))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTexture[i],
fill = soilTexture_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2#======================================================================
soilTextureTypes <- c("Sand", "Loamy sand", "Sandy loam",
"Loam", "Silt loam", "Silty clay",
"Clay", "Clay loam", "Sandy clay loam")
soilTextureTypes_colors <- c("#d73229", "#f16a15", "#9d6e48",
"#eded31", "#59b03c", "#54c4c4",
"#2d8dbe", "#345da9", "#a71b6a")
xy <- list(c(-0.08, 1), c(0.17, 1), c(0.55, 1),
c(-0.08, 0.7), c(0.17, 0.7), c(0.55, 0.7),
c(-0.08, 0.4), c(0.17, 0.4), c(0.55, 0.4))
png("FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTextureTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTextureTypes[i],
fill = soilTextureTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 1.8, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2#======================================================================
ecolComp <- c("100% grass", "100% brush", "100% wood")
ecolComp_colors <- c("red", "green", "blue")
xy <- list(c(0.2, 1.05),
c(0.2, 0.75),
c(0.2, 0.45))
png("FigTerrains_tempFiles/FigTerrains_legend_ecolComp.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(ecolComp))
{
legend(xy[[i]][1], xy[[i]][2],
legend = ecolComp[i],
fill = ecolComp_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2#======================================================================
coverTypes <- c("Wood-grass", "shrubland", "woodland")
coverTypes_colors <- c("#c1c127", "#366d24", "#437c47")
xy <- list(c(0.2, 1.05),
c(0.2, 0.75),
c(0.2, 0.45))
png("FigTerrains_tempFiles/FigTerrains_legend_coverType.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(coverTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = coverTypes[i],
fill = coverTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 1.8, pt.cex = figScale * 5)
}
rm(i)
dev.off()## png
## 2Create blank image:
png("FigTerrains_tempFiles/FigTerrains_blank.PNG",
width = 10,
height = 10)## Warning in png("FigTerrains_tempFiles/FigTerrains_blank.PNG", width = 10, :
## 'width=10, height=10' are unlikely values in pixelsdev.off()## png
## 2Compose final figure
Load input images:
elevTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,".PNG")
soilTextTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTexture.PNG")
soilTextTypeTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTextureType.PNG")
ecolCompTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_ecolComp.PNG")
coverTypeTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_coverType.PNG")
numberOfColumns = 10
numberOfRows = 11
# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("FigTerrains_tempFiles/FigTerrains_seed.PNG",
"FigTerrains_tempFiles/FigTerrains_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilTextType.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_ecolComp.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_coverType.PNG",
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[1], ".PNG"),
elevTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
ecolCompTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
coverTypeTerrains[1],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[2], ".PNG"),
elevTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
ecolCompTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
coverTypeTerrains[2],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[3], ".PNG"),
elevTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
ecolCompTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
coverTypeTerrains[3],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[4], ".PNG"),
elevTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
ecolCompTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
coverTypeTerrains[4],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[5], ".PNG"),
elevTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
ecolCompTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
coverTypeTerrains[5],
# ---
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_ecolComp.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_coverType.PNG"
)
grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)Create Figure:
widths = 85 * c(leftWidth/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 85 * c(topHeight/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
bottomHeight/figHeight) * figScale
lay <- matrix(1:(numberOfColumns * numberOfRows), ncol = numberOfColumns, byrow = TRUE)
png("figures/terrains.png", width = figWidth, height = figHeight)
gridExtra::grid.arrange(grobs = grobs,
layout_matrix = lay,
widths = unit(widths, "cm"),
heights = unit(heights, "cm"))
dev.off()## png
## 2knitr::include_graphics("figures/terrains.png")