Nothing
####
# Spring Phenology R functions v. 0.3
# AMW 8/16/2019
#
# This script contains functions describing numerous spring onset indices and false spring
# exposure indices based on several papers. The following specifically are included here:
# 1) Extended Spring Indices (SI-x) and supporting functions (Ault et al. 2015; Schwarz et al. 2006; Schwartz et al. 2013)
# 2) Early and Late False Spring (EFS and LFS) - Allstadt et al. 2015
# 3) False Spring Exposure Index (FSEI) - Peterson and Abatzoglou, 2014
# Early FSEI and Late FSEI are calculated using those indicators from the Allstadt et al. 2015 study.
#
# Notes by version
#
# v0.1
# - The chill hours and chill date functions are included based on SI-x, but it is my rough interpretation
# given a lack of documentation in literature. Use with caution.
# - The SI-x are created based on translation from the original Matlab code from Ault et al. (2015).
# However, there are a couple improvements made. First, the FLI and FBI calculations incorporate the bug
# fix by Allstadt et al. (2015). Second, the original Matlab code had a function which read in precalculated solar
# declination angle. In this R function, I've incorporated code for calculating solar declination from R's solrad package.
# - The FSEI formula in the functions is based on Peterson and Abatzoglou (2014), however, the definition
# of a false spring is based on Allstadt et al. (2015). The FSEI function can be used to determine false
# spring exposure with both EFS and LFS.
#
#
# v0.2 6/24/2019
# - incorporated the function soldec, which brings in the precalculated solar declination values from Ault et al (2015).
# As a result, the declination function is no longer called.
# - Added the function calc_si, which (like the original Matlab) serves as a wrapper and calculates all SI-x and FSEI values for one or multiple years.
# - There are situations where the last freeze during the year occurs in the winter following first leaf. This is another no damage situation, as plants would be dormant.
# As such, the damageindex function now contains an added conditional forcing the damage index to zero if the last freeze occurs in the winter following first leaf.
# - Tests with Livneh in the 3^5 domain demonstrated a situation where some locations will have years where tasmin is always >= the freeze value (frzval).
# In this situation, there is no damage or false springs, so the values for both are set to zero (i.e. no freeze damage, no false spring occurred),
# using an if statement at the end of the calc_si function.
#
#
# v0.3 8/16/2019
# - There is a chance that the a freeze will occur after the first bloom, but not occur between the first leaf and first bloom.
# In the prior version, the falsesprings function would have counted this occurrence as both a late false spring and an early false spring.
# Counting the event as both an early and late false spring would be logically inconsistent. A correction for this has been applied to the falsesprings function.
# Essentially the correction ensures that the falsesprings function will only mark that those freeze days occurring one week after first leaf, but before first bloom are considered.
# - In preparation for implementation in a formal R package, the unused declination function has been deleted.
# - No other changes from prior versions.
# - Cleaned up the daylength function to drop the unused yearlength argument
#
# Questions on code should be emailed to amwootte@ou.edu
#
#########################################
#########################################
# Calculate Spring Indices - Given appropriate inputs, this function will calculate all the spring indices
calc_si = function(TMAX,TMIN,lat){
# INPUTS
# TMAX, TMIN - Matrices of tasmax and tasmin values which are 366 rows x N columns (matching number of years).
# The matrix should always have 366 days, as the function can fill in gaps (such as missing leap days) on its own.
# Input units should also be degrees F, though it will warn and auto convert if TMAX and TMIN are provided in degrees K.
# lat - scalar latitude of location
#
# OUTPUTS - in a list format variable
# FLImat - Matrix of first leaf which is N rows (matching number of years) by 4 columns
# The 4 columns correspond to the mean first leaf and the first leaf for plants 1-3 [mean,plant 1, plant 2, plant 3]
# FBImat - Alike to the FLImat, but for first bloom.
# DMGmat - Alike to the FLImat, but for the damage index.
# lastfreeze - Vector containing the day of year value (1-366) for the day of last freeze in spring. Should be length of N years.
# FSmat - Matrix containing false spring indicators. Should be N rows (matching number of years) by 2 columns.
# The two columns represent each of the false spring indicators [early false spring, late false spring]
# A value of 1 in either column represents that a false spring occurred. A value of 0 represents that a false spring did not occur
# FSEImat - False Spring Exposure Index (FSEI) values. This should be a vector of length 2.
# There is one value of FSEI each for the early and late false springs. [early FSEI, late FSEI]
#
# Note: Thresholds for base temperature (baset) and freezing temperature (frzval) are fixed in this code. This matches Ault et al. (2015),
# but it could be changed to an extra argument for this function later.
# Note: FSmat and FSEImat are current calculated using the first values of FLImat and FBImat. This is also consistent with prior code, but
# could be altered to provide the FSmat and FSEImat based on all four values of FLImat and FBImat.
if(is.matrix(TMAX)==FALSE | is.matrix(TMIN)==FALSE){
stop("TMAX or TMIN are not matrices",.call=TRUE)
}
if(is.numeric(TMAX)==FALSE | is.numeric(TMIN)==FALSE | is.numeric(lat)==FALSE){
stop("TMAX, TMIN, or lat are not numeric",.call=TRUE)
}
if(all(TMAX[which(is.na(TMAX)==FALSE)]>200)==TRUE & all(TMIN[which(is.na(TMIN)==FALSE)]>200)==TRUE){
warning("TMAX and TMIN appear to be in degrees K, so I converted them to degrees F")
TMAX = (TMAX-273.15)*9/5+32
TMIN = (TMIN-273.15)*9/5+32
}
baset = 31 # baset is 31F, 0C, 273.15K
frzval = 28 # -2.2C, 28F - hard freeze temperature
# initialize matrices for outputs
FLImat = FBImat = DMGmat = matrix(NA,nrow=ncol(TMAX),ncol=4)
FSmat = matrix(NA,nrow=ncol(TMAX),ncol=2)
lastfreeze = c()
for(i in 1:ncol(TMAX)){
tasmax= as.vector(TMAX[,i]) # getting tasmin and tasmax into vectors
tasmin = as.vector(TMIN[,i])
# if values are missing from tasmax or tasmin, the average of the 30 days is used.
# this section fills in missing values, similarly to the code of Ault et al. 2015
for(x in seq(1,211,by=30)){
start = x
end = x+29
flagm = 0
tmaxNAidx = which(is.na(tasmax[start:end])==TRUE)
tminNAidx = which(is.na(tasmin[start:end])==TRUE)
cntmax = length(tasmax[start:end])-length(tmaxNAidx)
cntmin = length(tasmin[start:end])-length(tminNAidx)
tmax = tasmax[start:end]
tmin = tasmin[start:end]
idx = start:end
if(length(tmaxNAidx)>0){
avgmax = mean(tmax[-tmaxNAidx])
tasmax[idx[tmaxNAidx]]=avgmax
}
if(length(tminNAidx)>0){
avgmin = mean(tmin[-tminNAidx])
tasmin[idx[tminNAidx]]=avgmin
}
if(cntmax < 20 | cntmin <20) flagm = 1
if(flagm == 1){
warning(paste("There is a lot of missing data in TMAX and TMIN for column ",i,sep=""))
}
}
DOY = 1:length(tasmax) # the example I have here is with a full year of data, so you may which to calculate with less than one full year.
daystop = max(DOY) # for daylength, you can have daylength calculated until a particular DOY. I just have it run the whole year here.
dlen = daylength(daystop,lat=lat) # Calculating daylength
FLImat[i,2] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=1,type="leaf",plant=1) # Following with the Matlab code in Ault et al. 2015, the FLI numbers are ordered such that the mean FLI is the first number in the output vector
FLImat[i,3] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=1,type="leaf",plant=2) # and the three following numbers are the FLI for each plant used to calculate first leaf.
FLImat[i,4] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=1,type="leaf",plant=3)
FLImat[i,1] = round(mean(FLImat[i,2:4],na.rm=TRUE))
FBImat[i,2] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=FLImat[i,2],type="bloom",plant=1) # FBI has a similar vector out to FLI above. FBI also takes the FLI for each plant as it's reference date.
FBImat[i,3] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=FLImat[i,3],type="bloom",plant=2) # This is shown in the Matlab code of Ault et al. 2015
FBImat[i,4] = leafindex(tasmax=tasmax,tasmin=tasmin,daylen=dlen,baset=baset,refdate=FLImat[i,4],type="bloom",plant=3)
FBImat[i,1] = round(mean(FBImat[i,2:4],na.rm=TRUE))
freezeinfo = freezedates(tasmin=tasmin,frzval=frzval,DOY=DOY)
lastfreeze[i] = freezeinfo$lastfreeze
freezedata = subset(freezeinfo[[4]],DOY>=1 & DOY<=180) # table containing the tasmin, and day of year for the hard freeze.
if(is.na(lastfreeze[i])==FALSE){
DMGmat[i,2] = damageindex(leafidx=FLImat[i,2],lastfreeze=freezeinfo$lastfreeze) # per the Schwarz and Ault papers, the last freeze is all that is needed calculate damage.
DMGmat[i,3] = damageindex(leafidx=FLImat[i,3],lastfreeze=freezeinfo$lastfreeze) # also with a similar vector output to FLI and FBI.
DMGmat[i,4] = damageindex(leafidx=FLImat[i,4],lastfreeze=freezeinfo$lastfreeze)
DMGmat[i,1] = round(mean(DMGmat[i,2:4]))
fs = falsesprings(SI=c(FLImat[i,1],FBImat[i,1]),freezedata = freezedata) # takes the freezedata table to determine if there's an early false spring or a late false spring.
FSmat[i,1] = fs$EFS
FSmat[i,2] = fs$LFS
} else {
message("Freeze didn't happen this year, so no damage or false springs")
DMGmat[i,1:4] = 0
FSmat[i,1:2] = 0
}
message("Calcs finished for year: ",i," / ",ncol(TMAX))
}
FSEImat = apply(FSmat,2,FSEI)
list(FLImat=FLImat,FBImat=FBImat,DMGmat=DMGmat,lastfreeze=lastfreeze,FSmat=FSmat,FSEImat=FSEImat)
}
####
# Solar declination - from the original matlab package, precalculated
soldec = function(DOY){
# INPUTS
# DOY - scalar day of year, from 1 to 366 or 1 to 365 depending on the calendar of the date
#
# OUTPUTS
# spring climatological calendar day of year for the location and calendar day of year matching the appropriate solar declination
decvals=c(-23.16,-23.08,-23,-22.92,-22.83,-22.73,-22.62,-22.5,-22.38,-22.24,-22.11,-21.96,-21.81,-21.65,-21.48,-21.31,-21.13,-20.94,
-20.75,-20.55,-20.34,-20.13,-19.91,-19.68,-19.45,-19.21,-18.97,-18.72,-18.46,-18.2,-17.94,-17.68,-17.39,-17.11,-16.83,-16.53,-16.23,
-15.93,-15.62,-15.31,-15,-14.68,-14.36,-14.03,-13.7,-13.36,-13.03,-12.68,-12.34,-11.99,-11.64,-11.28,-10.93,-10.57,-10.2,-9.84,-9.47,
-9.1,-8.73,-8.35,-7.98,-7.6,-7.22,-6.83,-6.45,-6.06,-5.68,-5.29,-4.9,-4.51,-4.12,-3.72,-3.33,-2.93,-2.54,-2.15,-1.75,-1.36,-0.97,-0.57,
-0.17,0.22,0.62,1.01,1.41,1.8,2.19,2.58,2.98,3.37,3.76,4.14,4.53,4.92,5.3,5.68,6.06,6.44,6.82,7.19,7.57,7.94,8.31,8.67,9.04,9.4,9.76,
10.11,10.47,10.82,11.16,11.51,11.85,12.19,12.52,12.85,13.18,13.5,13.83,14.14,14.45,14.76,15.07,15.37,15.67,15.95,16.24,16.53,16.81,17.08,
17.35,17.61,17.87,18.13,18.38,18.62,18.87,19.09,19.32,19.54,19.76,19.97,20.18,20.38,20.57,20.76,20.94,21.12,21.29,21.46,21.61,21.77,21.91,
22.05,22.18,22.31,22.43,22.54,22.65,22.75,22.84,22.93,23.01,23.08,23.15,23.21,23.26,23.31,23.35,23.38,23.41,23.43,23.44,23.44,23.43,23.41,
23.38,23.35,23.31,23.26,23.21,23.15,23.08,23.01,22.93,22.84,22.75,22.65,22.54,22.43,22.31,22.18,22.05,21.91,21.77,21.61,21.46,21.29,21.12,
20.94,20.76,20.57,20.38,20.18,19.97,19.76,19.54,19.32,19.09,18.87,18.62,18.38,18.13,17.87,17.61,17.35,17.08,16.81,16.53,16.24,15.95,15.67,
15.37,15.07,14.76,14.45,14.14,13.83,13.5,13.18,12.85,12.52,12.19,11.85,11.51,11.16,10.82,10.47,10.11,9.76,9.4,9.04,8.67,8.31,7.94,7.57,7.19,
6.82,6.44,6.06,5.68,5.3,4.92,4.53,4.14,3.76,3.37,2.98,2.58,2.19,1.8,1.41,1.01,0.62,0.22,-0.17,-0.57,-0.97,-1.36,-1.75,-2.15,-2.54,-2.93,-3.33,
-3.72,-4.12,-4.51,-4.9,-5.29,-5.68,-6.06,-6.45,-6.83,-7.22,-7.6,-7.98,-8.35,-8.73,-9.1,-9.47,-9.84,-10.2,-10.57,-10.93,-11.28,-11.64,-11.99,-12.34,
-12.68,-13.03,-13.36,-13.7,-14.03,-14.36,-14.68,-15,-15.31,-15.62,-15.93,-16.23,-16.53,-16.83,-17.11,-17.39,-17.68,-17.94,-18.2,-18.46,-18.72,-18.97,
-19.21,-19.45,-19.68,-19.91,-20.13,-20.34,-20.55,-20.75,-20.94,-21.13,-21.31,-21.48,-21.65,-21.81,-21.96,-22.11,-22.24,-22.38,-22.5,-22.62,-22.73,
-22.83,-22.92,-23,-23.08,-23.16,-23.23,-23.29,-23.34,-23.38,-23.42,-23.45,-23.47,-23.48,-23.49,-23.5,-23.5,-23.49,-23.48,-23.47,-23.45,-23.42,
-23.38,-23.34,-23.29,-23.23)
dayvals = c(307:366,1:306)
dayvals[DOY]
}
###
# Daylength calculator - supporting function, but is called directly for leaf index calculations
daylength = function(daystop,lat){
# INPUTS
# daystop - scalar day of year at which to end calculations
# lat - latitude for location of interest in decimal degrees
#
# OUTPUTS
# DAYLEN - vector of length equal to input daystop describing the total hours of daylight per day.
DAYLEN = c()
for(i in 1:daystop){
if(lat < 40){ # high vs. low latitude calculations for daylength
DLL = 12.14+3.34*tan(lat*pi/180)*cos(0.0172*soldec(i)-1.95)
} else {
DLL = 12.25 + (1.6164+1.7643*(tan(lat*pi/180))^2)*cos(0.0172*soldec(i)-1.95)
}
if(DLL < 1) DLL=1 # for very high latitudes, these ifs allow for having at least 1 hour of day or night.
if(DLL> 23) DLL = 23
DAYLEN[i]=DLL
}
DAYLEN
}
###
# Chill Hours calculator - supporting function for SI-x chill indices, not called directly.
chillh = function(tmax,tmin,daylen,baset){
# INPUTS
# tmax - scalar daily high temperature for a given day
# tmin - scalar daily low temperature for a given day
# daylen - scalar daylength for a given day (calculated from daylength function)
# baset - base temperature for determining chill hours. Typically this is 7.2C, but can be plant specific.
#
# OUTPUTS
# CHOUR - total chill hours for a given day. That is, number of hours the temperature fell below baset.
DTR = tmax-tmin # diurnal temperature range
IDL = floor(daylen) # rounds off daylen to an integer value
###
# Estimates hourly temperature curves during daylight hours
Thr = c()
Thr[1] = tmin
Thr[2:(IDL+1)]=DTR*sin(pi/(daylen+4)*(1:IDL))+tmin
###
# Estimates hourly temperature curve during night hours and puts them together with daylight hours.
TS1 = DTR*sin(pi/(daylen+4)*daylen)+tmin
if(TS1 <=0) TS1 = 0.01
Thr[(IDL+2):24] = TS1-(TS1-tmin)/(log(24-daylen))*log((1:(23-IDL)))
###
# chill hours calculation
CHRS = baset-Thr
CHRSflag = ifelse(CHRS<0,0,1)
CHOUR = sum(CHRSflag)
#list(CHOUR,Thr) # you can uncomment this line and comment out the next to see the hourly temperature curve.
CHOUR
}
###
# Chill date calculator - provides output to calculate the Composite Chill Date for SI-x.
chilldate = function(tasmax,tasmin,DOY,daylen,baset,plant=1){
# INPUTS
# tasmax - vector of daily high temperatures for a given year
# tasmin - vector of daily low temperatures for a given year
# daylen - vector of daylengths for a given year (calculated from daylength function).
# DOY - vector of day of year values for a given year (1:365 or 1:366 depending on calendar)
# *** tasmax, tasmin, daylen, and DOY must be the same length.
# baset - base temperature for determining chill hours. Typically this is 7.2C, but can be plant specific.
# plant - signifies plant type 1, 2, or 3. plant = 1 for lilac, plant = 2 for arnold red, plant =3 for zabelli
#
# OUTPUTS
# chillDOY - scalar, estimated day of year the chill requirement for the specified plant is met.
#
# Many questions regarding how this is calculated, use at your own risk.
###
# adjusts the calendar to begin from July 1 and wrap to June 30.
chillframe = data.frame(tasmax,tasmin,daylen,DOY)
if(nrow(chillframe)==365){
chillframe$DOYadj = ifelse(DOY<=180,DOY+185,DOY-180)
} else {
chillframe$DOYadj = ifelse(DOY<=180,DOY+186,DOY-180)
}
chillframe = chillframe[order(chillframe$DOYadj),]
###
# Calculates chill hours for each day and accumulated chill hours
chillframe$CH = NA
chillframe$ACH = NA
for(i in 1:nrow(chillframe)){
chillframe$CH[i] = chillh(chillframe$tasmax[i],chillframe$tasmin[i],chillframe$daylen[i],baset=baset)
if(i==1) chillframe$ACH[i]=chillframe$CH[i]
if(i>1) chillframe$ACH[i]=chillframe$ACH[(i-1)]+chillframe$CH[i]
}
###
# identifies the first day the where the chill requirements are met.
if(plant==1) chillidx = which(chillframe$ACH >= 1375)
if(plant==2 | plant==3) chillidx = which(chillframe$ACH >= 1250)
chillDOY = chillframe$DOY[chillidx[1]] # chill DOY is returned in the normal calendar year DOY.
chillDOY
}
###
# Growing degree hours calculator - supporting function for SI-x FLI and FBI calcs, not called directly.
growdh = function(tmax,tmin,daylen,baset){
# INPUTS
# tmax - scalar daily high temperature for a given day
# tmin - scalar daily low temperature for a given day
# daylen - scalar daylength for a given day (calculated from daylength function)
# baset - base temperature for determining growing degree hours. Typically this is 0C, 31F, or 273.15K
#
# OUTPUTS
# GDHOUR - total growing degree hours for a given day. Similar calc to growing degree days.
if(tmin==0) tmin=0.01
if(tmax==tmin) tmax=tmax+0.01
DTR = tmax-tmin # diurnal temperature range
IDL = floor(daylen) # round daylength to nearest integer
###
# Estimates hourly temperature curves during daylight hours
Thr = c()
Thr[1] = tmin
Thr[2:(IDL+1)]=DTR*sin(pi/(daylen+4)*(1:IDL))+tmin
###
# Estimates hourly temperature curves during night hours and strings together with daylight hours.
TS1 = DTR*sin(pi/(daylen+4)*daylen)+tmin
if(TS1 <=0) TS1 = 0.01
Thr[(IDL+2):24] = TS1-(TS1-tmin)/(log(24-daylen))*log((1:(23-IDL)))
###
# growing degree hours calculation
GHRS = Thr-baset # like growing degree days, this uses the actual temperature difference
GHRS = ifelse(GHRS<0,0,GHRS) # values less than 0 are made zero
GDHOUR = sum(GHRS)
GDHOUR
}
###
# Synoptic Event Identifier - supporting function for SI-x FLI and FBI calcs, not called directly.
synval = function(GDH,lagGDH,type="leaf"){
# INPUTS
# GDH - scalar growing degree hours for a given day
# lagGDH - vector growing degree hours for the 7 previous days
# type - for FLI calcs, type="leaf"
# for FBI calcs, type="bloom"
#
# OUTPUTS - list with the following
# synflag - synoptic event flag. Has a value of 1 if a synoptic event happened, value of 0 otherwise.
# dde2 - sum of growing degree hours for current day and previous two days
# dd57 - sum of growing degree hours 5-7 prior to current day.
##
# Synoptic event definition depends on if one is interested in FLI or FBI
SYNFLAG=0
if(type == "leaf") LIMIT=637
if(type == "bloom") LIMIT=2001
if(type != "leaf" & type !="bloom") stop("Improper event type",.call=TRUE)
##
# Synoptic event occurs if the growing degrees of most recent three days are >= LIMIT
VALUE=GDH+lagGDH[1]+lagGDH[2]
if(VALUE >= LIMIT){
SYNFLAG=1
} else {
SNYFLAG=0
}
##
# gathering outputs
DDE2=GDH+lagGDH[1]+lagGDH[2]
DD57=sum(lagGDH[5:7])
output = list(synflag = SYNFLAG,ddE2=DDE2,dd57=DD57)
output
}
####
# FLI and FBI - primary workhorse for SI-x function, calls multiple supporting functions.
leafindex = function(tasmax,tasmin,daylen,baset,refdate,type="leaf",plant=1, verbose=FALSE){
# INPUTS
# tasmax - vector of daily high temperatures for a given year
# tasmin - vector of daily low temperatures for a given year
# daylen - vector of daylength for the given year
# *** tasmax, tasmin, and daylen must be the same length.
#
# baset - scalar, base temperature for growing degree hours calculation (typically 0C)
# refdate - scalar, reference day of year at which to begin calculations of FLI or FBI.
# for FLI, refdate = 1
# for FBI, refdate = FLI
# type - for FLI calcs, type="leaf"
# for FBI calcs, type="bloom"
# plant - scalar, signifies plant type 1, 2, or 3. plant = 1 for lilac, plant = 2 for arnold red, plant =3 for zabelli
# verbose - FALSE by default, set this to TRUE to help with debugging.
#
# OUTPUTS
# OUTDOY - The value of FLI or FBI (depending on the type input), for the specified plant
###
# legacy code intialization inherited from Matlab
ID2 = refdate
daymax = 240
Eflag = 0
AGDH = 0
SYNOP = 0
CNT = refdate-1
lagGDH = rep(0,7)
###
# coefficients for estimating FLI and FBI by plant
Alf=list(c(3.306,13.878,0.201,0.153),c(4.266,20.899,0.000,0.248),c(2.802,21.433,0.266,0.000))
Alb=list(c(-23.934,0.116),c(-24.825,0.127),c(-11.368,0.096))
###
# while loop goes through all days until the requirements for FLI or FBI are met.
parametersout = matrix(NA,nrow=366,ncol=6)
lagGDHvals = matrix(NA,nrow=366,ncol=7)
while(CNT<daymax & Eflag==0){
###
# growing degree hours calculation
CNT = CNT+1
#message("DOY ",CNT," - doing calcs for MDSUM1")
GDH = growdh(tasmax[CNT],tasmin[CNT],daylen[CNT],baset)
###
# initializing lagGDH with faux values for the first day
if(CNT==1 & type=="leaf"){
lagGDH[1] = GDH
lagGDH[2] = GDH
}
lagGDHvals[CNT,] = lagGDH
###
# Identify if synoptic event occurred
synlist = synval(GDH,lagGDH,type=type)
###
# start accumulating synoptic events and growing degree hours once past the reference date
if(CNT>=refdate){
AGDH = GDH+AGDH
if(synlist[[1]]==1) SYNOP=SYNOP+1
}
###
# calculates the indicator (MDSUM1) for if first leaf or first bloom is achieved
# functions inherited from Matlab code translated to R
if(CNT>=refdate){
MDS0 = CNT-refdate # number of days past the reference date
parametersout[CNT,1] = MDS0
parametersout[CNT,2] = SYNOP
parametersout[CNT,3] = synlist[[2]]
parametersout[CNT,4] = synlist[[3]]
parametersout[CNT,5] = AGDH
parametersout[CNT,6] = GDH
if(type=="leaf"){
if(plant==1) MDSUM1=(Alf[[1]][1]*MDS0)+(Alf[[1]][2]*SYNOP)+(Alf[[1]][3]*synlist[[2]])+(Alf[[1]][4]*synlist[[3]])
if(plant==2) MDSUM1=(Alf[[2]][1]*MDS0)+(Alf[[2]][2]*SYNOP)+(Alf[[2]][3]*synlist[[2]])+(Alf[[2]][4]*synlist[[3]])
if(plant==3) MDSUM1=(Alf[[3]][1]*MDS0)+(Alf[[3]][2]*SYNOP)+(Alf[[3]][3]*synlist[[2]])+(Alf[[3]][4]*synlist[[3]])
}
if(type=="bloom"){
if(plant==1) MDSUM1=(Alb[[1]][1]*MDS0)+(Alb[[1]][2]*AGDH)
if(plant==2) MDSUM1=(Alb[[2]][1]*MDS0)+(Alb[[2]][2]*AGDH)
if(plant==3) MDSUM1=(Alb[[3]][1]*MDS0)+(Alb[[3]][2]*AGDH)
}
} else {
MDSUM1 = 1
}
#message("MDSUM1 = ",MDSUM1)
###
# First leaf or first bloom has been achieved when MDSUM1 >= 1000
# the if statements below trigger the end of the while loop and determine OUTDOY.
if(MDSUM1>=999.5 & Eflag==0){
if(type=="leaf"){
if(plant==1 | plant==3) OUTDOY = CNT
if(plant==2) OUTDOY = CNT+1
}
if(type=="bloom"){
OUTDOY = CNT
}
Eflag=1
}
###
# holding on to up to 7 previous day so of growing degree hours
lagGDH[2:7]=lagGDH[1:6]
lagGDH[1] = GDH
}
if(CNT>=daymax){
OUTDOY=NA
}
if(verbose==TRUE){
list(OUTDOY=round(OUTDOY),parameters=parametersout,lagGDHvals = lagGDHvals) # this is FLI or FBI depending on the type argument in the function call.
} else {
round(OUTDOY)
}
}
####
# Freeze information calculator - supporting function, called directly and output used in other functions.
freezedates = function(tasmin,frzval,DOY){
# INPUTS
# tasmin - vector of daily low temperatures for a given year
# DOY - vector of day of year values for the given year
# *** tasmin and DOY must be the same length.
#
# frzval - scalar, freeze temperature of interest (hard freeze is tasmin < -2.2C)
#
# OUTPUTS - list including the following
# firstfreeze - scalar, first day of year a hard freeze occurs in the fall
# lastfreeze - scalar, last day of year a hard freeze occurs in the spring
# freeseperiod - scalar, range of days between firstfreeze and lastfreeze
# freezedata - data.frame, table with tasmin, DOY, and adjusted DOYs for all hard freeze days
###
# initialize data table, adjust calendar, re-order to July 1 - June 30 calendar
frzframe = data.frame(tasmin,DOY)
if(nrow(frzframe)==365){
frzframe$DOYadj = ifelse(DOY<=180,DOY+185,DOY-180)
} else {
frzframe$DOYadj = ifelse(DOY<=180,DOY+186,DOY-180)
}
frzframe = frzframe[order(frzframe$DOYadj),]
###
# identify which rows of the table are hard freeze days
idx = which(frzframe$tasmin<=frzval)
###
# determine firstfreeze, lastfreeze, and freeze period, and write outputs list
lastfreeze = frzframe[idx[length(idx)],2]
firstfreeze = frzframe[idx[1],2]
freezeperiod = frzframe[idx[length(idx)],3]- frzframe[idx[1],3]
if(length(lastfreeze)==0) lastfreeze=NA
freezelist = list(firstfreeze=firstfreeze,lastfreeze=lastfreeze,freezeperiod=freezeperiod,freezedata=frzframe[idx,])
freezelist
}
####
# Damage Index calculator - SI-x index, called directly
damageindex = function(leafidx,lastfreeze){
# INPUTS
# leafidx - scalar, the FLI value
# *** in the literature FLI is used, but FBI could be input here also.
# lastfreeze - scalar, the day of year for the last freeze
#
# OUTPUTS
# DMG - damage index, increasing values indicate greater damage potential
idx = lastfreeze - leafidx # damage index calculation
DMG = ifelse(idx<0 | idx>=180,0,idx) # setting values where last freeze happens before first leaf to zero, i.e. no damage potential
# Also, if the last freeze occurs in the winter following first leaf, this is set to zero. i.e. no damage potential that spring.
DMG
}
####
# False Spring Indicator - based on Allstadt et al. (2015)
falsesprings = function(SI=c(60,65),freezedata){
# INPUTS
# SI - vector with two arguments, the FLI and the FBI
# *** FLI should be first.
# freezedata - data.frame of freeze data from the freezedates function
#
# OUTPUTS - FSlist - list with the following
# EFS - Early False Spring Indicator
# EFS occurs when a hard freeze happens 7 or more days after FLI and before FBI
# EFS = 1 when an early false spring occurs, EFS = 0 otherwise
# LFS - Late False Spring Indicator
# LFS occurs when a hard freeze happens any day after FBI
# LFS = 1 when a late false spring occurs, LFS = 0 otherwise
###
# determines the number of days where EFS or LFS conditions are met
EFSidx = which(freezedata$DOY>=(SI[1]+7) & freezedata$DOY<=SI[2])
LFSidx = which(freezedata$DOY>SI[2])
###
# Checks length of the above and sets EFS or LFS to 1 if the length is greater than zero.
EFS = ifelse(length(EFSidx)>0,1,0)
LFS = ifelse(length(LFSidx)>0,1,0)
###
# gathering output
FSlist= list(EFS=EFS,LFS=LFS)
FSlist
}
####
# False Spring Exposure Index - from Peterson and Abatzoglou, 2014
FSEI = function(falsespringvector){
# INPUTS
# falsespringvector - vector with N years of the EFS or LFS from false springs
#
# OUTPUTS
# FSEI - scalar, False Spring Exposure Index value.
FSEI = (sum(falsespringvector,na.rm=TRUE)/length(falsespringvector))*100
FSEI
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.