Module 4: Chill Unit Models
Script 01: Steps 1 and 2
Max Hofmann
Welcome to the first script of Module 4! Here, we will start to work on the Chill unit models. You will learn how to:
- calculate hourly temperatures from daily temperature records
- calculate daily chill units
- calculate and visualize accumulated chill units
We follow the Utah chill unit approach.
Please read the following code instructions carefully, try to understand the code that follows each instruction, execute it and see what happens. Do not hesitate to change the code or write your own code, and execute it as well!
1. Preparations
As usual, we start by loading the packages we need, and by defining our working directory.
library(dplyr)
library(ggplot2)
library(pals)
library(chillR)
library(stringr)
dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_4_Chill_Unit_Models/M04_Part01_Chill-Unit-Models_Steps_1_2/"We load the data that we need for this script. While the coordinates are saved as a CSV file, the meteorological and phenological data are saved as R data files (.RDS). You can export many R object as such RDS files, which is very handy when the object is in a special format. For example, the meteorological and phenological data is a list, which you could not save as an ordinary CSV files. You can import RDS files using the function readRDS():
coords <- read.csv(paste0(dir, "data/coordinates.csv"))
meteo <- readRDS(paste0(dir, "data/weather_data.rds"))
pheno <- readRDS(paste0(dir, "data/phenological_data.rds"))The object “coords” contains the latitude and longitude coordinates of each agrometeorological station. We need the latitude later to calculate hourly temperatures from daily temperatures:
head(coords)| Station_Name | Lati_new | Longi_new | elevation |
|---|---|---|---|
| Bagratashen | 41.24528 | 44.82556 | 451 |
| Odzun | 41.06028 | 44.61028 | 1105 |
| Stepanavan | 41.00194 | 44.41278 | 1397 |
| Vanadzor | 40.83889 | 44.43444 | 1376 |
| Dilijan | 40.74111 | 44.86556 | 1256 |
| Ijevan | 40.87167 | 45.14722 | 732 |
The object “meteo” is a list, each list element is a table with the daily temperature record of one station. Some columns are in special formats that are required to calculate hourly temperatures later on:
names(meteo)## [1] "Ararat" "Areni" "Armavir" "Artashat" "Ashtarak"
## [6] "Bagratashen" "Dilijan" "Yeghvard" "Goris" "Ijevan"
## [11] "Martuni" "Masrik" "Meghri" "Odzun" "Sisian"
## [16] "Stepanavan" "Urtsadzor" "Vanadzor" "Yerevan_agro"head(meteo[[1]])| DATE | Year | Month | Day | Station | DOY | Tmax | Tmin | |
|---|---|---|---|---|---|---|---|---|
| 1 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Ararat | 1 | 3.4 | -0.8 |
| 12 | 1995-01-02 12:00:00 | 1995 | 1 | 2 | Ararat | 2 | 4.6 | -2.4 |
| 23 | 1995-01-03 12:00:00 | 1995 | 1 | 3 | Ararat | 3 | 2.3 | -4.9 |
| 26 | 1995-01-04 12:00:00 | 1995 | 1 | 4 | Ararat | 4 | -1.8 | -4.9 |
| 27 | 1995-01-05 12:00:00 | 1995 | 1 | 5 | Ararat | 5 | -1.8 | -5.3 |
| 28 | 1995-01-06 12:00:00 | 1995 | 1 | 6 | Ararat | 6 | 0.3 | -2.8 |
The object “pheno” is also a list, each list element is a table with the phenological observation record of one fruit type:
names(pheno)## [1] "apple" "apricot" "cornel" "peach" "pear" "quince"head(pheno[[1]])| year | station | crop | variety | X01_Bud_swelling | X02_Bud_bursting | X03_Emergence_1st_leaf | X04_Flower_bud_splitting | X05_Blossoming | X06_Cease_blossoming | X07_Fruit_formation | X08_Fruit_maturation | X09_Harvest | X10_Fall_leaf_colors | X11_Defoliation |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1995 | Odzun | apple | Bellefleur | 96 | 104 | 110 | 116 | 126 | 138 | 181 | 259 | 279 | 293 | 304 |
| 1996 | Odzun | apple | Bellefleur | 99 | 111 | 117 | 121 | 129 | 141 | 172 | 264 | NA | NA | 305 |
| 1997 | Odzun | apple | Bellefleur | 100 | 108 | 114 | 118 | 130 | 140 | 171 | 263 | 280 | 283 | 293 |
| 1998 | Odzun | apple | Bellefleur | 98 | 110 | 112 | 120 | 126 | 138 | 171 | 259 | 276 | 297 | 314 |
| 1999 | Odzun | apple | Bellefleur | 87 | 100 | 106 | 116 | 128 | 144 | 171 | 263 | 291 | 293 | 303 |
| 2000 | Odzun | apple | Bellefleur | 99 | 109 | 113 | 121 | 129 | 141 | 156 | 262 | NA | NA | 305 |
We will again work with the station sisian, and create a new object from “meteo”, accordingly. The object “sisian” contains daily measurements of minimum and maximum temperature for each day from January 01 1995 to September 30 2020. You can ignore the addendum “12:00:00 GMT”, this is just required to calculate hourly temperatures. Let us also extract the latitude of sisian station:
sisian <- meteo[[which(names(meteo)=="Sisian")]]
min(sisian$DATE)## [1] "1995-01-01 12:00:00 GMT"max(sisian$DATE)## [1] "2020-09-30 12:00:00 GMT"sisian_lat <- coords$Lati_new[coords$Station_Name=="Sisian"]2. Calculate accumulated chill units
First, let’s calculate hourly temperatures with the function stack_hourly_temps() from the package chillR. This function takes the latitude coordinate, because this determines daytime length and thereby daily air cooling and warming. We see that the new object “sisian_hourly” contains 24 entries for each day. We name the newly created column at the end accordingly:
sisian_hourly <- stack_hourly_temps(sisian, latitude=sisian_lat)$hourtemps
names(sisian_hourly)[11] <- "Temp_hourly"
head(sisian_hourly)| DATE | Year | Month | Day | Station | DOY | Tmax | Tmin | JDay | Hour | Temp_hourly | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 0 | -2.8 |
| 9406 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 1 | -2.8 |
| 18811 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 2 | -2.8 |
| 28216 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 3 | -2.8 |
| 37621 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 4 | -2.8 |
| 47026 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 5 | -2.8 |
To better understand the relationship between daily and hourly temperatures, let’s select the first 5 days (=120 rows) and plot three lines, one for each of the columns “Tmax”, “Tmin” and “Temp_hourly”:
sisian_hourly_10days <- sisian_hourly[c(1:120),]
sisian_hourly_10days$count <- c(1:120)
ggplot(sisian_hourly_10days, aes(x=count)) +
ggtitle(paste0("Temperature at Sisian station from Jan 1 to 5, 1995")) +
geom_line( aes(y=Tmax, colour="Daily Maximum Temperature", size="Daily Maximum Temperature")) +
geom_line( aes(y=Tmin, colour="Daily Minimum Temperature", size="Daily Minimum Temperature")) +
geom_line( aes(y=Temp_hourly, colour="Hourly Temperature", size="Hourly Temperature")) +
scale_y_continuous( name = "Temperature (°C)") +
theme(axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line")) +
xlab("date") + ylab("Temperature (°C)") +
scale_color_manual(name = "", values = c("Daily Maximum Temperature" = "red",
"Daily Minimum Temperature" = "blue",
"Hourly Temperature" = "black")) +
scale_size_manual(name = "", values = c("Daily Maximum Temperature" = 1,
"Daily Minimum Temperature" = 1,
"Hourly Temperature" = 1))
Calculate hourly chill units from hourly temperatures:
sisian_hourly$CU_hourly <- NA
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly < 1.5 ] <- 0
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 1.5 & sisian_hourly$Temp_hourly < 2.5 ] <- 0.5
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 2.5 & sisian_hourly$Temp_hourly < 9.2 ] <- 1
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 9.2 & sisian_hourly$Temp_hourly < 12.5 ] <- 0.5
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 12.5 & sisian_hourly$Temp_hourly < 16.0 ] <- 0
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 16.0 & sisian_hourly$Temp_hourly < 18.0 ] <- -0.5
sisian_hourly$CU_hourly[sisian_hourly$Temp_hourly >= 18.0 ] <- -1
head(sisian_hourly)| DATE | Year | Month | Day | Station | DOY | Tmax | Tmin | JDay | Hour | Temp_hourly | CU_hourly | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 0 | -2.8 | 0 |
| 9406 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 1 | -2.8 | 0 |
| 18811 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 2 | -2.8 | 0 |
| 28216 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 3 | -2.8 | 0 |
| 37621 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 4 | -2.8 | 0 |
| 47026 | 1995-01-01 12:00:00 | 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 1 | 5 | -2.8 | 0 |
Convert hourly chill units to daily chill units and change negative daily chill units to zero:
sisian_daily <- sisian_hourly %>% dplyr::group_by(Year, Month, Day, Station, DOY, Tmax, Tmin) %>% summarize(CU = sum(CU_hourly)) %>% as.data.frame()
sisian_daily$CU[sisian_daily$CU < 0] <- 0
head(sisian_daily)| Year | Month | Day | Station | DOY | Tmax | Tmin | CU |
|---|---|---|---|---|---|---|---|
| 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 2.5 |
| 1995 | 1 | 2 | Sisian | 2 | 6.5 | 1.6 | 19.5 |
| 1995 | 1 | 3 | Sisian | 3 | 9.6 | 1.8 | 15.0 |
| 1995 | 1 | 4 | Sisian | 4 | 9.3 | -4.6 | 9.0 |
| 1995 | 1 | 5 | Sisian | 5 | 10.3 | -3.1 | 12.0 |
| 1995 | 1 | 6 | Sisian | 6 | 4.7 | -0.7 | 8.0 |
Calculate cumulative daily chill units:
sisian_daily$CU_cum <- NA
for (i in 1:NROW(sisian_daily))
{ ## exception 1: for the first entry in "sisian_daily", the cumulative CU value is the CU value from that day
if (i==1) { sisian_daily$CU_cum[i] <- sisian_daily$CU[i]
} ## exception 2: on August 1st of each year, cumulative CU are reset to 0
else if ((sisian_daily$Month[i] == 8) & (sisian_daily$Day[i]== 1)) { sisian_daily$CU_cum[i] <- 0
}
## regular case: cumulative CU is calculated as the sum of all previous CU values
else { sisian_daily$CU_cum[i] <- sisian_daily$CU_cum[i-1] + sisian_daily$CU[i] }
}
head(sisian_daily)| Year | Month | Day | Station | DOY | Tmax | Tmin | CU | CU_cum |
|---|---|---|---|---|---|---|---|---|
| 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 2.5 | 2.5 |
| 1995 | 1 | 2 | Sisian | 2 | 6.5 | 1.6 | 19.5 | 22.0 |
| 1995 | 1 | 3 | Sisian | 3 | 9.6 | 1.8 | 15.0 | 37.0 |
| 1995 | 1 | 4 | Sisian | 4 | 9.3 | -4.6 | 9.0 | 46.0 |
| 1995 | 1 | 5 | Sisian | 5 | 10.3 | -3.1 | 12.0 | 58.0 |
| 1995 | 1 | 6 | Sisian | 6 | 4.7 | -0.7 | 8.0 | 66.0 |
As mentioned above, the cumulation restarts on August 1st. That means that one crop cycle spans the period from August 1 to July 31 of the next year. Let’s create a crop cycle ID for that. We start with the year 1995 and add 1 year whenever we pass August 1st.
sisian_daily$CropCycle <- NA
cropcycle <- 1995
for (i in 1:NROW(sisian_daily))
{ if ((sisian_daily$Month[i] == 8) & (sisian_daily$Day[i]== 1)) { cropcycle <- cropcycle+1 }
sisian_daily$CropCycle[i] <- cropcycle
}
head(sisian_daily)| Year | Month | Day | Station | DOY | Tmax | Tmin | CU | CU_cum | CropCycle |
|---|---|---|---|---|---|---|---|---|---|
| 1995 | 1 | 1 | Sisian | 1 | 2.1 | -2.8 | 2.5 | 2.5 | 1995 |
| 1995 | 1 | 2 | Sisian | 2 | 6.5 | 1.6 | 19.5 | 22.0 | 1995 |
| 1995 | 1 | 3 | Sisian | 3 | 9.6 | 1.8 | 15.0 | 37.0 | 1995 |
| 1995 | 1 | 4 | Sisian | 4 | 9.3 | -4.6 | 9.0 | 46.0 | 1995 |
| 1995 | 1 | 5 | Sisian | 5 | 10.3 | -3.1 | 12.0 | 58.0 | 1995 |
| 1995 | 1 | 6 | Sisian | 6 | 4.7 | -0.7 | 8.0 | 66.0 | 1995 |
Let’s edit our data a bit so we can create nice plots in the next step. We delete the crop cycles 1995 and 2021 because they are incomplete. We add the variable “DOS” (= day of season), calculate daily average temperature, convert the crop cycle ID to a factor variable and create a new date variable:
sisian_daily <- sisian_daily[which(sisian_daily$CropCycle %in% c(1996:2020)),]
sisian_daily$DOS <- 1
for (i in 2:NROW(sisian_daily))
{ sisian_daily$DOS[i] <- sisian_daily$DOS[i-1]+1
if ((sisian_daily$Month[i] == 8) & (sisian_daily$Day[i]== 1)) { sisian_daily$DOS[i] <- 1 }
}
sisian_daily$Tavg <- (sisian_daily$Tmin + sisian_daily$Tmax)/2
sisian_daily$CropCycle_fac <- as.factor(sisian_daily$CropCycle)
sisian_daily$date <- paste(str_pad(sisian_daily$Month, 2, pad="0"), str_pad(sisian_daily$Day, 2, pad="0"), sep="-")
head(sisian_daily)| Year | Month | Day | Station | DOY | Tmax | Tmin | CU | CU_cum | CropCycle | DOS | Tavg | CropCycle_fac | date | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 213 | 1995 | 8 | 1 | Sisian | 213 | 25.6 | 16.4 | 0 | 0 | 1996 | 1 | 21.00 | 1996 | 08-01 |
| 214 | 1995 | 8 | 2 | Sisian | 214 | 25.5 | 10.1 | 0 | 0 | 1996 | 2 | 17.80 | 1996 | 08-02 |
| 215 | 1995 | 8 | 3 | Sisian | 215 | 27.0 | 10.8 | 0 | 0 | 1996 | 3 | 18.90 | 1996 | 08-03 |
| 216 | 1995 | 8 | 4 | Sisian | 216 | 26.9 | 8.6 | 0 | 0 | 1996 | 4 | 17.75 | 1996 | 08-04 |
| 217 | 1995 | 8 | 5 | Sisian | 217 | 29.3 | 7.6 | 0 | 0 | 1996 | 5 | 18.45 | 1996 | 08-05 |
| 218 | 1995 | 8 | 6 | Sisian | 218 | 26.3 | 13.0 | 0 | 0 | 1996 | 6 | 19.65 | 1996 | 08-06 |
3. Plot chill unit data
Let’s plot the trajectories of cumulative chill units for every year. We see that the maximum amount of chilling reached in summer differs from year to year:
ggplot(sisian_daily, aes(x = DOS, y = CU_cum, colour=CropCycle_fac)) +
geom_line(size=1) +
ggtitle("Accumulated Utah chill units at Sisian station") +
theme(axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line"),
legend.title = element_blank(),
legend.position="bottom") +
xlab("date") +
guides(fill=guide_legend(ncol=10)) +
ylab("Accumulated chill units") +
scale_x_continuous(labels = sisian_daily$date[c(1,32,62,93,123,154,185,214,245,275,306,336,366)],
breaks=c(1,32,62,93,123,154,185,214,245,275,306,336,366)) +
scale_colour_manual(values = as.vector(kovesi.rainbow(25)))
We can also plot the accumulated chill units together with the daily chill units and daily mean temperature, e.g. for the year 2008:
sisian_daily_2008 <- filter(sisian_daily, CropCycle==2008)
ggplot(sisian_daily_2008, aes(x=DOS)) +
ggtitle("Temperature and Utah chill units at Sisian station") +
geom_line( aes(y=CU, colour="Daily Chill Units", size="Daily Chill Units")) +
geom_line( aes(y=Tavg, colour="Average Temperature", size="Average Temperature")) +
geom_line( aes(y=CU_cum / 50, colour="Accumulated Chill Units", size="Accumulated Chill Units")) +
scale_y_continuous( name = "Average temperature / Daily Chill Units",
sec.axis = sec_axis(~.*50, name="Accumulated Chill Units")) +
theme(axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line") ) +
xlab("date") + ylab("Accumulated Positive Chill Units") +
scale_x_continuous(labels = sisian_daily$date[c(1,32,62,93,123,154,185,214,245,275,306,336,366)],
breaks=c(1,32,62,93,123,154,185,214,245,275,306,336,366)) +
scale_color_manual(name = "", values = c("Daily Chill Units" = "red",
"Accumulated Chill Units" = "red",
"Average Temperature" = "black")) +
scale_size_manual(name = "", values = c("Daily Chill Units" = 0.5,
"Accumulated Chill Units" = 1,
"Average Temperature" = 0.5))
Let’s now have a look at our phenological observations and work with apple here. We are interested in the bud bursting observations, we we select them. We exclude all NAs and the year 1995 because we do not have a complete time series of CU accumulation for that year:
apple <- pheno[[1]][,c(1:4,6)]
apple <- apple[complete.cases(apple),]
apple <- filter(apple, station == "Sisian", year != 1995)
head(apple)| year | station | crop | variety | X02_Bud_bursting |
|---|---|---|---|---|
| 1996 | Sisian | apple | Bellefleur | 121 |
| 1997 | Sisian | apple | Bellefleur | 124 |
| 1998 | Sisian | apple | Bellefleur | 114 |
| 1999 | Sisian | apple | Bellefleur | 116 |
| 2000 | Sisian | apple | Bellefleur | 113 |
| 2001 | Sisian | apple | Bellefleur | 98 |
The table “apple” now contains for each year the DOY during which bud bursting was observed at Sisian station. Let us know find the respective amount of accumulated chill units for each of these dates:
apple$CU_cum <- NA
for (i in 1:NROW(apple))
{ apple$CU_cum[i] <- sisian_daily$CU_cum[which ( (sisian_daily$CropCycle == apple$year[i]) & (sisian_daily$DOY == apple$X02_Bud_bursting[i])) ]
}Let us plot the data distribution of in the column “CU_cum” in “apple” as a boxplot:
ggplot(data=apple, aes(x=crop, y=CU_cum, fill = crop, group = crop)) +
geom_jitter(aes(color=crop), size=1.5, alpha=0.4) +
geom_boxplot(alpha=0.7) +
theme(axis.text.x=element_text(angle=45, hjust=1), legend.position="none") +
ggtitle("CU at Bud Bursting") + xlab("") + ylab("Accumulated Chill Units at Bud Bursting")
And finally, we extract the minimum, 1st quartile, 3rd quartile and maximum values of accumulated chill units at bud bursting. In the next script, we will need these to determine optimal, sub-optimal and critical CU ranges (= suitability classes), so let’s save them as a CSV file:
stats <- summary(apple$CU_cum)[c(1,2,5,6)]
write.csv(as.vector(stats), paste0(dir, "data/stats.csv"), row.names=F)Congratulations, you made it to the end of the first script! You can now go on with script 02, or solve the exercises 01 to 04.