如何在 R 中编写 Walsh & Lawler 降雨季节性指数的计算代码

Question

我正在努力循环需要计数器来计算降雨季节性的多个变量。由于计数器无法正常工作或顺序不正确，我最终收到了错误。我想计算每年许多不同 GCM 的指数，以便稍后评估降雨季节性趋势。

我的代码尝试使用数据库 DF 计算降雨季节性，其结构如下：

structure(list(Date = structure(c(3651, 3651, 3651, 3651, 3651, 

3651, 3651, 3652, 3652, 3652), 类 = "日期"), 站点 = 结构(c(1L, 1L、1L、1L、1L、1L、1L、1L、1L、1L)，级别 =“Table_Mountain_NP”，类别 =“因子”)， GCM = 结构(c(13L, 15L, 7L, 1L, 5L, 4L, 14L, 13L, 15L, 7L)，级别 = c("ACCESS-CM2"，"AWI-CM-1-1-MR"，"CMCC-ESM2"， “CNRM-CM6-1”、“CNRM-ESM2-1”、“EC-Earth3-CC”、“EC-Earth3-Veg-LR”、 “GFDL-ESM4”、“INM-CM4-8”、“INM-CM5-0”、“IPSL-CM6A-LR”、“KIOST-ESM”、 “MIROC-ES2L”、“MIROC6”、“MRI-ESM2-0”、“NESM3”、“NorESM2-MM”、 “BCC-CSM2-MR”、“CanESM5”、“CESM2”、“CMCC-CM2-SR5”、“FGOALS-g3” ), 类 = "系数"), MaxTemp = c(16.75, 22.05, 22.45, 20.45, 22.05, 22.95, 25.25, 19.15, 17.55, 20.15), 最低温度 = c(15.25, 16.95、20.45、18.15、17.25、19.35、18.15、15.05、16.55、18.15 ), 雨 = c(5.72918, 1.07309, 0.14049, 0, 0.03815, 0.04369, 0, 1.37117, 2.41574, 0.40211), AveTemp = c(16, 19.5, 21.45, 19.3, 19.65, 21.15, 21.7, 17.1, 17.05, 19.15), 年份 = c(1979, 1979, 1979, 1979, 1979, 1979, 1979, 1980, 1980, 1980), 月份 = c(12, 12, 12, 12, 12, 12, 12, 1, 1, 1), 季节 = 结构(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), 级别 = c("秋季", "春天", "夏天", "冬天"), class = "因素")), row.names = c(NA, 10L), 类 = "data.frame")

DF <- dat_base 
VarName <- "Rain"
  DF$Year <- as.numeric(format(DF$Date, '%Y'))
  DF$Month <- as.numeric(format(DF$Date, '%m'))
  VarIndex <- which(colnames(DF) == VarName)

# Calculate the Rainfall for each year
YearlyTotal <- aggregate(DF[, VarIndex], by = list(Year = as.numeric(format(DF$Date,'%Y')), GCM = DF$GCM), FUN = sum, na.rm = T)
 
# Calculate the total monthly rainfall for each year and month
MonthlyTotals <- aggregate(DF[, VarIndex], by = list(Year = as.numeric(format(DF$Date,'%Y')), 
                                                       Month =  as.numeric(format(DF$Date,'%m')), GCM = DF$GCM), FUN = sum, na.rm = TRUE)

# Merge the YearlyTotal with MonthlyTotals
Season_data <- merge(MonthlyTotals, YearlyTotal, by = c("Year","GCM"), suffixes = c("_Month", "_Year"))
  
  #Excel formula (see [https://www.researchgate.net/post/Can-anyone-help-me-in-calculating-seasonal-index-Walsh-Lawler-1981-using-excel-spreadsheet]) 
  #SI = 1/B15*(ABS(B3-B15/12)+ABS(B4-B15/12)+ABS(B5-B15/12)+ABS(B6-B15/12)+ABS(B7-B15/12)+ABS(B8-B15/12)+ABS(B9-B15/12)+ABS(B10-B15/12)+ABS(B11-B15/12)+ABS(B12-B15/12)+ABS(B13-B15/12)+ABS(B14-B15/12))
  
#Calculate the Seasonality Index
  
#SI <- 1 / YearlyTotal[,3] * rowSums(abs(MonthlyTotals[,4] - YearlyTotal[,3] / 12))

SI <- as.data.frame(NA)  

i <- unique(Season_data$Year)[1]

output <- length(unique(Season_data$GCM))*length(unique(Season_data$Year))

for (m in 1: output){
  
  for (i in i:Season_data$Year[length(Season_data$Year)]) {
 current_data <- Season_data[c(Season_data$Year == i),] 
 gcm_levels <- levels(na.omit(current_data$GCM))
  
    SI[m,1] <- i
    
    for (gcm_level in gcm_levels){
    
    subset_data <- current_data[current_data$GCM == gcm_level, ]
  
      if (length(subset_data$Month) == 12) {
      
      for (k in 1:12) { 
        
         SI[m,2] <- gcm_level  
         SI[m,3] <- 1 / subset_data$x_Year[1] * sum(abs(subset_data$x_Month[k] - subset_data$x_Year[1] / 12))
         
      }
        }
       
       else {SI[m,3] <- "NA"} 
    
  }
      }
}

我认为我需要另一个计数器来跟踪 SI 中的当前行，因为循环现在嵌套在年份和 GCM 级别上。每年应该有一个 SI 输出。

如果有一个现有的函数或公式可以进行此计算而无需推导，那就太好了！

Answer 1

感谢那些回复的人。我能够使用 Zoo 和 HydroTSM 包结合 for 循环和 if 语句解决问题，如下所示：

dat$Year <- as.numeric(format(dat$Date,'%Y'))
dat$Month <- as.numeric(format(dat$Date,'%m'))

SI <-data.frame(matrix(NA, ncol = 3)) #create database for the output
names(SI) <- c("Year","GCM","SI")
k <- 1

for (i in (dat$Year[1]+1):2099){ #loop through all the years in the main dataset

for (j in 1: length(levels(dat$GCM))){ #loop through the sub-groups within year

GCM <- levels(dat$GCM)[j]

if  (nrow(dat[c(dat$GCM == GCM & dat$Year == i),])== 0) { #if no data exists for a subgroup, then just record NA for the SI index

SI[k,1] <- i
SI[k,2] <- GCM
SI[k,3] <- NA

} else { #if data does exist for the group, then conver the current data parcel to zoo format and calculate the SI using the hydroSTM package. 

datx <- dat[c(dat$GCM == GCM & dat$Year == i),] 
datx <- zoo(datx$Rain, seq(from = datx$Date[1], to = datx$Date[length(datx$Date)], by = 1))

SI[k,1] <- i
SI[k,2] <- GCM
SI[k,3] <- si(datx) #Save the result in the SI database

} 

k <- k + 1 #loop to the next subgroup

} #loop to the next year

}