CODEWITHSUNDEEP
rfunction
Fabián
Fabián

Reputation: 21

R function for running sum

Since there is no R package, I have a question concerning the creation of my 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

Upvotes: 0

Views: 159

Answers (4)

Patrick
Patrick

Reputation: 32244

You can use the filter() function for this.

WSD <- function(P, PET, n){
  # This generates WSD with a backwards looking window of n observations
  # Result is a time-series object
  wat_bal <- filter(P - PET, filter = rep(1, n), method = "convolution", sides = 1)

  # Convert to vector, if so desired
  wat_bal <- as.vector(wat_bal)
  
  # Produce the desired result
  data.frame(wat_bal = wat_bal, condition = ifelse(wat_bal > 0, "wet", "dry"))
}

In addition to doing what you want, it can also apply weights (all 1's in the example above just to sum) to the filter depending on the distance from the value to be computed, something which is typical in water-balance research, so this function may be of wider interest to you.

Upvotes: 0

Fabi&#225;n
Fabi&#225;n

Reputation: 21

Could this be the solution:

WSD <- function(P, PET, n){
    
    lag_sum = 0
    wat_bal = P - PET
 
  for(i in 1:n){ 
    
  lag_sum = lag_sum + lag(wat_bal, i)
   
    }
   return(lag_sum)
}

?

Upvotes: 0

Fabi&#225;n
Fabi&#225;n

Reputation: 21

@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))
}

Upvotes: 2

Mikko Marttila
Mikko Marttila

Reputation: 11878

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

wsd <- cumsum(p - pet)

Upvotes: 1

Related Questions