R function - Index-CodePudding

since there is no R package, I have a question concerning the creation of an own function to calculate the SAPEI index - standardized antecedent precipitation evapotranspiration index ("A standardized index for assessing sub-monthly compound dry and hot conditions with application in China" by Li et al., 2021). As I understood, to calculate the SAPEI index, I need to calculate first the accumulated daily difference between precipitation and PET (such as 3-month scale) for each calendar day.

The equation is as follows:

(1) Antecedent water surplus or deficit (WSD)

What I did:

WSD <- function(P, PET, n){
 
  wat_bal <- P - PET
  
  for(i in (n 1):length(wat_bal)){ 
   
    condition = ifelse(wat_bal > 0, "wet", "dry")

    return(data.frame(wat_bal, condition))
  }
}

Unfortunately, I am not getting the expected result. I know that my function is not complete, but I also do not know how to proceed further. Especially the n (number of previous days) is a problem. Could anyone help me out?

Thank you very much in advance Fabian

CodePudding user response：

Based on your description of the algorithm, it sounds like you may be able to just do:

wsd <- cumsum(p - pet)

CodePudding user response：

@MikkoMarttila Unfortunately, the cumsum() function did not work out and just accumulated the values day by day. Maybe I will quote a part of the article: "[...] The daily difference between precipitation and potential evapotranspiration was the calculated to estimate the water balance. To reflect dry and wet conditions of a given day, the antecedent water surplus or deficit (WSD) was calculated through the following equation: (see above), where n is the number of previous days, PET represents the potential evapotranspiration, and P represents precipitation. The WSD values can be aggregated at different timescales, such as 3, 6, 9 months, and so on. [...]" (Li et al., 2021).

  WSD <- function(P, PET, n){
 
  wat_bal <- cumsum(P - PET)
  
  for(i in (n 1):length(wat_bal)){ 
   
    condition = ifelse(wat_bal > 0, "wet", "dry")

    }
   return(data.frame(wat_bal, condition))
}