Welcome to the second script of Module 4! Here, we will continue to work on the chill unit models. You will learn how to:

calculate daily chill units from ERA5-Land temperature rasters
calculate and visualize accumulated chill units in maps
reclassify rasters

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(stringr)
library(rgdal)

dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_4_Chill_Unit_Models/M04_Part02_Chill-Unit-Models_Steps_3_4/"

Here, we work with modelled temperature data from ERA5-Land. This data has the advantage that it contains hourly time steps, so we do not have to infer hourly temperatures from daily temperatures. Let us only work with the data from two years, i.e. 2007 and 2008, to keep it simple:

era2007 <- readRDS(paste0(dir, "data/era_cropped_2007.rds"))
era2008 <- readRDS(paste0(dir, "data/era_cropped_2008.rds"))
era2008

## class      : RasterBrick 
## dimensions : 26, 34, 884, 8784  (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1  (x, y)
## extent     : 43.35, 46.75, 38.75, 41.35  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : X2008.01.01.00.00.00, X2008.01.01.01.00.00, X2008.01.01.02.00.00, X2008.01.01.03.00.00, X2008.01.01.04.00.00, X2008.01.01.05.00.00, X2008.01.01.06.00.00, X2008.01.01.07.00.00, X2008.01.01.08.00.00, X2008.01.01.09.00.00, X2008.01.01.10.00.00, X2008.01.01.11.00.00, X2008.01.01.12.00.00, X2008.01.01.13.00.00, X2008.01.01.14.00.00, ... 
## min values :        -2.254799e+01,        -2.236040e+01,        -2.303648e+01,        -2.415963e+01,        -2.499773e+01,        -2.543749e+01,        -2.235431e+01,        -1.726602e+01,        -1.581152e+01,        -1.462015e+01,        -1.387098e+01,        -1.379911e+01,        -1.413167e+01,        -1.434972e+01,        -1.606977e+01, ... 
## max values :          1.051582367,          0.963874241,          0.783585316,          0.906620326,          0.705622537,          0.845711905,          1.750811037,          2.361113413,          3.235758335,          4.111621425,          4.527016854,          4.642742853,          4.412509023,          3.936205173,          3.189467935, ...

We quickly plot one layer of “era2007” and add the outline of Armenia to get a feeling what the data looks like. The spatial resolution is 0.1 degrees:

armenia <- readOGR(paste0(dir, "data/shapefiles/gadm41_ARM_0.shp"), verbose=F)
plot(era2007[[2]], main= "Temperature on Jan 01 2007 at 1 AM (°C)")
plot(armenia, add=T)

2. Calculate accumulated chill units

We know want to repeat with the ERA5-Land data what we did in script 1, i.e. calculate daily chill units from hourly chill units, and then calculate accumulated chill units over time. We need a time series of each pixel to do so, so we first have to convert the rasters to a table. We can do that with the following steps:

## create a lookup table where we log the row and column positions of each cell in "era2007"
lookup <- data.frame(cellID = c(1:ncell(era2007)), ROW = rep(c(1:nrow(era2007)), each=ncol((era2007))), COL = rep(c(1:ncol(era2007))))

## convert "era2007" to an array
era2007_array <- as.array(era2007)

## create an empty table with one row for each day and one column for each cell
era2007_hourly <- as.data.frame(matrix(nrow=nlayers(era2007), ncol=ncell(era2007)), row.names=names(era2007))

## for each cell, paste the respective time series vector from "era2007_array" into "era2007_hourly"
for (j in 1:ncell(era2007))
{ era2007_hourly[,j] <- era2007_array[lookup$ROW[j],lookup$COL[j],] }


## repeat for "era2008" and bind resulting dataframes for 2007 and 2008 data together:
era2008_array <- as.array(era2008)
era2008_hourly <- as.data.frame(matrix(nrow=nlayers(era2008), ncol=ncell(era2008)), row.names=names(era2008))
for (j in 1:ncell(era2008))
{ era2008_hourly[,j] <- era2008_array[lookup$ROW[j],lookup$COL[j],] }

era_hourly <- rbind(era2007_hourly, era2008_hourly)

Let’s have a quick look at “era_hourly”, where each column contains the time series of one ERA5-Land cell:

head(era_hourly[c(1:7)])

	V1	V2	V3	V4	V5	V6	V7
X2007.01.01.00.00.00	-13.17662	-13.50797	-13.74307	-14.80653	-16.84209	-19.60977	-19.42583
X2007.01.01.01.00.00	-13.70409	-14.02082	-14.30952	-15.43633	-17.49138	-20.19206	-20.07877
X2007.01.01.02.00.00	-14.04762	-14.30587	-14.61650	-15.80543	-17.95063	-20.76947	-20.78043
X2007.01.01.03.00.00	-14.15969	-14.26079	-14.54219	-15.77132	-18.04564	-21.14223	-21.41144
X2007.01.01.04.00.00	-14.50321	-14.44961	-14.65548	-15.84685	-18.19670	-21.63559	-21.97302
X2007.01.01.05.00.00	-14.66766	-14.44108	-14.53488	-15.64829	-18.12970	-22.27512	-22.44689

We can now easily convert hourly temperatures to hourly chill units according to the same rules that we applied before:

CU_hourly <- era_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 a column for year, month and day to CU_hourly and summarize hourly chill units to daily chill Units:

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_daily <- CU_hourly %>% dplyr::group_by(year, month, day) %>% summarize(across(everything(),list(sum))) %>% as.data.frame()

Change negative daily chill units to zero:

CU_daily[CU_daily < 0 ] <- 0

Calculate cumulative daily chill Units. Remember that the cumulation restarts on August 1! We integrate two lines of code at the end of the loop that print the computation progress for us:

CU_daily_cum <- CU_daily
CU_daily_cum[,c(4:887)] <- NA

for (j in 4: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 if ((CU_daily_cum$month[i] == "08") & (CU_daily_cum$day[i]== "01")) { CU_daily_cum[i,j] <- 0   ## restart on August 1
      } else { CU_daily_cum[i,j] <- CU_daily[i,j] + CU_daily_cum[i-1,j] }
  }
  if(j %in% seq(from=100, to=length(CU_daily_cum), by=100)) {
  print(paste(j, length(CU_daily_cum), sep=" / ")) }
}

## [1] "100 / 887"
## [1] "200 / 887"
## [1] "300 / 887"
## [1] "400 / 887"
## [1] "500 / 887"
## [1] "600 / 887"
## [1] "700 / 887"
## [1] "800 / 887"

3. Plot maximum accumulated chill units

The object “CU_daily_cum” now contains the daily accumulated chill units from 2007 to 2008, with one column for each ERA5-Land cell. Let’s select the crop cycle that starts on August 1st 2007 and ends on July 31st 2008:

CU_daily_cum_sel <- CU_daily_cum[c(213:578),]

We convert the last line of “CU_daily_cum_sel” back to a raster. This line is the accumulated amount of chill units at the end of the cropcycle 2007/2008 and represents the maximum amount of chilling that was reached in each ERA5-Land cell:

## convert last row of "CU_daily_cum_sel" to a vector
CU_max <- as.vector(unlist(CU_daily_cum_sel[366,c(4:887)]))

## convert "CU_max" into a matrix with 34 rows and 26 columns. We transpose this matrix and convert it to a raster.
CU_max_raster <- raster(t(matrix(nrow=34, ncol=26, CU_max)))

## our raster lacks a proper extent, resolution and coordinate system. We apply the ones from the original ERA5-Land data
CU_max_raster

## class      : RasterLayer 
## dimensions : 26, 34, 884  (nrow, ncol, ncell)
## resolution : 0.02941176, 0.03846154  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 1115.5, 2925.5  (min, max)

extent(CU_max_raster) <- extent(era2007)
res(CU_max_raster) <- res(era2007)
crs(CU_max_raster) <- crs(era2007)
CU_max_raster

## class      : RasterLayer 
## dimensions : 26, 34, 884  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 43.35, 46.75, 38.75, 41.35  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1115.5, 2925.5  (min, max)

Plot “CU_max_raster” - This map shows how many chill units were accumulated from 2007 to 2008:

plot(CU_max_raster, col=cubicl(26), main=paste0("Maximum Accumulated Chill Units in 2008"))
plot(armenia, add=T)

4. Plot suitability map

Finally, we reclassify this raster according to the information in “stats.csv” using the function reclassify(). The object “reclass” contains the lower and upper limits that translate the continuous values from “CU_max_raster” into five new discrete values (1, 2, 3, 4, 5):

stats <- read.csv(paste0(dir, "/data/stats.csv"))

reclass <- as.matrix(data.frame(from=c(0,stats[,1]), to=c(stats[,1],5000), becomes=c(1:5)))
reclass

##        from     to becomes
## [1,]    0.0 1228.5       1
## [2,] 1228.5 1630.0       2
## [3,] 1630.0 1941.5       3
## [4,] 1941.5 2340.5       4
## [5,] 2340.5 5000.0       5

raster_reclass <- reclassify(CU_max_raster, reclass)
plot(raster_reclass)

Let us plot our final raster “raster_reclass” with ggplot! To do so, we have to convert the data to a “SpatialPixelsDataFrame” and then supply it to the function geom_tile():

raster_reclass_conv <- as.data.frame(as(raster_reclass, "SpatialPixelsDataFrame"))
raster_reclass_conv$layer <- as.factor(raster_reclass_conv$layer)

ggplot() +  
  geom_tile(data=raster_reclass_conv, aes(x=x, y=y, fill=layer), alpha=1) + 
  ggtitle("Suitability for Apple in 2008") +
  geom_polygon(data=armenia, aes(x=long, y=lat, group=group), fill=NA, color="black", size=0.25) +
  scale_fill_discrete("The amount of CU is:", type = as.vector(cubicl(5)),drop=FALSE, labels = c("insufficient, beyond current         observations", "below average but sufficient", "optimal", "above average", "above average, beyond current observations")) +
  coord_equal() + 
  xlab("") + ylab("")

Congratulations, you made it to the end of the second script! You can now solve the exercises 05 to 07.

Module 4: Chill Unit Models
Script 02: Steps 3 and 4

Max Hofmann

1. Preparations

2. Calculate accumulated chill units

3. Plot maximum accumulated chill units

4. Plot suitability map

Module 4: Chill Unit Models Script 02: Steps 3 and 4

Max Hofmann

1. Preparations

2. Calculate accumulated chill units

3. Plot maximum accumulated chill units

4. Plot suitability map

Module 4: Chill Unit Models
Script 02: Steps 3 and 4