循环代码从 R 中的每日最大值和最小值模拟每小时温度

Question

我的最终目标是使用 R 来模拟从 1986 年到 2017 年的每日温度最大值和最小值的每小时温度。我已经成功地为单个日期的数据编写了代码，但是我无法跨应用此代码许多日期。

我从国家资源保护局 (NRCS) 获得了我的重点站点的每日温度数据： https://wcc.sc.egov.usda.gov/nwcc/site?sitenum=526

关注此处发布的模型：

Reicosky, D.S., Winkelman, L.J., Baker, J.M., Baker, D.G. 1989. 根据每日最小值和最大值计算的每小时气温的准确性。农业和森林气象学。 46:193-209

我编写了以下代码，非常适合对一天的每小时温度数据进行建模：

#create df for SINGLE DATE. 
#The actual data frame that I wish to model temperatures from will be exactly like this 
#but with 11,689 rows.

d8a <- data.frame(
  Day.of.Year = 213, 
  Date = as.Date("01-Aug-2011",format = "%d-%b-%Y"), 
  SunRise_decimal = 4.9, 
  Air.Temperature.Minimum..degC. = 8.0, 
  Air.Temperature.Maximum..degC. = 22.1
) 

#create matrix to serve as repository for modeled hourly temp data

OneDay <- data.frame(OneDay <- matrix(0, ncol = 0, nrow = 24))

hour <- OneDay$hour <- c(0:23)
rise <- OneDay$sunrise <- d8a$SunRise_decimal
tmax <- OneDay$tmax <- d8a$Air.Temperature.Maximum..degC.
tmin <- OneDay$tmin <- d8a$Air.Temperature.Minimum..degC.
tavg <- OneDay$tavg <- (OneDay$tmax + OneDay$tmin) / 2
peakhour <- OneDay$peakhour <- 14
amp <- OneDay$amp <- (OneDay$tmax - OneDay$tmin)/2

#Now for the actual modelling:

OneDay$tmod <- ifelse(hour < rise, tavg + amp * cos(pi * (hour + 10) / (10 + rise)), 
           ifelse(hour > peakhour, tavg + amp * cos(pi * (hour - peakhour) / (10 + rise)),
                    ifelse(hour >= rise, tavg - amp * cos(pi * (hour - rise) / (peakhour - rise)),
                                                               99999)))

plot(tmod ~ hour, data = OneDay, pch = 19, cex = 1.5, ylim = c(8,23), 
    main = "01 August 2011", las = 1, ylab = "Temp (C)", xlab = "Hour of Day")
lines(tmod ~ hour, data = OneDay)

最后，我的问题：

如何在由许多日期组成的数据框中的每个日期迭代此代码（或此代码的更高效版本）？

我知道最终的数据集会很大。（（31 年 * 每年 365 天 * 每天 24 小时）= 280,320 行）

Answer 1

一个非常简单的方法是 for 循环，我想你也可以用 apply 做一些事情，但我想这里一个循环就足够了，尤其是因为它只有 11000 次计算 (...)。

假设您的数据保存在数据帧 d8a 中

    OneDay<-list()
for(i in 1:nrow(d8a)){
OneDay[[i]] <- data.frame(OneDay[[i]] <- matrix(0, ncol = 8, nrow = 24))

hour <- OneDay[[i]][,1] <- c(0:23)
rise <- OneDay[[i]][,2] <- d8a$SunRise_decimal[i]
tmax <- OneDay[[i]][,3] <- d8a$Air.Temperature.Maximum..degC.[i]
tmin <- OneDay[[i]][,4] <- d8a$Air.Temperature.Minimum..degC.[i]
tavg <- OneDay[[i]][,5] <- (OneDay[[i]][,3] + OneDay[[i]][,4]) / 2
peakhour <- OneDay[[i]][,6] <- 14
amp <- OneDay[[i]][,7] <- (OneDay[[i]][,3] - OneDay[[i]][,4])/2

#Now for the actual modelling:

OneDay[[i]][,8] <- ifelse(hour < rise, tavg + amp * cos(pi * (hour + 10) / (10 + rise)), 
           ifelse(hour > peakhour, tavg + amp * cos(pi * (hour - peakhour) / (10 + rise)),
                    ifelse(hour >= rise, tavg - amp * cos(pi * (hour - rise) / (peakhour - rise)),
                                                               99999)))
}

这可能会让您更好地理解代码，因为它本质上是带有包装循环的代码。现在每天都会保存在一个单独的列表中，您可以稍后将它们合并或保持原样。

Answer 2

似乎 data.table 可以让这一切变得简单！

首先，将您的建模逻辑包含在一个函数中：

ModelHourly <- function(hour, rise, tmax, tmin) {
  peakhour <- 14
  tavg <- (tmax + tmin) / 2
  amp <- (tmax - tmin) / 2
  tmod <- ifelse(hour < rise, tavg + amp * cos(pi * (hour + 10) / (10 + rise)), 
         ifelse(hour > peakhour, tavg + amp * cos(pi * (hour - peakhour) / (10 + rise)),
                ifelse(hour >= rise, tavg - amp * cos(pi * (hour - rise) / (peakhour - rise)),
                       99999))) 
  return(tmod)
}

现在设置一个两天的示例数据集。

d8a <- data.frame(
  Day.of.Year = 213, 
  Date = as.Date("01-Aug-2011",format = "%d-%b-%Y"), 
  SunRise_decimal = 4.9, 
  Air.Temperature.Minimum..degC. = 8.0, 
  Air.Temperature.Maximum..degC. = 22.1
) 
d9a <- 
  data.frame(
    Day.of.Year = 214, 
    Date = as.Date("02-Aug-2011",format = "%d-%b-%Y"), 
    SunRise_decimal = 5.0, 
    Air.Temperature.Minimum..degC. = 7.0, 
    Air.Temperature.Maximum..degC. = 25.1
  ) 

dat <- rbind(d8a, d9a)

把它变成data.table

library('data.table')
dat <- as.data.table(dat)

现在我们需要将每一行复制 24 次并用 0:23 填充它。这似乎是概念上最简单的方法，但可能还有更巧妙的方法：

hourly <- dat[, .(hour=0:23), .(Date)]
dat <- merge(hourly, dat, by='Date')

如果您不熟悉 data.table，我所做的是创建一个新的 table (hourly)，其中有一个名为 "hour" 的列0:23，我每个 Date 都这样做。然后我们将它合并回原始数据 table Date 列。

现在只需调用您的函数即可！

dat[, modeled := ModelHourly(hour, SunRise_decimal, Air.Temperature.Maximum..degC., Air.Temperature.Minimum..degC.)]

如果你plot(dat$modeled)你会看到两条正弦曲线

循环代码从 R 中的每日最大值和最小值模拟每小时温度

Looping code to model hourly temperature from daily maxima and minima in R

iteration

modeling

loops

r

temperature