#library(lubridate) library(geosphere) # for the temperature pick the observations at midday, 13:00 # because the normal temperature is at midday TempBlindern <- TempBlindern[as.character(index(TempBlindern),"%H")=="13"] # wind at the same time as temperature FFBlindern <- FFBlindern[as.character(index(FFBlindern),"%H")=="13"] # pick the rainfall observed at the end of the day RR12Blindern <- RR12Blindern[as.character(index(RR12Blindern),"%H")=="19"] # pick the cloud coverage observed in the morning NNBlindern <- NNBlindern[as.character(index(NNBlindern),"%H")=="07"] # now convert the posix dates to regular date classes, so can do econometrics OT24Blindern <- xts(as.matrix(OT24Blindern),as.Date(index(OT24Blindern))) names (OT24Blindern) <- "OT24" RR12Blindern <- xts(as.matrix(RR12Blindern),as.Date(index(RR12Blindern))) names (RR12Blindern) <- "RR12" NNBlindern <- xts(as.matrix(NNBlindern), as.Date(index(NNBlindern))) names(NNBlindern) <- "NN" FFBlindern <- xts(as.matrix(FFBlindern), as.Date(index(FFBlindern))) names(FFBlindern) <- "FF" TempBlindern <- xts(as.matrix(TempBlindern),as.Date(index(TempBlindern))) names(TempBlindern) <- "TA" data <- merge(TempBlindern,normtempsBlindern) DiffMeanTemp <- data$TA - data$NormTemp names(DiffMeanTemp) <- "DiffMean" # the OT stuff needs to be converted a bit # first we lag by one day, since the observation is at 7 in the morning # the following day OT24Blindern <- lag(OT24Blindern,-1) # then we construct the variable that measures # the fraction of the day with sunshine length <- daylength(59,as.Date(index(OT24Blindern))) # Blindern is at 59 north FracDaySunny <- OT24Blindern/length names(FracDaySunny) <- "FracDaySunny" data <- merge(TempBlindern,FFBlindern,normtempsBlindern,all=FALSE) temp <- data$TA wind <- data$FF * (1000/3600) normtemp <- data$NormTemp #convert from m/s to km/h windchill <- 13.12 + 0.6215 * temp -11.37 *wind^0.16 + 0.3965 * temp * wind^0.16 # calculate an windchill adjusted normal temperature # assuming wind equal average of wind #meanwind <- mean(wind) meanwind <- 0 windchill_normal <- 13.12 + 0.6215 * normtemp - 11.37 * (meanwind)^0.16 + 0.3965 * normtemp * (meanwind)^0.16 DiffWindChill <- windchill - windchill_normal names(DiffWindChill) <- c("DiffWindChill") head(DiffWindChill)