【问题标题】:What is the most elegant way to calculate seasonal means with R?用 R 计算季节性均值的最优雅方法是什么?
【发布时间】:2013-09-29 11:18:33
【问题描述】:

我有平均空间时间序列和每日平均观测数据。

如何以最简单的方式计算季节性平均值?季节应遵循 DJF(=冬季:12 月、1 月、2 月)、MAM、JJA 和 SON 的气象命名法。

这意味着 12 月的值来自 x-1 年。

这里很好地展示了每月平均值的计算: How to calculate a monthly mean?

在计算季节性平均值时可以遵循这个想法。但是,有几个注意事项使它不是很透明,必须小心!

我还在以前的帖子中处理过这个问题的一小部分:How to switch rows in R?

现在是完整的故事:

0:制作随机时间序列

ts.pdsi <- data.frame(date = seq(
                from=as.Date("1901-01-01"), 
                to=as.Date("2009-12-31"), 
                by="day"))
ts.pdsi$scPDSI <- rnorm(dim(ts.foo)[1],  mean=1, sd=1)    # add some data

第一个:使用 seas 包并将季节添加到您的时间序列中,必须将其格式化为 data.frame。

library(seas)
# add moth/seasons
ts.pdsi$month  <- mkseas(ts.pdsi,"mon")   # add months
ts.pdsi$seas <- mkseas(ts.pdsi,"DJF")     # add seasons
ts.pdsi$seasyear <- paste(format(ts.pdsi[,1],"%Y"), 
                          ts.pdsi$seas ,sep="")   # add seasyears, e.g. 1950DJF

这给了

> head(ts.pdsi)
    date      scPDSI month seas seasyear
1 1901-01-01 -0.10881074   Jan  DJF  1901DJF
2 1901-02-01 -0.22287750   Feb  DJF  1901DJF
3 1901-03-01 -0.12233192   Mär  MAM  1901MAM
4 1901-04-01 -0.04440915   Apr  MAM  1901MAM
5 1901-05-01 -0.36334082   Mai  MAM  1901MAM
6 1901-06-01 -0.52079030   Jun  JJA  1901JJA

第二次:然后您可以按照上述方法使用 $seasyear 列计算季节性平均值

> MEAN <- tapply(pdsi$scPDSI, ts.pdsi$seasyear, mean, na.rm = T)
> head(MEAN)
1901DJF     1901JJA     1901MAM     1901SON     1902DJF     1902JJA 
-0.45451556 -0.72922229 -0.17669396 -1.12095590 -0.86523850 -0.04031273 

注意:春季 (MAM) 和夏季 (JJA) 因严格的字母排序而切换。

第三次:切换回来

foo <- MEAN
for(i in 1:length(MEAN)) {
    if (mod (i,4) == 2) {
        foo[i+1] <- foo[i]    #switch 2nd 3rd row (JJA <-> MAM)
        foo[i] <- MEAN[i+1]
    }
}
# and generate new names for the array
d <- data.frame(date=seq(from=as.Date("1901-01-01"), to=as.Date("2009-12-31"), by="+3 month"))
d$seas <- mkseas(d,"DJF") 
d$seasyear <- paste(format(d[,1],"%Y"), d$seas ,sep="")
names(foo)<-d$seasyear  # add right order colnames
MEAN <-foo

最后,这会产生季节性均值的时间序列。好吧,我觉得它太复杂了,我想周围有更简单的解决方案。

此外,这个解决方案在冬季 DJF 中还有一个非常大的问题:12 月到目前为止还没有从前一年中选择。这很容易解决(我猜),但使给定的方式更加复杂。

我真的希望周围有更好的想法!

【问题讨论】:

  • 此代码 sn-p 可能会有所帮助:dd &lt;- c(Sys.Date(), as.Date(c("2013-11-30", "2013-12-01"))); season_year &lt;- as.numeric(format(dd + 31, "%Y")).
  • 好点,snipplet 可能会有所帮助
  • 为了解决冬季的问题(在DJF中,D应该是n-1年的D),一个想法是创建一个“假”的年份列,其中包含当前年份的值每个月,除了 12 月,您使用 n+1。

标签: r time-series mean


【解决方案1】:

这就是你想要的?

# # create some data: daily values for three years
df <- data.frame(date = seq(from = as.Date("2007-01-01"),
                            to = as.Date("2009-12-31"),
                            by = "day"))
df$vals <- rnorm(nrow(df))

# add year
df$year <- format(df$date, "%Y")

# add season
df$seas <- mkseas(x = df, width = "DJF")

# calculate mean per season within each year
df2 <- aggregate(vals ~ seas + year, data = df, mean)

df2
#    seas year         vals
# 1   DJF 2007 -0.048407610
# 2   MAM 2007  0.086996842
# 3   JJA 2007  0.013864555
# 4   SON 2007 -0.081323367
# 5   DJF 2008  0.170887946
# 6   MAM 2008  0.147830260
# 7   JJA 2008  0.003008866
# 8   SON 2008 -0.057974215
# 9   DJF 2009 -0.043437437
# 10  MAM 2009 -0.048345979
# 11  JJA 2009  0.023860506
# 12  SON 2009 -0.060076870

因为mkseas 将日期转换为具有所需顺序的级别的季节性因素,所以在按年份和季节进行聚合后,顺序也是正确的。

【讨论】:

  • Henrik,这个看起来真的很漂亮/优雅!诚然,mkseas 也为 DJF 保留了正确的顺序。
  • 我终于添加了一个日期变量来绘制时间序列 >>> df2$date
  • 这不适用于月度数据(DJF 涵盖两年),作为答案添加了月度解决方案。
【解决方案2】:

如果您使用数字而不是字符串来表示月份和季节,这可能会更容易,至少一开始是这样。您可以通过简单的算术运算获得所需的季节,包括将 12 月作为下一年的一部分。

pdsi <- data.frame(date = seq(
            from=as.Date("1901-01-01"), 
            to=as.Date("2009-12-31"), 
            by="day"))
pdsi$scPDSI <- rnorm(nrow(pdsi),  mean=1, sd=1)
pdsi$mon<-mon(pdsi$date)+1
pdsi$seas<-floor((pdsi$mon %% 12)/3)+1
pdsi$year<-year(pdsi$date)+1900
pdsi$syear<-pdsi$year
pdsi$syear[pdsi$mon==12]<-pdsi$syear[pdsi$mon==12]+1

要计算季节性平均值,您可以这样做:

meanArray<-tapply(pdsi$scPDSI,list(year=pdsi$syear,seas=pdsi$seas),mean)

现在你有了

>head(meanArray)
      seas
year           1         2         3         4
  1901 1.0779676 1.0258306 1.1515175 0.9682434
  1902 0.9900312 0.8964994 1.1028336 1.0074296
  1903 0.9912233 0.9858088 1.1346901 1.0569518
  1904 0.7933653 1.1566892 1.1223454 0.8914211
  1905 1.1441863 1.1824074 0.9044940 0.8971485
  1906 0.9900826 0.9933909 0.9185972 0.8922987

如果你想要它作为一个平面数组,并有适当的名称,你首先进行转置,然后将数组展平,并添加名称

colnames(meanArray)<-c("DJF","MAM","JJA","SON")
meanArray<-t(meanArray)
MEAN<-array(meanArray)
names(MEAN)<-paste(colnames(meanArray)[col(meanArray)],rownames(meanArray)[row(meanArray)],sep="")

这会让你得到你想要的结果

> head(MEAN)
  1901DJF   1901MAM   1901JJA   1901SON   1902DJF   1902MAM 
1.0779676 1.0258306 1.1515175 0.9682434 0.9900312 0.8964994  

【讨论】:

    【解决方案3】:

    如前所述,可以有非常简单的解决方案(也发布了here)。我会使用 zooseas 包的组合按季节聚合,看起来像这样:

    library(zoo); library(seas)
    
    seasTS <- aggregate(dataTS, mkseas(x=time(dataTS),width="DJF"), sum)
    

    要每年执行此操作,只需按年循环 mkseas()。请给我加点语法糖的咖啡。

    干杯,

    亚当

    【讨论】:

      【解决方案4】:

      我遇到了同样的问题,但是对于月度数据,aggregate 多年来不允许 DJF 拆分。为了解决这个问题,您可以添加一个合成年份列,将 12 月的值分配给下一年。

      library(dplyr)
      library(seas)
      library(lubridate)
      
      df <- data.frame(yearmonth = c("187601", "187602", "187603", "187604", "187605", "187606", "187607","187608", "187609", "187610", "187611", "187612", "187701", "187702", "187703", "187704", "187705", "187706", "187707", "187708", "187709", "187710", "187711", "187712", "187801", "187802", "187803", "187804", "187805", "187806", "187807", "187808", "187809", "187810", "187811", "187812", "187901", "187902", "187903", "187904", "187905", "187906", "187907", "187908", "187909", "187910", "187911", "187912"), 
                       SOI = rnorm(n = 48, mean = 0, sd = 4))
      
      
      df %>% 
        mutate(yearmonth = lubridate::ymd(yearmonth, truncated = 1),
               year = year(yearmonth),
               month = month(yearmonth),
               seas = mkseas(yearmonth, width = "DJF"),
               year2 = ifelse(test = month == 12,
                              yes = year + 1,
                              no = year)) %>% 
        group_by(year2, seas) %>% 
        summarise(meanSOI = mean(SOI))
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 2020-01-31
        • 1970-01-01
        • 1970-01-01
        • 2013-07-25
        • 2010-12-15
        • 2011-03-22
        相关资源
        最近更新 更多