Predicting Peak Cherry Bloom: Part 1
Vaughn Shirey, Leslie Ries, and Gina Wimp
4/1/2022
Your name(s): Madeline Wang Members of your group: Kaitlynn Himelright, Caroline Castellino, Harnoor Sachar
MAKE SURE ALL NECESSARY PACKAGES ARE INSTALLED (IT SHOULD ASK ON LAUNCH - say YES!) If you don’t have it: >install.package(“cowplot”) or similar for other packages
Section 2.1: Custom Functions and Conditional Statements
Functions are tools that we can use in programming to help us implement custom algorithms in order to solve some of the computational challenges that we may face, especially with the use of larger datasets. You’ve already used some pre-written functions last week (and in weeks prior). One of the most commonly used functions in R is the vector function, written as “c()” where you fill in values between the parentheses to create a vector. Typically, most functions in R include parentheses and the values that you enter into them are called “arguments.”
Let’s explore some custom functions with these ideas in mind:
#--------------------------------------------------------------------------------------------------
# Suppose we want to write a custom function that will convert Fahrenheit to Celsius. We can -
# easily do something like this in R if we know the mathematics behind the conversion. -
#--------------------------------------------------------------------------------------------------
# First, we define the function just like we would create a new variable using the <- (arrow).
# Next, we need to tell the function that it can accept a number for the Fahrenheit values we
# want to convert. Finally, we need to include {} (open and closing brackets) behind function().
# The code inside these brackets will define how our function works and translates into english as CREATE A FUNCTION THAT TAKES TEMPERATURE IN F AND CALL IT tempConv
tempConv <- function(tempF){}
# Okay, so we have our first function, but right now it doesn't do anything since we haven't filled
# out the math between the brackets. Let's do that now. The conversion from Fahrenheit to Celsius
# is C=(F-32)*(5/9), try coding that into the function below, storing the new Celsius value into a variable called "tempC."
tempConv <- function(tempF){
tempC <- (tempF-32)*(5/9) # takes the value entered for "tempF" (when the tempConv function is implemented) and returns the value in the variable "tempC"
return(tempC) # This will return your calculated Celsius value when you call the function.
}
# Now let's call that function we just created to see if it works.
tempConv(tempF=32)## [1] 0QUESTION 2.1.1: Using your function from above, what is 41 degrees Fahrenheit in degrees Celsius?
tempConv(tempF=41)## [1] 5Awesome! We’ve made our own function using R. You can image how helpful this will be over typing out the formula to convert degrees Fahrenheit to degrees Celsius over and over again.
Functions usually include lots of different types of commands, including logic tests (true/false) or conditional tests (if variable meets some requirement, then do one thing, else do something else). Obviously, we can’t go into any depth here. But here are just a few illustrations.
Below are a few examples of logic statments that use the == statement syntax which means does this equal that. Predict with your partner whether each statment will be TRUE or FALSE?
756 == 756 # Statement 1 will be TRUE
8560 == 6590 # Statement 2 will be FALSE
1e3 == 1000 # Statement 3 will be TRUE
#Now run the logic tests and see how you did.
756 == 756 # Statement 1## [1] TRUE8560 == 6590 # Statement 2## [1] FALSE1e3 == 1000 # Statement 3## [1] TRUEBelow are a few examples of logic functions that use the is.__ syntax which allows you to test whether an object is a certain TYPE of object (which is often necessary to run certain procedures with it.
that. Predict with your partner whether each statment will be TRUE or FALSE?
is.numeric(100) # Statement 1 will be TRUE
is.character(“100”) # Statement 2 will be TRUE
This next logic statement tests three values: 1,2,NA so it will return three answers.
First, we need to concatonate the three values into one object c(1,2,NA) then we can ask if each value is an NA: is.na(c(1,2,NA)). Testing whether there is a value or if a cell is empty (in other words = NA) is one of the most common types of tests!
is.na(c(1,2,NA)) # Statement 3 will be FALSE, FALSE, TRUE
# Below are a few examples of logic functions that use the is.__ syntax which allows you to test whether an object is a certain TYPE of object (which is often necessary to run certain procedures with it. Predict with your partner
# whether each statement will be TRUE or FALSE.
is.numeric(100) # Statment 1## [1] TRUEis.character("100") # Statment 2## [1] TRUEis.na(c(1,2,NA)) # Statement 3## [1] FALSE FALSE TRUENow let’s look at a different conditional test and put it together with an if, else statement. Here, we are going to ask if a number that we enter into an object (called “value”) is even or odd.
# First, we'll start by creating the value that we want to test. You can set this value to any integer you want
value <- 5
# Here is how we would code a conditional statement that will tell us if the number stored in value is an even or odd number. We will divide by 2 and see if the remainder is equal to (==) 0. Note that if we just use one = sign, it will actually set the value equal to 0. To calculate the remainder we use %%. Then, we will use if/else statements to print out our determination.
if((value%%2)==0){
print("The value is even.")
} else {
print("The value is odd.")
}## [1] "The value is odd."# Go ahead and run this code block to see how conditional statements work. You can change the value to see if the code works properly.Alright, let’s make sure that if someone enters a string or character into our tempConv function that R won’t have a melt down. I’ve gone ahead and done this for you already:
# Here, I've recreated your tempConv function to "catch" character values and ask the user to
# input a valid numeric value instead.
tempConv <- function(tempF){
if(is.character(tempF)){
print("Please enter a numeric value.")
} else{
tempC <- (tempF-32)*(5/9)
return(tempC) # This will return your calculated Celsius value when you call the function.
}
}QUESTION 2.1.2: What will the function tempConv return when the input value is “thirty-two” vs. 32? When the value is “thirty-two,” the function tempConv will return “please enter a numeric value”. When the input value is 32, the function will return the temperature in degrees Celsuism in this case 0. Why did you get this output? We got this output because the first value was not a numeric value and the second one was. The code only reutrns the desired output when a numeric value is inputted.
# Test the value "thirty-two" (including quotes) below:
tempConv(tempF="thirty-two")## [1] "Please enter a numeric value."# Test the value 32 below:
tempConv(tempF=32)## [1] 0Awesome, you’re well on your way to becoming a master of functions and logic/conditional statements in R. This will give a little context as we step you through the different steps we take to test the predictions of the cherry blossom peak bloom model.
Section 2.2: Data Cleaning using Functions
One of the largest tourist attractions in Washington, DC is the bloom of the Yoshino cherry trees (Prunus x yedoensis) each year along the tidal basin. In this longstanding tradition of spring, predicting “peak bloom” has been a hobby for cherry blossom enthusiasts around the world. Today, we will use temperature measurements from DC’s national airport to predict when we might see peak bloom of Yoshino cherries around the tidal basin.
You should now see a new data table called “climDat” pop up in your environment panel to the upper-right. Clicking on this will allow you to explore the dataset. At first glance it looks like there are a lot of NA values…
Let’s do some quick data checks before we start our modeling. Note that when we download climate data from NOAA, it is returned in farenheit, but we need these to be in celcius. Since only 2023-2025 has been added since we last did this exercise, none of those rows have been converted. We can see this by checking for missing temperature values in one of the Celcius fields.
# We'll use some conditional and logical statments to handle this task! R let's us easily count how many TRUE and FALSE values we have by using the sum() function.
# Count how many temperature values in degrees Celsius are empty
sum(is.na(climDat$TmC))## [1] 730YIKES! Okay, we have some empty values here. What will we do? Fortunately, we’ve already written a function to convert Fahrenheit into Celsius in the first part of this lab. Let’s deploy that function now!
# We will use the tidyverse pipe syntax to do this easily. Tidyverse is a package that makes it easier to do data manipulations. Here, it is going to look for the empty (NA) celcius fields (is.na) then use the "tempConv" function to convert F to C. Note, we use a function here called "ifelse" to easily test a conditional statement and then deploy an algorithm depending on if the conditional returns TRUE or FALSE:
climDat <- climDat %>%
dplyr::mutate(TxC=ifelse(is.na(TxC), tempConv(TxF), TxC),
TnC=ifelse(is.na(TnC), tempConv(TnF), TnC),
TmC=ifelse(is.na(TmC), tempConv(TmF), TmC))
# Count how many temperature values in degrees Celsius are empty (this should be zero now)
sum(is.na(climDat$TmC))## [1] 0Great, now our data has been back filled for missing degree Celsius values.In the real world, data are often a lot more messy that what we just dealt with. It will ALWAYS be good data science practice to examine the dataset in full before deploying any sort of analysis.
Section 2.3: Assessing Model Performance
Let’s move onto coding in our actual model. Below is a large block of code that encodes our chill tests for Yoshino cherries. Remember, these values were derived from actual research on Yoshino and Kwanzan cherries! Don’t worry about the code which will be quite a bit more complicated than what we did above. But now you still will be able to recognize some of the elements.
First, we’ll define some constant values to use throughout the code below. The constant values come from the Chung et al. paper that are shown in Table 1.
#--------------------------------------------------------------------------------------------------
Rc <- -78.9 # This is the accumulated number of daily chill units (Cd) needed to release from dormancy.
Rh <- 221.2 # This is the accumulated number of daily heat units to get peak flowering.
Tc <- 4.3 # This is our threshold temperature.
#--------------------------------------------------------------------------------------------------
# A custom function, chill.test(), is used to compute all of our chill test results using ifelse logic. We went over their complicated system for calculating chill degrees. You will recognize the different tests to decide on how to calculate chill days from the powerpoint slide in today's presentation. -
#--------------------------------------------------------------------------------------------------
chill.test <- function(climateData){
dat <- climateData # Copy over our input climate data into a new variable
# Calculate our different chill test values.
dat <- dat %>% dplyr::mutate(CMN1=0,
CMN2=((TmC-TnC)-((TxC-Tc)/2)),
CMN3=-(TmC-TnC),
CMN4=-((TxC/(TxC-TnC))*(TxC/2)),
CMN5=-((TxC/(TxC-TnC))*(TxC/2)*((TxC-Tc)/2))
)
# Using if/else statements, find the true Cd value for a given case
dat <- dat %>% dplyr::mutate(Cd = ifelse(TnC >= 4.3, 0,
ifelse(TnC >= 0,
ifelse(TxC >= 4.3, CMN2, CMN3),
ifelse(TxC >= 4.3, CMN5, CMN4))))
return(dat)
}
#--------------------------------------------------------------------------------------------------
# A custom function, calc.Dc(), to sum up all of our Cd values to calculate Dc. -
#--------------------------------------------------------------------------------------------------
calc.Dc <- function(chillTestData, yearOfInterest=2012){
# copy over our input data into a new variable, filter for the flowering year of interest
dat <- dplyr::filter(chillTestData, BloomYear==yearOfInterest)
# iterate through each column, adding a new column Dc with the cumulative sum of Cd
dat <- dplyr::mutate(dat, Dc=0) # initialize the new column Dc
for(i in 2:nrow(dat)){
dat[i, ]$Dc <- dat[i-1, ]$Dc + dat[i, ]$Cd
}
return(dat)
}
chill.end<-function(x){
dat <- x %>% dplyr::mutate(DcRc = Dc/Rc)
#if(dat$DcRc[yday(dat$Date)==15]>=1){
# (start<-as.Date(paste0(dat$BloomYear[1],"-01-15"), format = "%Y-%m-%d"))
#} else {
if(dat$DcRc[yday(dat$Date)==60]>=1){
(start<-dplyr::filter(dat, DcRc>=1)[1,]$Date)
} else {
(start<-as.Date(paste0(dat$BloomYear[1],"-03-01"), format = "%Y-%m-%d"))
}
return(start)}Okay, we’ve coded in our custom functions for the chill tests. Now, for a given year (2021), let’s use these functions to calculate when the Yoshino cherries will have broken dormancy.
# Call our custom function and store the result in DC_chill
DC_chill <- chill.test(climDat)
# Call our other custom function and store the result in DC_chill.2021
DC_chill.2021 <- calc.Dc(DC_chill, yearOfInterest=2021)
# Finally, calculate Dc/Rc for our year of interest and extract the date the cherries broke
# through dormancy.
DC_chill.2021 <- DC_chill.2021 %>% dplyr::mutate(DcRc = Dc/Rc)
if(DC_chill.2021$DcRc[yday(DC_chill.2021$Date)==60]>=1){
(start<-dplyr::filter(DC_chill.2021, DcRc>=1)[1,]$Date)
} else {
print("Did not reach chill day threshold. Start set to ")
(start<-as.Date(paste0(DC_chill.2021$BloomYear[1],"-03-01"), format = "%Y-%m-%d"))
}## [1] "2021-01-14"Okay great, chills tests are done! Now let’s move onto the heat tests, where we can begin to predict the actual date of peak cherry blossom bloom. We’ll use a longer code block for this as well.
#--------------------------------------------------------------------------------------------------
# A custom function, heat.test(), is used to compute all of our heat test results by picking out the right equation using ifelse logic. -
#--------------------------------------------------------------------------------------------------
heat.test <- function(climateData){
dat <- climateData # copy over our input climate data into a new variable
# Calculate our different heat test values.
dat <- dat %>% mutate(HMN1=TmC-Tc,
HMN2=(TxC-Tc)/2,
HMN3=0,
HMN4=0,
HMN5=(TxC-Tc)/2
)
# Using if/else statements, find the true Ca value for a given case
dat <- dat %>% mutate(Ca = ifelse(TnC >= 4.3, HMN1,
ifelse(TnC >= 0,
ifelse(TxC >= 4.3, HMN2, HMN3),
ifelse(TxC >= 4.3, HMN5, HMN4))))
return(dat)
}
#--------------------------------------------------------------------------------------------------
# A custom function, calc.Dc(), to sum up all of our Cd values to calculate Dc. Note that this - function also requires the chill data! -
#--------------------------------------------------------------------------------------------------
calc.Dh <- function(heatTestData, yearOfInterest=2012, start){
# Copy over our input data into a new variable, filter for the flowering year of interest
dat <- filter(heatTestData, BloomYear==yearOfInterest & Date >= start)
# Iterate through each column, adding a new column Dc with the cumulative sum of Cd
dat <- mutate(dat, Dh=0) # Initialize the new column Dc
dat$Dh[1]<-dat$Ca[1]
for(i in 2:nrow(dat)){
dat[i, ]$Dh <- dat[i-1, ]$Dh + dat[i, ]$Ca
}
return(dat)
}Alright, let’s use these new functions to calculate the predicted peak bloom date in 2021
# Call our custom function and store the result in DC_heat
DC_heat <- heat.test(climDat)
# Call our other custom function and store the result in DC_heat.2021, remember to include the
# chill test results from earlier!
start.2021 <- chill.end(DC_chill.2021)
DC_heat.2021 <- calc.Dh(DC_heat, yearOfInterest=2021, start.2021)
# Finally, calculate Dh/Rh for our year of interest and extract the date the cherries broke
# through dormancy.
DC_heat.2021 <- DC_heat.2021 %>% dplyr::mutate(DhRh = Dh/Rh)
filter(DC_heat.2021, DhRh>=1)[1,]$Date## [1] "2021-03-25"You can find out the real peak bloom date on this website: https://cherryblossomwatch.com/2021-cherry-blossoms/
QUESTION 2.3.2: On what date did the Yoshino cherry trees projected to reach peak bloom in 2021? How close were we to the right date? The projected date was 25 March 2021. The real date was 28 March 2021, so our prediction was 3 days early.
Section 2.4: Assessing Model Performance
Okay, now we’ve made projections about breaking dormancy and peak bloom for a single year. Let’s see how accurate our model of cherry blossom phenology is when compared to actual data.
First, let’s run our chill and heat test functions several times to grab the predicted dates of peak bloom for 2011 through 2025.
## [1] "Processing for bloom year: 2012"
## [1] "Processing for bloom year: 2013"
## [1] "Processing for bloom year: 2014"
## [1] "Processing for bloom year: 2015"
## [1] "Processing for bloom year: 2016"
## [1] "Processing for bloom year: 2017"
## [1] "Processing for bloom year: 2018"
## [1] "Processing for bloom year: 2019"
## [1] "Processing for bloom year: 2020"
## [1] "Processing for bloom year: 2021"
## [1] "Processing for bloom year: 2022"
## [1] "Processing for bloom year: 2023"
## [1] "Processing for bloom year: 2024"
## [1] "Processing for bloom year: 2025"## [1] "2012-03-23" "2013-04-12" "2014-03-29" "2015-04-06" "2016-03-29"
## [6] "2017-04-05" "2018-03-11" "2019-03-29" "2020-02-23" "2021-03-25"
## [11] "2022-03-18" "2023-02-19" "2024-03-25" "2025-03-18"QUESTION 2.4.1: What are the predicted peak bloom dates for all years? 23 March 2012, 12 April 2013, 29 March 2014, 6 April 2015, 29 March 2016, 5 April 2017, 11 March 2018, 29 March 2019, 23 February 2020, 25 March 2021, 18 March 2022, 19 February 2023, 25 March 2024, 18 March 2025
Now let’s grab the actual peak bloom values and make a plot comparing predicted to actual dates.
Now let’s plot the actual peak bloom compared to our predicted peak bloom dates.
QUESTION 2.4.2: How do the model predictions do for 2012-2025 compared to the original data set in the Chung paper (Fig. 2B)? Describe how they differ The model predictions were quite different than the original data set in the Chung paper, and a lot of the model predictions had peak bloom occurring a lot earlier than the data set in the Chung paper. There was a lot more variability in the model predictions than the Chung paper data set.
QUESTION 2.4.3: Discuss with your group some possible explanations for why the performance is different? The predicted peak bloom dates show more variability because the model relies on simplified assumptions and temperature inputs that can fluctuate widely from year to year. In contrast, actual bloom dates are constrained by biological and environmental factors that limit how much the timing can vary in reality.
You are done with this exercise! To turn it in - make sure to put everyone’s name that you are turning this in for. Then, knit the .rmd file. This will create an .html file in your project folder (files tab to the right). Now, check off that .html file, chose the “more” gear icon and click “export”. Export to a local drive location, then email to everyone who is turning in this shared report. Finally, you should all upload to canvas separately.