Module 4: Chill Unit Models
Exercises: Solutions
Max Hofmann
Welcome to the exercises of module 4! Please try to solve the following exercises, and if you do not remember how to write the code for certain things, have a look again at the previous scripts. Programming is about learning by doing! So, do not hesitate to look for help in the internet, too. Remember that there is almost always more than just one way to solve a problem in R!
Part A: Script 1 (weather station measurements)
Exercise 01 (10 points)
Calculate hourly temperatures for the agrometeorological station assigned to your group. Plot the first 10 days of hourly temperature data in a plot and save it as PNG.
library(chillR)
library(ggplot2)
library(dplyr)
library(pals)
library(stringr)
library(raster)
library(rgdal)
dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_4_Chill_Unit_Models/solutions/"
## load data
meteo <- readRDS(paste0(dir, "data/weather_data.rds"))
coords <- read.csv(paste0(dir, "data/coordinates.csv"))
## choose "Odzun" station:
odzun <- meteo[[which(names(meteo)=="Odzun")]]
odzun_lat <- coords$Lati_new[coords$Station_Name=="Odzun"]
## calculate hourly temperatures
odzun_hourly <- stack_hourly_temps(odzun, latitude=odzun_lat)$hourtemps
names(odzun_hourly)[11] <- "Temp_hourly"
## create new object with first 10 days (=240 hours)
odzun_hourly_10days <- odzun_hourly[c(1:240),]
odzun_hourly_10days$count <- c(1:240)
## plot
plot1 <- ggplot(odzun_hourly_10days, aes(x = count)) +
ggtitle(paste0("Hourly temperature at Odzun station from Jan 1 to 10, 1995")) +
geom_line( aes(y=Temp_hourly)) +
scale_y_continuous( name = "Temperature (°C)") +
theme(text = element_text(size=20),
axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line")) + ylab("Temperature (°C)") + xlab("hours")
plot1
png("plot1.png", height = 600, width = 800)
plot1
dev.off()
## png
## 2Exercise 02 (10 points)
Calculate accumulated Utah chill units for the entire record of your station, taking into account that a new crop cycle starts each year on August 1st. Plot the trajectories of accumulated chill units for all crop cycles as lines, and save the plot as PNG.
## calculate hourly chill units
odzun_hourly$CU_hourly <- NA
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly < 1.5 ] <- 0
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 1.5 & odzun_hourly$Temp_hourly < 2.5 ] <- 0.5
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 2.5 & odzun_hourly$Temp_hourly < 9.2 ] <- 1
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 9.2 & odzun_hourly$Temp_hourly < 12.5 ] <- 0.5
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 12.5 & odzun_hourly$Temp_hourly < 16.0 ] <- 0
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 16.0 & odzun_hourly$Temp_hourly < 18.0 ] <- -0.5
odzun_hourly$CU_hourly[odzun_hourly$Temp_hourly >= 18.0 ] <- -1
## summarize hourly chill units to daily chill units
odzun_daily <- odzun_hourly %>% dplyr::group_by(Year, Month, Day, Station, DOY, Tmax, Tmin) %>% summarize(CU = sum(CU_hourly)) %>% as.data.frame()
## `summarise()` has grouped output by 'Year', 'Month', 'Day', 'Station', 'DOY', 'Tmax'. You can override using the `.groups` argument.
odzun_daily$CU[odzun_daily$CU < 0] <- 0
## create new columns in "odzun_daily"
odzun_daily$CU_cum <- NA
odzun_daily$CropCycle <- NA
cropcycle_temp <- 1995
## for each day, calculate accumulated chill units and assign crop cycle ID
for (i in 1:NROW(odzun_daily))
{ ## exception 1: for the first entry in "odzun_daily", the cumulative CU value is the CU value from that day
if (i==1) { odzun_daily$CU_cum[i] <- odzun_daily$CU[i]
} ## exception 2: on August 1st of each year, cumulative CU are reset to 0
else if ((odzun_daily$Month[i] == 8) & (odzun_daily$Day[i]== 1)) {
odzun_daily$CU_cum[i] <- 0
odzun_daily$CropCycle[i] <- cropcycle_temp
cropcycle_temp <- cropcycle_temp + 1
}
## regular case: cumulative CU is calculated as the sum of all previous CU values
else {
odzun_daily$CU_cum[i] <- odzun_daily$CU_cum[i-1] + odzun_daily$CU[i]
odzun_daily$CropCycle[i] <- cropcycle_temp
}
}
## do some data editing
odzun_daily <- odzun_daily[which(odzun_daily$CropCycle %in% c(1996:2020)),]
odzun_daily$DOS <- 1
for (i in 2:NROW(odzun_daily))
{ odzun_daily$DOS[i] <- odzun_daily$DOS[i-1]+1
if ((odzun_daily$Month[i] == 8) & (odzun_daily$Day[i]== 1)) { odzun_daily$DOS[i] <- 1 }
}
odzun_daily$Tavg <- (odzun_daily$Tmin + odzun_daily$Tmax)/2
odzun_daily$CropCycle_fac <- as.factor(odzun_daily$CropCycle)
odzun_daily$date <- paste(str_pad(odzun_daily$Month, 2, pad="0"), str_pad(odzun_daily$Day, 2, pad="0"), sep="-")
## plot
plot2 <- ggplot(odzun_daily, aes(x = DOS, y = CU_cum, colour=CropCycle_fac)) +
geom_line(size=1) +
ggtitle("Accumulated Utah chill units at Odzun station") +
theme(text = element_text(size=20),
axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line"),
legend.title = element_blank()) +
xlab("date") +
guides(fill=guide_legend(ncol=10)) +
ylab("Accumulated chill units") +
scale_x_continuous(labels = odzun_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)))
plot2
png("plot2.png", height = 600, width = 800)
plot2
dev.off()
## png
## 2Exercise 03 (10 points)
For each year and the respective DOY when bud bursting was observed for the crop assigned to your group, determine the amount of accumulated chill units at your station. Plot these values in a boxplot and save it as PNG. From this data, determine the minimum, 1st quartile, 3rd quartile and maximum amount of chill units during bud bursting of your crop and save these four values as CSV.
## load phenological data and select apple
pheno <- readRDS(paste0(dir, "data/phenological_data.rds"))
apple_raw <- pheno[[1]][,c(1:4,6)]
apple_raw <- apple_raw[complete.cases(apple_raw),]
apple <- filter(apple_raw, station == "Odzun", year != 1995)
## find CU values for bud bursting dates in "odzun_daily", and write them to "apple"
apple$CU_cum <- NA
for (i in 1:NROW(apple))
{ apple$CU_cum[i] <- odzun_daily$CU_cum[which ( (odzun_daily$CropCycle == apple$year[i]) & (odzun_daily$DOY == apple$X02_Bud_bursting[i])) ]
}
## plot data
plot3 <- ggplot(data=apple, aes(x=crop, y=CU_cum, fill = crop, group = crop)) +
geom_jitter(aes(color=crop), size=2.5, alpha=0.4) +
geom_boxplot(alpha=0.7) +
theme(text = element_text(size=20),
axis.text.x=element_blank(),
legend.position="none") +
ylab("CU at Bud Bursting") + xlab("") + ggtitle("Accumulated Chill Units, Bud Bursting of Apple, Odzun station")
plot3
png("plot3.png", height = 600, width = 800)
plot3
dev.off()
## png
## 2
## calculate stats and export as CSV
stats3 <- as.matrix(summary(apple$CU_cum)[c(1,2,5,6)])
stats3
## [,1]
## Min. 1257.5
## 1st Qu. 1750.0
## 3rd Qu. 2152.0
## Max. 2454.5
write.csv(stats3, paste0(dir, "stats3.csv"))Exercise 04 (20 points)
Calculate the amount of accumulated chill units at bud bursting for the crop assigned to you, for all stations. From this data, determine the minimum, 1st quartile, 3rd quartile and maximum amount of chill units during bud bursting of your crop and save these four values as CSV.
## identify all stations for which your crop has data
stations <- sort(unique(apple_raw$station))
## select all stations for which your crop has data from "meteo"
stations_data <- meteo[which(names(meteo) %in% stations)]
## select latitude coordinates. Make sure they are in the same order as the stations are organized in "stations_data"!
names(stations_data)
## [1] "Ashtarak" "Dilijan" "Yeghvard" "Goris" "Martuni"
## [6] "Masrik" "Odzun" "Sisian" "Stepanavan"
stations_lat <- NA
for (i in 1:length(names(stations_data)))
{ stations_lat[i] <- coords$Lati_new[coords$Station_Name == names(stations_data)[i]] }
stations_data_hourly <- list()
for (i in 1:length(stations_data))
{ stations_data_hourly[[i]] <- stack_hourly_temps(stations_data[[i]], latitude=stations_lat[i])$hourtemps }
View(stations_data_hourly[[2]])
## from here, you can continue as above, but write a loop around everything and loop over each station, i.e. over each element in the list "stations_data". In the end, you would have to use "rbind()" to get all CU values at bud bursting into one final dataframe, and calculate the stats of this final dataframe.Part B: Script 2 (ERA5-Land data)
Exercise 05 (20 points)
Download ERA5-Land data from https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=form. Check “2m temperature”, the two years assigned to your group, all months, all days and all day times. Under “Geographical area”, choose a suitable sub-region defined by coordinates that entail all of Armenia. Download the file as NetCDF (it might take a while until your order is processed!) and import it to R Studio as a Raster Stack. Display the first layer in a map, add the shapefile of Armenia and save your map as PNG.
## See the PDF "How to download ERA5 data.pdf" for some more instructions.
## import ERA5 download as raster stack
era5 <- stack(paste0(dir, "data/era5data.nc"))
## Loading required namespace: ncdf4
## convert from Kelvin to Celsius
era5 <- era5 - 273.15
## import Armenia border
arm <- readOGR(paste0(dir, "shapefile/gadm41_ARM_0.shp"), verbose = F)
## plot first layer of era5 with Armenia border and export as PNG
png("plot5.png", height = 600, width = 800)
plot(era5[[1]], main = "Temperature on Jan 01 2005 at 00:00 AM")
plot(arm, add=T)
dev.off()
## png
## 2
plot(era5[[1]], main = "Temperature on Jan 01 2005 at 00:00 AM")
plot(arm, add=T)
Exercise 06 (10 points)
Calculate accumulated Utah chill units for the crop cycle contained in your two years of ERA5-Land data. Choose any cell and plot the trajectory of accumulated chill units in this cell as a line, and save the plot as PNG.
## discard all layers that are before Aug 1 2005 or after July 31 2006
era5_sel <- era5[[5089:13848]]
head(names(era5_sel))
## [1] "X2005.08.01.00.00.00" "X2005.08.01.01.00.00" "X2005.08.01.02.00.00"
## [4] "X2005.08.01.03.00.00" "X2005.08.01.04.00.00" "X2005.08.01.05.00.00"
tail(names(era5_sel))
## [1] "X2006.07.31.18.00.00" "X2006.07.31.19.00.00" "X2006.07.31.20.00.00"
## [4] "X2006.07.31.21.00.00" "X2006.07.31.22.00.00" "X2006.07.31.23.00.00"
## create a lookup table where we log the row and column positions of each cell in "era5_sel"
lookup <- data.frame(cellID = c(1:ncell(era5_sel)), ROW = rep(c(1:nrow(era5_sel)), each=ncol((era5_sel))), COL = rep(c(1:ncol(era5_sel))))
## convert "era5_sel" to an array
era5_sel_array <- as.array(era5_sel)
## create an empty table with one row for each day and one column for each cell
era5_sel_hourly <- as.data.frame(matrix(nrow=nlayers(era5_sel), ncol=ncell(era5_sel)), row.names=names(era5_sel))
## for each cell, paste the respective time series vector from "era5_sel_array" into "era5_sel_hourly"
for (j in 1:ncell(era5_sel))
{ era5_sel_hourly[,j] <- era5_sel_array[lookup$ROW[j],lookup$COL[j],] }
## convert hourly temperatures to hourly CU:
CU_hourly <- era5_sel_hourly
CU_hourly[CU_hourly < 1.5 ] <- 0
CU_hourly[CU_hourly >= 1.5 & CU_hourly < 2.5 ] <- 0.5
CU_hourly[CU_hourly >= 2.5 & CU_hourly < 9.2 ] <- 1
CU_hourly[CU_hourly >= 9.2 & CU_hourly < 12.5 ] <- 0.5
CU_hourly[CU_hourly >= 12.5 & CU_hourly < 16.0 ] <- 0
CU_hourly[CU_hourly >= 16.0 & CU_hourly < 18.0 ] <- -0.5
CU_hourly[CU_hourly >= 18.0 ] <- -1
## add some columns for year, month and day
CU_hourly$year <- substr(row.names(CU_hourly),2,5)
CU_hourly$month <- substr(row.names(CU_hourly),7,8)
CU_hourly$day <- substr(row.names(CU_hourly),10,11)
CU_hourly$date <- as.Date(substr(row.names(CU_hourly),2,11), format="%Y.%m.%d")
## summarize hourly CU to daily CU
CU_daily <- CU_hourly %>% dplyr::group_by(date, year, month, day) %>% summarize(across(everything(),list(sum))) %>% as.data.frame()
## `summarise()` has grouped output by 'date', 'year', 'month'. You can override using the `.groups` argument.
## Change negative daily chill units to zero
CU_daily[CU_daily < 0 ] <- 0
## calculate accumulated chill units
CU_daily_cum <- CU_daily
CU_daily_cum[,c(5:1685)] <- NA
for (j in 5:length(CU_daily_cum))
{ for (i in 1:NROW(CU_daily_cum))
{ if (i==1) { CU_daily_cum[i,j] <- CU_daily[i,j] ## take 1st value from CU_daily
} else { CU_daily_cum[i,j] <- CU_daily[i,j] + CU_daily_cum[i-1,j]
}
}
}
## plot data of cell no. 1000
plot6 <- ggplot(CU_daily_cum, aes(x = date, y = V1_1)) +
geom_line(size=1) +
ggtitle("Accumulated Utah chill units of cell no. 1000") +
theme(text = element_text(size=20),
axis.text.x=element_text(angle=45, hjust=1),
legend.key.size = unit(2,"line"),
legend.title = element_blank()) +
xlab("") + ylab("Accumulated chill units")
plot6 
png("plot6.png", height = 600, width = 800)
plot6
dev.off()
## png
## 2Exercise 07 (20 points)
Convert the accumulated chill units at the end of the crop cycle to a raster. Reclassify this raster according to the values determined in exercise 4 (alternatively, take the values from exercise 3), and plot the reclassified raster with ggplot(). Save your map as PNG.
CU_max <- as.vector(unlist(CU_daily_cum[365,c(5:1685)]))
## 41 rows and columns, because these are the dimensions of the era5 raster we downloaded!
## (depends on the bounding box that you define during the download request)
CU_max_raster <- raster(t(matrix(nrow=41, ncol=41, CU_max)))
extent(CU_max_raster) <- extent(era5)
res(CU_max_raster) <- res(era5)
crs(CU_max_raster) <- crs(era5)
## plot
plot(CU_max_raster, col=cubicl(26), main=paste0("Maximum Accumulated Chill Units in 2005/2006"))
plot(arm, add=T)
## define reclassification table
reclass <- as.matrix(data.frame(from=c(0,stats3[,1]), to=c(stats3[,1],5000), becomes=c(1:5)))
raster_reclass <- reclassify(CU_max_raster, reclass)
raster_reclass_conv <- as.data.frame(as(raster_reclass, "SpatialPixelsDataFrame"))
raster_reclass_conv$layer <- as.factor(raster_reclass_conv$layer)
plot(raster_reclass)
levels(raster_reclass_conv$layer) <- factor(c("2", "3", "4", "5", "1"))
plot7 <- ggplot() +
geom_tile(data=raster_reclass_conv, aes(x=x, y=y, fill=layer), alpha=1) +
ggtitle("Suitability for Apple in 2005/2006") +
geom_polygon(data=arm, aes(x=long, y=lat, group=group), fill=NA, color="black", size=0.5) +
scale_fill_discrete("The amount of CU is:", type = as.vector(cubicl(5)), drop = FALSE,
labels = c("insufficient",
"below average but sufficient",
"optimal",
"above average",
"above average, beyond current observations")) +
theme(text = element_text(size=20)) +
coord_equal() +
xlab("") + ylab("")
## Regions defined for each Polygons
plot7
png("plot7.png", height = 600, width = 800)
plot7
dev.off()
## png
## 2