Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

notes/warnings fixes in data.atmosphere package #3218

Merged
merged 30 commits into from
Mar 12, 2024
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
af9685e
fix notes/warnings in data.atmosphere
moki1202 Sep 15, 2023
7997829
mend
moki1202 Oct 15, 2023
f51d9df
revert accidental remote log file changes
moki1202 Nov 4, 2023
600c619
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Jan 2, 2024
3aa034e
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Jan 11, 2024
77fa8d2
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Jan 30, 2024
5dccb4a
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Feb 13, 2024
3acbf5a
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Feb 19, 2024
00cdd56
typo
infotroph Mar 4, 2024
c1203b6
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Mar 4, 2024
fd75c44
.data
infotroph Mar 4, 2024
0b0d243
remove dead code
infotroph Mar 4, 2024
5e9e282
only one fn in file now, so rename to match
infotroph Mar 4, 2024
a355e7c
Update modules/data.atmosphere/R/temporal.downscaling.R
mdietze Mar 4, 2024
0f7f810
stop importing all of dplyr and tidyselect
infotroph Mar 6, 2024
dcb4c60
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Mar 7, 2024
be499cd
Merge branch 'develop' into fix_warnings_atmosphere
infotroph Mar 7, 2024
913bd9e
fix always-true date comparison in met.process
infotroph Mar 7, 2024
7a30f98
not used
infotroph Mar 9, 2024
f571b0c
not sure why tests were failing to create outdir, but this is cleaner…
infotroph Mar 9, 2024
15181b9
eliminate direct import of magrittr
infotroph Mar 9, 2024
9f6b510
test updates in preparation to excise data.table
infotroph Mar 9, 2024
ab7285c
air_pressure not present in output
infotroph Mar 9, 2024
4feb208
unlisted import of .env
infotroph Mar 9, 2024
32b7afb
stop using data.table in load.cfmet
infotroph Mar 9, 2024
22d8e88
stop using data.table in cfmet downscaling
infotroph Mar 9, 2024
3a73092
update Rcheck_reference.log
infotroph Mar 9, 2024
a787406
biocro assumes met will be data.table; deferring any changes there fo…
infotroph Mar 9, 2024
691aa6b
not sure where temp_min and temp_max come from, but it is not df
infotroph Mar 9, 2024
f6b0857
try again to remove data.remote log changes
infotroph Mar 11, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
180 changes: 90 additions & 90 deletions modules/data.atmosphere/R/debias_met_regression.R

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion modules/data.atmosphere/R/download.Fluxnet2015.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
##' @title download.Fluxnet2015
##' @export
##' @param sitename the FLUXNET ID of the site to be downloaded, used as file name prefix.
##' The 'SITE_ID' field in \href{http://fluxnet.fluxdata.org//sites/site-list-and-pages/}{list of Ameriflux sites}
##' The 'SITE_ID' field in \href{https://fluxnet.org/sites/site-list-and-pages/}{list of Ameriflux sites}
##' @param outfolder location on disk where outputs will be stored
##' @param start_date the start date of the data to be downloaded. Format is YYYY-MM-DD (will only use the year part of the date)
##' @param end_date the end date of the data to be downloaded. Format is YYYY-MM-DD (will only use the year part of the date)
Expand Down
2 changes: 1 addition & 1 deletion modules/data.atmosphere/R/download.NEONmet.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
##'
##' @export
##' @param sitename the NEON ID of the site to be downloaded, used as file name prefix.
##' The 4-letter SITE code in \href{http://www.neonscience.org/science-design/field-sites/list}{list of NEON sites}
##' The 4-letter SITE code in \href{https://www.neonscience.org/science-design/field-sites/list}{list of NEON sites}
##' @param outfolder location on disk where outputs will be stored
##' @param start_date the start date of the data to be downloaded. Format is YYYY-MM-DD (will only use the year and month of the date)
##' @param end_date the end date of the data to be downloaded. Format is YYYY-MM-DD (will only use the year and month part of the date)
Expand Down
2 changes: 1 addition & 1 deletion modules/data.atmosphere/R/download.US_Wlef.R
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ download.US_Wlef <- function(start_date, end_date, timestep = 1) {
url <- paste0(base_url, year,"/flux_", year, ".txt") #Build proper url
PEcAn.logger::logger.info(paste0("Reading data for year ", year))
print(url)
influx <- tryCatch(read.table(url, header = T, sep = ""), error=function(e) {NULL}, warning=function(e) {NULL})
influx <- tryCatch(utils::read.table(url, header = T, sep = ""), error=function(e) {NULL}, warning=function(e) {NULL})
if (is.null(influx)) { #Error encountered in data fetching.
PEcAn.logger::logger.warn(paste0("Data not avaliable for year ", year, ". All values for ", year, " will be NA."))
# Determine the number of days in the year
Expand Down
4 changes: 2 additions & 2 deletions modules/data.atmosphere/R/extract_ERA5.R
Original file line number Diff line number Diff line change
Expand Up @@ -113,13 +113,13 @@ extract.nc.ERA5 <-
# send out as xts object
xts::xts(all.data.point, order.by = timestamp)
}) %>%
setNames(paste0("ERA_ensemble_", ensemblesN))
stats::setNames(paste0("ERA_ensemble_", ensemblesN))

#Merge mean and the speard
return(point.data)

}) %>%
setNames(years)
stats::setNames(years)


# The order of one.year.out is year and then Ens - Mainly because of the spead / I wanted to touch each file just once.
Expand Down
12 changes: 6 additions & 6 deletions modules/data.atmosphere/R/met.process.R
Original file line number Diff line number Diff line change
Expand Up @@ -229,24 +229,24 @@ met.process <- function(site, input_met, start_date, end_date, model,
cf.id <- raw.id <- db.file
}else{
# I did this bc dbfile.input.check does not cover the between two time periods situation
mimetypeid <- get.id(table = "mimetypes", colnames = "type_string",
mimetypeid <- PEcAn.DB::get.id(table = "mimetypes", colnames = "type_string",
values = "application/x-netcdf", con = con)

formatid <- get.id(table = "formats", colnames = c("mimetype_id", "name"),
formatid <- PEcAn.DB::get.id(table = "formats", colnames = c("mimetype_id", "name"),
values = c(mimetypeid, "CF Meteorology"), con = con)

machine.id <- get.id(table = "machines", "hostname", PEcAn.remote::fqdn(), con)
machine.id <- PEcAn.DB::get.id(table = "machines", "hostname", PEcAn.remote::fqdn(), con)
# Fidning the tiles.
raw.tiles <- tbl(con, "inputs") %>%
filter(
site_id == register$ParentSite,
.data$site_id == register$ParentSite,
start_date >= start_date,
end_date <= end_date,
infotroph marked this conversation as resolved.
Show resolved Hide resolved
format_id == formatid
.data$format_id == formatid
) %>%
filter(grepl(met, "name")) %>%
inner_join(tbl(con, "dbfiles"), by = c('id' = 'container_id')) %>%
filter(machine_id == machine.id) %>%
filter(.data$machine_id == machine.id) %>%
collect()

cf.id <- raw.id <- list(input.id = raw.tiles$id.x, dbfile.id = raw.tiles$id.y)
Expand Down
4 changes: 4 additions & 0 deletions modules/data.atmosphere/R/met_functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,10 @@ get.latlonbox <- function(lati, loni, Lat = Lat, Lon = Lon) {

get.cruncep <- function(lat, lon, start.date, end.date) {
result <- load.cfmet(lat, lon)
Lat <- ncdf4::ncvar_get(nc = met.nc, varid = "lat")
Lon <- ncdf4::ncvar_get(nc = met.nc, varid = "lon")
lati <- which.min(abs(Lat - lat))
infotroph marked this conversation as resolved.
Show resolved Hide resolved

hourly.result <- cruncep_hourly(result, lat = Lat[lati])
weather <- cruncep_dt2weather(hourly.result)
return(weather)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ met_temporal_downscale.Gaussian_ensemble <- function(in.path, in.prefix, outfold
if (highday > reso_len) {
highday <- reso_len
}
dwnsc_day <- rand_vect_cont(div, source$precipitation_flux[x], sd = sd(train$precipitation_flux[lowday:highday]))
dwnsc_day <- rand_vect_cont(div, source$precipitation_flux[x], sd = stats::sd(train$precipitation_flux[lowday:highday]))
precip <- append(precip, dwnsc_day)
}
df$precipitation_flux <- precip
Expand All @@ -196,7 +196,7 @@ met_temporal_downscale.Gaussian_ensemble <- function(in.path, in.prefix, outfold
}
dwnsc_day <- vector()
for (n in seq_len(div)) {
dwnsc_day[n] <- rnorm(1, mean = sour[x], sd = sd(a[lowday:highday]))
dwnsc_day[n] <- stats::rnorm(1, mean = sour[x], sd = stats::sd(a[lowday:highday]))
}
train_vec <- append(train_vec, dwnsc_day)
}
Expand Down Expand Up @@ -330,7 +330,7 @@ met_temporal_downscale.Gaussian_ensemble <- function(in.path, in.prefix, outfold
if (precip[i] == 0) {
p <- 2
}
hdry[i] <- A * p * (1 - exp(-1 * bmlist[i] * ((temp_max[i] - temp_min[i])^C)))
hdry[i] <- A * p * (1 - exp(-1 * bmlist[i] * ((df$temp_max[i] - df$temp_min[i])^C)))
infotroph marked this conversation as resolved.
Show resolved Hide resolved
}
hdry[hdry < 0] <- 0
swflux <- hdry * I
Expand Down
18 changes: 9 additions & 9 deletions modules/data.atmosphere/R/metgapfill.NOAA_GEFS.R
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ metgapfill.NOAA_GEFS <- function(in.prefix, in.path, outfolder, start_date, end_
for (i in 1:nrow(allvars)) {
for (j in k:ncol(allvars)) {
if (is.na(allvars[i,j])) {
allvars[i,j] = sample(na.omit(allvars[i,seq(j, 1, by = -4)]), 1)
allvars[i,j] = sample(stats::na.omit(allvars[i,seq(j, 1, by = -4)]), 1)
}
}
}
Expand Down Expand Up @@ -114,22 +114,22 @@ metgapfill.NOAA_GEFS <- function(in.prefix, in.path, outfolder, start_date, end_
# Unfortunately, R is picky, and the data.frame[['var_as_string']] notation doesn't work
# for the lm function; only the $ notation does, hence this if/else if section.
if (dependent_vars[i] == "specific_humidity") {
reg <- lm(fitted.data$specific_humidity ~.,fitted.data)
reg <- stats::lm(fitted.data$specific_humidity ~.,fitted.data)
} else if (dependent_vars[i] == "surface_downwelling_longwave_flux_in_air") {
reg <- lm(fitted.data$surface_downwelling_longwave_flux_in_air ~.,fitted.data)
reg <- stats::lm(fitted.data$surface_downwelling_longwave_flux_in_air ~.,fitted.data)
} else if (dependent_vars[i] == "surface_downwelling_shortwave_flux_in_air") {
reg <- lm(fitted.data$surface_downwelling_shortwave_flux_in_air ~.,fitted.data)
reg <- stats::lm(fitted.data$surface_downwelling_shortwave_flux_in_air ~.,fitted.data)
} else if (dependent_vars[i] == "air_pressure") {
reg <- lm(fitted.data$air_pressure ~.,fitted.data)
reg <- stats::lm(fitted.data$air_pressure ~.,fitted.data)
} else if (dependent_vars[i] == "eastward_wind") {
reg <- lm(fitted.data$eastward_wind ~.,fitted.data)
reg <- stats::lm(fitted.data$eastward_wind ~.,fitted.data)
} else if (dependent_vars[i] == "northward_wind") {
reg <- lm(fitted.data$northward_wind ~.,fitted.data)
reg <- stats::lm(fitted.data$northward_wind ~.,fitted.data)
}else if (dependent_vars[i] == "wind_speed") {
reg <- lm(fitted.data$wind_speed ~.,fitted.data)
reg <- stats::lm(fitted.data$wind_speed ~.,fitted.data)
}

prediction <- predict(reg, fitted.data)
prediction <- stats::predict(reg, fitted.data)

# Update the values in the data frame
for (j in 1:length(prediction)) {
Expand Down
2 changes: 1 addition & 1 deletion modules/data.atmosphere/R/metutils.R
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ qcshum <- function(x) {
##' converting specific humidity into relative humidity
##' NCEP surface flux data does not have RH
##' from Bolton 1980 Teh computation of Equivalent Potential Temperature
##' \url{http://www.eol.ucar.edu/projects/ceop/dm/documents/refdata_report/eqns.html}
##' \url{https://archive.eol.ucar.edu/projects/ceop/dm/documents/refdata_report/eqns.html}
##' @title qair2rh
##' @param qair specific humidity, dimensionless (e.g. kg/kg) ratio of water mass / total air mass
##' @param temp degrees C
Expand Down
19 changes: 12 additions & 7 deletions modules/data.atmosphere/R/tdm_lm_ensemble_sims.R
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ lm_ensemble_sims <- function(dat.mod, n.ens, path.model, direction.filter, lags.
# shortwave is different because we only want to model daylight
if (v == "surface_downwelling_shortwave_flux_in_air") {
# Finding which days have measurable light
thresh.swdown <- quantile(dat.train$surface_downwelling_shortwave_flux_in_air[dat.train$surface_downwelling_shortwave_flux_in_air > 0], 0.05)
thresh.swdown <- stats::quantile(dat.train$surface_downwelling_shortwave_flux_in_air[dat.train$surface_downwelling_shortwave_flux_in_air > 0], 0.05)


hrs.day <- unique(dat.train$time[dat.train$time$DOY == day.now &
Expand Down Expand Up @@ -227,8 +227,13 @@ lm_ensemble_sims <- function(dat.mod, n.ens, path.model, direction.filter, lags.
}

# Load the saved model
load(file.path(path.model, v, paste0("model_", v, "_", day.now,
".Rdata")))
model.file <- file.path(path.model, v, paste0("model_", v, "_", day.now,
".Rdata"))
if(file.exists(model.file)) {
env = new.env()
load(model.file, env = env)
infotroph marked this conversation as resolved.
Show resolved Hide resolved
mod.save <- env$mod.save
}

# Pull coefficients (betas) from our saved matrix

Expand Down Expand Up @@ -340,18 +345,18 @@ lm_ensemble_sims <- function(dat.mod, n.ens, path.model, direction.filter, lags.


filter.mean <- mean(dat.filter, na.rm=T)
filter.sd <- sd(dat.filter, na.rm=T)
filter.sd <- stats::sd(dat.filter, na.rm=T)
} else {

if(v %in% vars.sqrt){
filter.mean <- mean(c(dat.pred^2, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- sd(c(dat.pred^2, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- stats::sd(c(dat.pred^2, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
} else if(v %in% vars.log){
filter.mean <- mean(c(exp(dat.pred), utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- sd(c(exp(dat.pred), utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- stats::sd(c(exp(dat.pred), utils::stack(dat.sim[[v]])[,1]), na.rm=T)
} else {
filter.mean <- mean(c(dat.pred, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- sd(c(dat.pred, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
filter.sd <- stats::sd(c(dat.pred, utils::stack(dat.sim[[v]])[,1]), na.rm=T)
}
}

Expand Down
58 changes: 29 additions & 29 deletions modules/data.atmosphere/R/tdm_model_train.R
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL
dat.subset$year <- as.ordered(dat.subset$year)
if (v == "air_temperature") {

mod.doy <- lm(air_temperature ~ as.ordered(hour) * air_temperature_max.day *
mod.doy <- stats::lm(air_temperature ~ as.ordered(hour) * air_temperature_max.day *
(lag.air_temperature + lag.air_temperature_min + air_temperature_min.day) +
as.ordered(hour) * air_temperature_min.day * next.air_temperature_max -
1 - as.ordered(hour) - lag.air_temperature - lag.air_temperature_min -
Expand All @@ -41,15 +41,15 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL
hrs.day <- unique(dat.subset[dat.subset$surface_downwelling_shortwave_flux_in_air >
threshold, "hour"])

mod.doy <- lm(surface_downwelling_shortwave_flux_in_air ~
mod.doy <- stats::lm(surface_downwelling_shortwave_flux_in_air ~
as.ordered(hour) * surface_downwelling_shortwave_flux_in_air.day -
1 - surface_downwelling_shortwave_flux_in_air.day -
as.ordered(hour), data = dat.subset[dat.subset$hour %in%
hrs.day, ]) ###
}

if (v == "surface_downwelling_longwave_flux_in_air") {
mod.doy <- lm(sqrt(surface_downwelling_longwave_flux_in_air) ~
mod.doy <- stats::lm(sqrt(surface_downwelling_longwave_flux_in_air) ~
as.ordered(hour) * surface_downwelling_longwave_flux_in_air.day *
(lag.surface_downwelling_longwave_flux_in_air + next.surface_downwelling_longwave_flux_in_air) -
as.ordered(hour) - 1 - lag.surface_downwelling_longwave_flux_in_air -
Expand All @@ -66,25 +66,25 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL
# fraction of precip occuring in each hour we're going to estimate the
# probability distribution of rain occuring in a given hour
dat.subset$rain.prop <- dat.subset$precipitation_flux/(dat.subset$precipitation_flux.day)
mod.doy <- lm(rain.prop ~ as.ordered(hour) - 1 , data = dat.subset)
mod.doy <- stats::lm(rain.prop ~ as.ordered(hour) - 1 , data = dat.subset)
}

if (v == "air_pressure") {
mod.doy <- lm(air_pressure ~ as.ordered(hour) * (air_pressure.day +
mod.doy <- stats::lm(air_pressure ~ as.ordered(hour) * (air_pressure.day +
lag.air_pressure + next.air_pressure) - as.ordered(hour) -
1 - air_pressure.day - lag.air_pressure - next.air_pressure,
data = dat.subset)
}

if (v == "specific_humidity") {
mod.doy <- lm(log(specific_humidity) ~ as.ordered(hour) *
mod.doy <- stats::lm(log(specific_humidity) ~ as.ordered(hour) *
specific_humidity.day * (lag.specific_humidity + next.specific_humidity +
air_temperature_max.day) - as.ordered(hour) - 1 - air_temperature_max.day,
data = dat.subset)
}

if (v == "wind_speed") {
mod.doy <- lm(sqrt(wind_speed) ~ as.ordered(hour) * wind_speed.day *
mod.doy <- stats::lm(sqrt(wind_speed) ~ as.ordered(hour) * wind_speed.day *
(lag.wind_speed + next.wind_speed) - as.ordered(hour) -
1 - wind_speed.day - lag.wind_speed - next.wind_speed -
wind_speed.day * lag.wind_speed - wind_speed.day * next.wind_speed,
Expand All @@ -104,13 +104,13 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL
# coefficients that we can pull from without needing to do this step
# every day
if(n.beta>1){
mod.coef <- coef(mod.doy)
mod.cov <- vcov(mod.doy)
mod.coef <- stats::coef(mod.doy)
mod.cov <- stats::vcov(mod.doy)
piv <- as.numeric(which(!is.na(mod.coef)))
Rbeta <- MASS::mvrnorm(n = n.beta, mod.coef[piv], mod.cov[piv,piv])
} else {
Rbeta <- matrix(coef(mod.doy), nrow=1)
colnames(Rbeta) <- names(coef(mod.doy))
Rbeta <- matrix(stats::coef(mod.doy), nrow=1)
colnames(Rbeta) <- names(stats::coef(mod.doy))
}

list.out <- list(model = mod.doy, betas = Rbeta)
Expand All @@ -120,63 +120,63 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL
if (resids == TRUE) {
if (v == "air_temperature") {
dat.subset[!is.na(dat.subset$lag.air_temperature) & !is.na(dat.subset$next.air_temperature_max),
"resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * (air_temperature_max.day *
"resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * (air_temperature_max.day *
air_temperature_min.day) - 1, data = dat.subset[!is.na(dat.subset$lag.air_temperature),
])
}

if (v == "surface_downwelling_shortwave_flux_in_air") {
dat.subset[dat.subset$hour %in% hrs.day, "resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * surface_downwelling_shortwave_flux_in_air.day -
dat.subset[dat.subset$hour %in% hrs.day, "resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * surface_downwelling_shortwave_flux_in_air.day -
1, data = dat.subset[dat.subset$hour %in% hrs.day,
])
}

if (v == "surface_downwelling_longwave_flux_in_air") {
dat.subset[!is.na(dat.subset$lag.surface_downwelling_longwave_flux_in_air) &
!is.na(dat.subset$next.surface_downwelling_longwave_flux_in_air),
"resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * surface_downwelling_longwave_flux_in_air.day -
"resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * surface_downwelling_longwave_flux_in_air.day -
1, data = dat.subset[, ])
}

if (v == "precipitation_flux") {
dat.subset[, "resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * precipitation_flux.day -
dat.subset[, "resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * precipitation_flux.day -
1, data = dat.subset[, ])
}

if (v == "air_pressure") {
dat.subset[!is.na(dat.subset$lag.air_pressure) & !is.na(dat.subset$next.air_pressure),
"resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * air_pressure.day -
"resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * air_pressure.day -
1, data = dat.subset[, ])
}

if (v == "specific_humidity") {
dat.subset[!is.na(dat.subset$lag.specific_humidity) &
!is.na(dat.subset$next.specific_humidity), "resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * specific_humidity.day -
!is.na(dat.subset$next.specific_humidity), "resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * specific_humidity.day -
1, data = dat.subset[, ])
}

if (v == "wind_speed") {
dat.subset[!is.na(dat.subset$lag.wind_speed) & !is.na(dat.subset$next.wind_speed),
"resid"] <- resid(mod.doy)
resid.model <- lm(resid ~ as.ordered(hour) * wind_speed.day -
"resid"] <- stats::resid(mod.doy)
resid.model <- stats::lm(resid ~ as.ordered(hour) * wind_speed.day -
1, data = dat.subset[, ])
}

if(n.beta>1){
res.coef <- coef(resid.model)
res.cov <- vcov(resid.model)
res.coef <- stats::coef(resid.model)
res.cov <- stats::vcov(resid.model)
res.piv <- as.numeric(which(!is.na(res.coef)))

beta.resid <- MASS::mvrnorm(n = n.beta, res.coef[res.piv], res.cov)
} else {
beta.resid <- matrix(coef(resid.model), nrow=1)
colnames(beta.resid) <- names(coef(mod.doy))
beta.resid <- matrix(stats::coef(resid.model), nrow=1)
colnames(beta.resid) <- names(stats::coef(mod.doy))
}


Expand Down
Loading
Loading