diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv index 2ffb26841a7..1f05a2b92f1 100644 --- a/docker/depends/pecan_package_dependencies.csv +++ b/docker/depends/pecan_package_dependencies.csv @@ -33,7 +33,6 @@ "data.table","*","base/visualization","Imports",FALSE "data.table","*","models/biocro","Imports",FALSE "data.table","*","models/ldndc","Imports",FALSE -"data.table","*","modules/data.atmosphere","Imports",FALSE "data.table","*","modules/data.remote","Suggests",FALSE "dataone","*","modules/data.land","Suggests",FALSE "datapack","*","modules/data.land","Imports",FALSE @@ -55,10 +54,10 @@ "dplyr","*","models/stics","Imports",FALSE "dplyr","*","modules/assim.sequential","Imports",FALSE "dplyr","*","modules/benchmark","Imports",FALSE -"dplyr","*","modules/data.atmosphere","Imports",FALSE "dplyr","*","modules/data.land","Imports",FALSE "dplyr","*","modules/data.remote","Suggests",FALSE "dplyr","*","modules/uncertainty","Imports",FALSE +"dplyr",">= 0.8.1","modules/data.atmosphere","Imports",FALSE "dplyr",">= 1.1.2","base/db","Imports",FALSE "ellipse","*","modules/assim.batch","Imports",FALSE "emdbook","*","modules/assim.sequential","Suggests",FALSE @@ -157,7 +156,6 @@ "magrittr","*","models/ed","Imports",FALSE "magrittr","*","modules/assim.sequential","Imports",FALSE "magrittr","*","modules/benchmark","Imports",FALSE -"magrittr","*","modules/data.atmosphere","Imports",FALSE "magrittr","*","modules/data.land","Imports",FALSE "magrittr","*","modules/data.remote","Imports",FALSE "markdown","*","modules/allometry","Suggests",FALSE @@ -178,7 +176,6 @@ "methods","*","modules/allometry","Imports",FALSE "methods","*","modules/assim.batch","Imports",FALSE "methods","*","modules/assim.sequential","Suggests",FALSE -"methods","*","modules/data.atmosphere","Depends",FALSE "methods","*","modules/emulator","Imports",FALSE "mgcv","*","modules/data.atmosphere","Imports",FALSE "minpack.lm","*","modules/rtm","Suggests",FALSE @@ -646,6 +643,7 @@ "withr","*","models/ed","Suggests",FALSE "withr","*","models/sibcasa","Suggests",FALSE "withr","*","modules/allometry","Suggests",FALSE +"withr","*","modules/data.atmosphere","Suggests",FALSE "XML","*","base/db","Imports",FALSE "XML","*","base/workflow","Imports",FALSE "XML","*","models/biocro","Imports",FALSE diff --git a/models/biocro/R/met2model.BIOCRO.R b/models/biocro/R/met2model.BIOCRO.R index 638a1c7110b..dee2e84566f 100644 --- a/models/biocro/R/met2model.BIOCRO.R +++ b/models/biocro/R/met2model.BIOCRO.R @@ -142,6 +142,11 @@ met2model.BIOCRO <- function(in.path, in.prefix, outfolder, overwrite = FALSE, ##' @author David LeBauer cf2biocro <- function(met, longitude = NULL, zulu2solarnoon = FALSE) { + if (!data.table::is.data.table(met)) { + met <- data.table::copy(met) + data.table::setDT(met) + } + if ((!is.null(longitude)) & zulu2solarnoon) { solarnoon_offset <- PEcAn.utils::ud_convert(longitude/360, "day", "minute") met[, `:=`(solardate = met$date + lubridate::minutes(solarnoon_offset))] diff --git a/modules/data.atmosphere/DESCRIPTION b/modules/data.atmosphere/DESCRIPTION index e1a4dd75f41..faa34789aa3 100644 --- a/modules/data.atmosphere/DESCRIPTION +++ b/modules/data.atmosphere/DESCRIPTION @@ -20,22 +20,18 @@ Description: The Predictive Ecosystem Carbon Analyzer (PEcAn) is a scientific package converts climate driver data into a standard format for models integrated into PEcAn. As a standalone package, it provides an interface to access diverse climate data sets. -Depends: - methods Imports: abind (>= 1.4.5), amerifluxr, arrow, curl, - data.table, - dplyr, + dplyr (>= 0.8.1), geonames (> 0.998), ggplot2, glue, httr, jsonlite, lubridate (>= 1.6.0), - magrittr, MASS, mgcv, ncdf4 (>= 1.15), @@ -71,7 +67,8 @@ Suggests: parallel, PEcAn.settings, progress, - reticulate + reticulate, + withr Remotes: github::adokter/suntools, github::chuhousen/amerifluxr, diff --git a/modules/data.atmosphere/NAMESPACE b/modules/data.atmosphere/NAMESPACE index 9dd14faa440..ad08797ea15 100644 --- a/modules/data.atmosphere/NAMESPACE +++ b/modules/data.atmosphere/NAMESPACE @@ -110,8 +110,7 @@ export(temporal_downscale_half_hour) export(upscale_met) export(wide2long) export(write_noaa_gefs_netcdf) -import(dplyr) -import(tidyselect) -importFrom(magrittr,"%>%") +importFrom(dplyr,"%>%") importFrom(rlang,.data) +importFrom(rlang,.env) importFrom(sf,st_crs) diff --git a/modules/data.atmosphere/R/ERA5_met_process.R b/modules/data.atmosphere/R/ERA5_met_process.R index 9de54deb5c2..de4835796c7 100644 --- a/modules/data.atmosphere/R/ERA5_met_process.R +++ b/modules/data.atmosphere/R/ERA5_met_process.R @@ -10,7 +10,7 @@ #' @export #' #' @author Dongchen Zhang -#' @importFrom magrittr %>% +#' @importFrom dplyr %>% #' ERA5_met_process <- function(settings, in.path, out.path, write.db=FALSE, write = TRUE){ #Initialize the multicore computation. diff --git a/modules/data.atmosphere/R/GEFS_helper_functions.R b/modules/data.atmosphere/R/GEFS_helper_functions.R index c63768daf11..754580ae0da 100644 --- a/modules/data.atmosphere/R/GEFS_helper_functions.R +++ b/modules/data.atmosphere/R/GEFS_helper_functions.R @@ -533,7 +533,6 @@ process_gridded_noaa_download <- function(lat_list, #' @param overwrite, logical stating to overwrite any existing output_file #' @param hr time step in hours of temporal downscaling (default = 1) #' @importFrom rlang .data -#' @import tidyselect #' #' @author Quinn Thomas #' diff --git a/modules/data.atmosphere/R/debias_met_regression.R b/modules/data.atmosphere/R/debias_met_regression.R index 9589777bf59..bae4dd6f036 100644 --- a/modules/data.atmosphere/R/debias_met_regression.R +++ b/modules/data.atmosphere/R/debias_met_regression.R @@ -66,8 +66,6 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU outfolder, yrs.save=NULL, ens.name, ens.mems=NULL, force.sanity=TRUE, sanity.tries=25, sanity.sd=8, lat.in, lon.in, save.diagnostics=TRUE, path.diagnostics=NULL, parallel = FALSE, n.cores = NULL, overwrite = TRUE, verbose = FALSE) { - library(MASS) - library(mgcv) set.seed(seed) @@ -190,9 +188,9 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU } # end y loop # Hard-coding in some sort of max for precipitaiton - rain.max <- max(train.data$precipitation_flux) + sd(train.data$precipitation_flux) - rainless.min <- ifelse(min(rainless)-sd(rainless)>=0, min(rainless)-sd(rainless), max(min(rainless)-sd(rainless)/2, 0)) - rainless.max <- ifelse(max(rainless)+sd(rainless)<=365, max(rainless)+sd(rainless), min(max(rainless)+sd(rainless)/2, 365)) + rain.max <- max(train.data$precipitation_flux) + stats::sd(train.data$precipitation_flux) + rainless.min <- ifelse(min(rainless)-stats::sd(rainless)>=0, min(rainless)-stats::sd(rainless), max(min(rainless)-stats::sd(rainless)/2, 0)) + rainless.max <- ifelse(max(rainless)+stats::sd(rainless)<=365, max(rainless)+stats::sd(rainless), min(max(rainless)+stats::sd(rainless)/2, 365)) } # ------------- @@ -216,7 +214,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # adjustment & anomaly as fraction of total annual precipitation if(v == "precipitation_flux"){ # Find total annual preciptiation - precip.ann <- aggregate(met.train$Y, by=met.train[,c("year", "ind")], FUN=sum) + precip.ann <- stats::aggregate(met.train$Y, by=met.train[,c("year", "ind")], FUN=sum) names(precip.ann)[3] <- "Y.tot" met.train <- merge(met.train, precip.ann, all=T) @@ -224,7 +222,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU } # Aggregate to get rid of years so that we can compare climatic means; bring in covariance among climatic predictors - dat.clim <- aggregate(met.train[,"Y"], by=met.train[,c("doy", "ind")], FUN=mean) + dat.clim <- stats::aggregate(met.train[,"Y"], by=met.train[,c("doy", "ind")], FUN=mean) # dat.clim[,v] <- 1 names(dat.clim)[3] <- "Y" # ----- @@ -241,7 +239,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # met.src[,v] <- if(v=="precipitation_flux"){ - src.ann <- aggregate(met.src$X, by=met.src[,c("year", "ind.src")], FUN=sum) + src.ann <- stats::aggregate(met.src$X, by=met.src[,c("year", "ind.src")], FUN=sum) names(src.ann)[3] <- "X.tot" met.src <- merge(met.src, src.ann, all.x=T) @@ -281,7 +279,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU met.src[,v] <- 0 # Aggregate to get rid of years so that we can compare climatic means - clim.src <- aggregate(met.src[met.src$year %in% yrs.overlap,c("X", vars.debias)], + clim.src <- stats::aggregate(met.src[met.src$year %in% yrs.overlap,c("X", vars.debias)], by=met.src[met.src$year %in% yrs.overlap,c("doy", "ind", "ind.src")], FUN=mean, na.rm=T) clim.src[,vars.debias[!vars.debias %in% names(dat.out)]] <- 0 @@ -350,29 +348,29 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # summary(mod.bias) # Saving the mean predicted & residuals - dat.clim[dat.clim$ind == ind, "pred"] <- predict(mod.bias) - dat.clim[dat.clim$ind == ind, "resid"] <- resid(mod.bias) + dat.clim[dat.clim$ind == ind, "pred"] <- stats::predict(mod.bias) + dat.clim[dat.clim$ind == ind, "resid"] <- stats::resid(mod.bias) # summary(dat.clim) # Storing the model residuals to add in some extra error - resid.bias <- resid(mod.bias) + resid.bias <- stats::resid(mod.bias) # # Checking the residuals to see if we can assume normality # plot(resid ~ pred, data=dat.clim); abline(h=0, col="red") # plot(resid ~ doy, data=dat.clim); abline(h=0, col="red") # hist(dat.clim$resid) - met.src [met.src $ind == ind, "pred"] <- predict(mod.bias, newdata=met.src [met.src $ind == ind, ]) - met.train[met.train$ind == ind, "pred"] <- predict(mod.bias, newdata=met.train[met.train$ind == ind, ]) + met.src [met.src $ind == ind, "pred"] <- stats::predict(mod.bias, newdata=met.src [met.src $ind == ind, ]) + met.train[met.train$ind == ind, "pred"] <- stats::predict(mod.bias, newdata=met.train[met.train$ind == ind, ]) # For Precip we need to bias-correct the total annual preciptiation + seasonal distribution if(v == "precipitation_flux"){ - mod.ann <- lm(Y.tot ~ X.tot , data=dat.ann[dat.ann$ind==ind,]) + mod.ann <- stats::lm(Y.tot ~ X.tot , data=dat.ann[dat.ann$ind==ind,]) # summary(mod.ann) - dat.ann[dat.ann$ind==ind,"pred.ann"] <- predict(mod.ann) - dat.ann[dat.ann$ind==ind,"resid.ann"] <- resid(mod.ann) + dat.ann[dat.ann$ind==ind,"pred.ann"] <- stats::predict(mod.ann) + dat.ann[dat.ann$ind==ind,"resid.ann"] <- stats::resid(mod.ann) - met.src[met.src$ind==ind,"pred.ann"] <- predict(mod.ann, newdata=met.src[met.src$ind==ind,]) + met.src[met.src$ind==ind,"pred.ann"] <- stats::predict(mod.ann, newdata=met.src[met.src$ind==ind,]) } # --------- @@ -394,8 +392,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU met.train[met.train$ind==ind,"anom.train"] <- met.train[met.train$ind==ind,"X"] met.src[met.src$ind==ind, "anom.raw"] <- met.src[met.src$ind==ind, "X"] } else { - met.train[met.train$ind==ind,"anom.train"] <- resid(anom.train) - met.src[met.src$ind==ind, "anom.raw"] <- met.src[met.src$ind==ind, "X"] - predict(anom.src, newdata=met.src[met.src$ind==ind, ]) + met.train[met.train$ind==ind,"anom.train"] <- stats::resid(anom.train) + met.src[met.src$ind==ind, "anom.raw"] <- met.src[met.src$ind==ind, "X"] - stats::predict(anom.src, newdata=met.src[met.src$ind==ind, ]) } # par(mfrow=c(2,1)) # plot(anom.train~doy, data=met.train) @@ -414,8 +412,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU anom.train2 <- mgcv::gam(Q ~ s(doy, k=6), data=met.train[met.train$ind==ind,]) anom.src2 <- mgcv::gam(Q ~ s(doy, k=6), data=met.src[met.src$year %in% yrs.overlap & met.src$ind==ind,]) - met.train[met.train$ind==ind, paste0(j, ".anom")] <- resid(anom.train2) - met.src[met.src$ind==ind, paste0(j, ".anom")] <- met.src[met.src$ind==ind,"Q"] - predict(anom.src2, newdata=met.src[met.src$ind==ind,]) + met.train[met.train$ind==ind, paste0(j, ".anom")] <- stats::resid(anom.train2) + met.src[met.src$ind==ind, paste0(j, ".anom")] <- met.src[met.src$ind==ind,"Q"] - stats::predict(anom.src2, newdata=met.src[met.src$ind==ind,]) rm(anom.train2, anom.src2) } @@ -445,7 +443,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # abline(lm(anom.train ~ anom.raw, data=dat.anom), col="red", lty="dashed") # Modeling in the predicted value from mod.bias - dat.anom$pred <- predict(mod.bias, newdata=dat.anom) + dat.anom$pred <- stats::predict(mod.bias, newdata=dat.anom) if (v %in% c("air_temperature", "air_temperature_maximum", "air_temperature_minimum")){ # ** We want to make sure we do these first ** @@ -506,16 +504,16 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # summary(mod.anom) # plot(mod.anom, pages=1) # pred.anom <- predict(mod.anom) - resid.anom <- resid(mod.anom) + resid.anom <- stats::resid(mod.anom) # --------- # -------- # Predicting a bunch of potential posteriors over the full dataset # -------- # Get the model coefficients - coef.gam <- coef(mod.bias) - coef.anom <- coef(mod.anom) - if(v == "precipitation_flux") coef.ann <- coef(mod.ann) + coef.gam <- stats::coef(mod.bias) + coef.anom <- stats::coef(mod.anom) + if(v == "precipitation_flux") coef.ann <- stats::coef(mod.ann) # Setting up a case where if sanity checks fail, we pull more ensemble members n.new <- 1 @@ -531,11 +529,11 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # I think the anomalies might be problematic, so lets get way more betas than we need and trim the distribution # set.seed=42 if(n.ens==1 | uncert.prop=="mean"){ - Rbeta <- matrix(coef(mod.bias), ncol=length(coef(mod.bias))) + Rbeta <- matrix(stats::coef(mod.bias), ncol=length(stats::coef(mod.bias))) } else { - Rbeta <- matrix(MASS::mvrnorm(n=n.new, coef(mod.bias), vcov(mod.bias)), ncol=length(coef(mod.bias))) + Rbeta <- matrix(MASS::mvrnorm(n=n.new, stats::coef(mod.bias), stats::vcov(mod.bias)), ncol=length(stats::coef(mod.bias))) } - dimnames(Rbeta)[[2]] <- names(coef(mod.bias)) + dimnames(Rbeta)[[2]] <- names(stats::coef(mod.bias)) # # Filter our betas to remove outliers # ci.beta <- matrix(apply(Rbeta, 2, quantile, c(0.01, 0.99)), nrow=2) @@ -553,11 +551,11 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # Generate a random distribution of betas using the covariance matrix # I think the anomalies might be problematic, so lets get way more betas than we need and trim the distribution if(n.ens==1){ - Rbeta.anom <- matrix(coef(mod.anom), ncol=length(coef(mod.anom))) + Rbeta.anom <- matrix(stats::coef(mod.anom), ncol=length(stats::coef(mod.anom))) } else { - Rbeta.anom <- matrix(MASS::mvrnorm(n=n.new, coef(mod.anom), vcov(mod.anom)), ncol=length(coef(mod.anom))) + Rbeta.anom <- matrix(MASS::mvrnorm(n=n.new, stats::coef(mod.anom), stats::vcov(mod.anom)), ncol=length(stats::coef(mod.anom))) } - dimnames(Rbeta.anom)[[2]] <- names(coef(mod.anom)) + dimnames(Rbeta.anom)[[2]] <- names(stats::coef(mod.anom)) # # Filter our betas to remove outliers # ci.anom <- matrix(apply(Rbeta.anom, 2, quantile, c(0.01, 0.99)), nrow=2) # @@ -575,9 +573,9 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU if(v == "precipitation_flux"){ if(n.ens==1){ - Rbeta.ann <- matrix(coef(mod.ann), ncol=length(coef.ann)) + Rbeta.ann <- matrix(stats::coef(mod.ann), ncol=length(coef.ann)) } else { - Rbeta.ann <- matrix(MASS::mvrnorm(n=n.new, coef(mod.ann), vcov(mod.ann)), ncol=length(coef(mod.ann))) + Rbeta.ann <- matrix(MASS::mvrnorm(n=n.new, stats::coef(mod.ann), stats::vcov(mod.ann)), ncol=length(stats::coef(mod.ann))) } # ci.ann <- matrix(apply(Rbeta.ann, 2, quantile, c(0.01, 0.99)), nrow=2) # Rbeta.ann <- Rbeta.ann[which(apply(Rbeta.ann, 1, function(x) all(x > ci.ann[1,] & x < ci.ann[2,]))),] @@ -585,16 +583,16 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU } # Create the prediction matrix - Xp <- predict(mod.bias, newdata=met.src[met.src$ind==ind,], type="lpmatrix") - Xp.anom <- predict(mod.anom, newdata=met.src[met.src$ind==ind,], type="lpmatrix") + Xp <- stats::predict(mod.bias, newdata=met.src[met.src$ind==ind,], type="lpmatrix") + Xp.anom <- stats::predict(mod.anom, newdata=met.src[met.src$ind==ind,], type="lpmatrix") if(v == "precipitation_flux"){ # Linear models have a bit of a difference in how we get the info out # Xp.ann <- predict(mod.ann, newdata=met.src, type="lpmatrix") met.src[met.src$ind==ind,"Y.tot"] <- met.src[met.src$ind==ind,"pred.ann"] - mod.terms <- terms(mod.ann) - m <- model.frame(mod.terms, met.src[met.src$ind==ind,], xlev=mod.ann$xlevels) - Xp.ann <- model.matrix(mod.terms, m, constrasts.arg <- mod.ann$contrasts) + mod.terms <- stats::terms(mod.ann) + m <- stats::model.frame(mod.terms, met.src[met.src$ind==ind,], xlev=mod.ann$xlevels) + Xp.ann <- stats::model.matrix(mod.terms, m, constrasts.arg <- mod.ann$contrasts) } # ----- @@ -644,7 +642,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # sim1b.norm <- apply(sim1b, 1, mean) # What we need is to remove the mean-trend from the anomalies and then add the trend (with uncertinaties) back in # Note that for a single-member ensemble, this just undoes itself - anom.detrend <- met.src[met.src$ind==ind,"anom.raw"] - predict(mod.anom) + anom.detrend <- met.src[met.src$ind==ind,"anom.raw"] - stats::predict(mod.anom) # NOTE: This section can probably be removed and simplified since it should always be a 1-column array now if(length(cols.redo)>1){ @@ -689,8 +687,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # max air temp = 70 C; hottest temperature from sattellite; very ridiculous # min air temp = -95 C; colder than coldest natural temperature recorded in Antarctica cols.redo <- which(apply(sim1, 2, function(x) min(x) < 273.15-95 | max(x) > 273.15+70 | - min(x) < mean(met.train$X) - sanity.sd*sd(met.train$X) | - max(x) > mean(met.train$X) + sanity.sd*sd(met.train$X) + min(x) < mean(met.train$X) - sanity.sd*stats::sd(met.train$X) | + max(x) > mean(met.train$X) + sanity.sd*stats::sd(met.train$X) )) } #"specific_humidity", @@ -699,8 +697,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # Also, the minimum humidity can't be 0 so lets just make it extremely dry; lets set this for 1 g/Mg cols.redo <- which(apply(sim1, 2, function(x) min(x^2) < 1e-6 | max(x^2) > 40e-3 | - min(x^2) < mean(met.train$X^2) - sanity.sd*sd(met.train$X^2) | - max(x^2) > mean(met.train$X^2) + sanity.sd*sd(met.train$X^2) + min(x^2) < mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2) | + max(x^2) > mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2) )) } #"surface_downwelling_shortwave_flux_in_air", @@ -709,8 +707,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # Lets round 1360 and divide that by 2 (because it should be a daily average) and conservatively assume albedo of 20% (average value is more like 30) # Source http://eesc.columbia.edu/courses/ees/climate/lectures/radiation/ cols.redo <- which(apply(sim1, 2, function(x) max(x^2) > 1360/2*0.8 | - min(x^2) < mean(met.train$X^2) - sanity.sd*sd(met.train$X^2) | - max(x^2) > mean(met.train$X^2) + sanity.sd*sd(met.train$X^2) + min(x^2) < mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2) | + max(x^2) > mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2) )) } if(v == "air_pressure"){ @@ -718,8 +716,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # - Lets round up to 1100 hPA # Also according to Wikipedia, the lowest non-tornadic pressure ever measured was 870 hPA cols.redo <- which(apply(sim1, 2, function(x) min(x) < 870*100 | max(x) > 1100*100 | - min(x) < mean(met.train$X) - sanity.sd*sd(met.train$X) | - max(x) > mean(met.train$X) + sanity.sd*sd(met.train$X) + min(x) < mean(met.train$X) - sanity.sd*stats::sd(met.train$X) | + max(x) > mean(met.train$X) + sanity.sd*stats::sd(met.train$X) )) } if(v == "surface_downwelling_longwave_flux_in_air"){ @@ -728,16 +726,16 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # ED2 sanity checks boudn longwave at 40 & 600 cols.redo <- which(apply(sim1, 2, function(x) min(x^2) < 40 | max(x^2) > 600 | - min(x^2) < mean(met.train$X^2) - sanity.sd*sd(met.train$X^2) | - max(x^2) > mean(met.train$X^2) + sanity.sd*sd(met.train$X^2) + min(x^2) < mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2) | + max(x^2) > mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2) )) } if(v == "wind_speed"){ # According to wikipedia, the hgihest wind speed ever recorded is a gust of 113 m/s; the maximum 5-mind wind speed is 49 m/s cols.redo <- which(apply(sim1, 2, function(x) max(x^2) > 50/2 | - min(x^2) < mean(met.train$X^2) - sanity.sd*sd(met.train$X^2) | - max(x^2) > mean(met.train$X^2) + sanity.sd*sd(met.train$X^2) + min(x^2) < mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2) | + max(x^2) > mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2) )) } if(v == "precipitation_flux"){ @@ -745,8 +743,8 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # https://www.wunderground.com/blog/weatherhistorian/what-is-the-most-rain-to-ever-fall-in-one-minute-or-one-hour.html # 16/2 = round number; x25.4 = inches to mm; /(60*60) = hr to sec cols.redo <- which(apply(sim1, 2, function(x) max(x) > 16/2*25.4/(60*60) | - min(x) < min(met.train$X) - sanity.sd*sd(met.train$X) | - max(x) > max(met.train$X) + sanity.sd*sd(met.train$X) + min(x) < min(met.train$X) - sanity.sd*stats::sd(met.train$X) | + max(x) > max(met.train$X) + sanity.sd*stats::sd(met.train$X) )) } n.new = length(cols.redo) @@ -769,12 +767,12 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # for known problem variables, lets force sanity as a last resort if(v %in% c("air_temperature", "air_temperature_maximum", "air_temperature_minimum")){ warning(paste("Forcing Sanity:", v)) - if(min(sim1) < max(184, mean(met.train$X) - sanity.sd*sd(met.train$X))) { - qtrim <- max(184, mean(met.train$X) - sanity.sd*sd(met.train$X)) + 1e-6 + if(min(sim1) < max(184, mean(met.train$X) - sanity.sd*stats::sd(met.train$X))) { + qtrim <- max(184, mean(met.train$X) - sanity.sd*stats::sd(met.train$X)) + 1e-6 sim1[sim1 < qtrim] <- qtrim } - if(max(sim1) > min(331, mean(met.train$X) + sd(met.train$X^2))) { - qtrim <- min(331, mean(met.train$X) + sanity.sd*sd(met.train$X)) - 1e-6 + if(max(sim1) > min(331, mean(met.train$X) + stats::sd(met.train$X^2))) { + qtrim <- min(331, mean(met.train$X) + sanity.sd*stats::sd(met.train$X)) - 1e-6 sim1[sim1 > qtrim] <- qtrim } } else if(v == "surface_downwelling_shortwave_flux_in_air"){ @@ -786,12 +784,12 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # max(x) > mean(met.train$X) + sanity.sd*sd(met.train$X) # )) warning(paste("Forcing Sanity:", v)) - if(min(sim1^2) < max(mean(met.train$X^2) - sanity.sd*sd(met.train$X^2))) { - qtrim <- max(mean(met.train$X^2) - sanity.sd*sd(met.train$X^2)) + if(min(sim1^2) < max(mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- max(mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 < qtrim] <- sqrt(qtrim) } - if(max(sim1^2) > min(1500*0.8, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2))) { - qtrim <- min(1500*0.8, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2)) + if(max(sim1^2) > min(1500*0.8, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- min(1500*0.8, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 > qtrim] <- sqrt(qtrim) } @@ -800,43 +798,43 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU # ED2 sanity checks boudn longwave at 40 & 600 warning(paste("Forcing Sanity:", v)) - if(min(sim1^2) < max(40, mean(met.train$X^2) - sanity.sd*sd(met.train$X^2))) { - qtrim <- max(40, mean(met.train$X^2) - sanity.sd*sd(met.train$X^2)) + if(min(sim1^2) < max(40, mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- max(40, mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 < qtrim] <- sqrt(qtrim) } - if(max(sim1^2) > min(600, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2))) { - qtrim <- min(600, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2)) + if(max(sim1^2) > min(600, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- min(600, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 > qtrim] <- sqrt(qtrim) } } else if(v=="specific_humidity"){ warning(paste("Forcing Sanity:", v)) # I'm having a hell of a time trying to get SH to fit sanity bounds, so lets brute-force fix outliers - if(min(sim1^2) < max(1e-6, mean(met.train$X^2) - sanity.sd*sd(met.train$X^2))) { - qtrim <- max(1e-6, mean(met.train$X^2) - sanity.sd*sd(met.train$X^2)) + if(min(sim1^2) < max(1e-6, mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- max(1e-6, mean(met.train$X^2) - sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 < qtrim] <- sqrt(qtrim) } - if(max(sim1^2) > min(3.2e-2, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2))) { - qtrim <- min(3.2e-2, mean(met.train$X^2) + sanity.sd*sd(met.train$X^2)) + if(max(sim1^2) > min(3.2e-2, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2))) { + qtrim <- min(3.2e-2, mean(met.train$X^2) + sanity.sd*stats::sd(met.train$X^2)) sim1[sim1^2 > qtrim] <- sqrt(qtrim) } } else if(v=="air_pressure"){ warning(paste("Forcing Sanity:", v)) - if(min(sim1)< max(45000, mean(met.train$X) - sanity.sd*sd(met.train$X))){ - qtrim <- min(45000, mean(met.train$X) - sanity.sd*sd(met.train$X)) + if(min(sim1)< max(45000, mean(met.train$X) - sanity.sd*stats::sd(met.train$X))){ + qtrim <- min(45000, mean(met.train$X) - sanity.sd*stats::sd(met.train$X)) sim1[sim1 < qtrim] <- qtrim } - if(max(sim1) < min(11000000, mean(met.train$X) + sanity.sd*sd(met.train$X))){ - qtrim <- min(11000000, mean(met.train$X) + sanity.sd*sd(met.train$X)) + if(max(sim1) < min(11000000, mean(met.train$X) + sanity.sd*stats::sd(met.train$X))){ + qtrim <- min(11000000, mean(met.train$X) + sanity.sd*stats::sd(met.train$X)) sim1[sim1 > qtrim] <- qtrim } } else if(v=="wind_speed"){ warning(paste("Forcing Sanity:", v)) - if(min(sim1)< max(0, mean(met.train$X) - sanity.sd*sd(met.train$X))){ - qtrim <- min(0, mean(met.train$X) - sanity.sd*sd(met.train$X)) + if(min(sim1)< max(0, mean(met.train$X) - sanity.sd*stats::sd(met.train$X))){ + qtrim <- min(0, mean(met.train$X) - sanity.sd*stats::sd(met.train$X)) sim1[sim1 < qtrim] <- qtrim } - if(max(sim1) < min(sqrt(85), mean(met.train$X) + sanity.sd*sd(met.train$X))){ - qtrim <- min(sqrt(85), mean(met.train$X) + sanity.sd*sd(met.train$X)) + if(max(sim1) < min(sqrt(85), mean(met.train$X) + sanity.sd*stats::sd(met.train$X))){ + qtrim <- min(sqrt(85), mean(met.train$X) + sanity.sd*stats::sd(met.train$X)) sim1[sim1 > qtrim] <- qtrim } } else { @@ -880,14 +878,14 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU } # n.now = number of rainless days for this sim - n.now <- round(rnorm(1, mean(rainless, na.rm=T), sd(rainless, na.rm=T)), 0) + n.now <- round(stats::rnorm(1, mean(rainless, na.rm=T), stats::sd(rainless, na.rm=T)), 0) if(n.now < rainless.min) n.now <- rainless.min # Make sure we don't have negative or no rainless days if(n.now > rainless.max) n.now <- rainless.max # Make sure we have at least one day with rain # We're having major seasonality issues, so lets randomly redistribute our precip # Pull ~twice what we need and randomly select from that so that we don't have such clean cuttoffs # set.seed(12) - cutoff <- quantile(sim1[rows.yr, j], min(n.now/366*2.5, max(0.75, n.now/366)), na.rm=T) + cutoff <- stats::quantile(sim1[rows.yr, j], min(n.now/366*2.5, max(0.75, n.now/366)), na.rm=T) if(length(which(sim1[rows.yr,j]>0)) < n.now){ # if we need to re-distribute our rain (make more rainy days), use the inverse of the cutoff # cutoff <- 1-cutoff @@ -940,7 +938,7 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU } # end z # If we have a worryingly high number of consequtive wet days (outside of 6 sd); try a new dry - if(max(ens.wet) > max(cons.wet)+sd(cons.wet) ){ + if(max(ens.wet) > max(cons.wet)+stats::sd(cons.wet) ){ # print("redistributing dry days") # If we have a wet period that's too long, lets find the random dry that's # closest to the midpoint of the longest @@ -1021,16 +1019,16 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU dat.pred$Date <- as.POSIXct(dat.pred$Date) dat.pred$obs <- apply(source.data[[v]], 1, mean, na.rm=T) dat.pred$mean <- apply(dat.out[[v]], 1, mean, na.rm=T) - dat.pred$lwr <- apply(dat.out[[v]], 1, quantile, 0.025, na.rm=T) - dat.pred$upr <- apply(dat.out[[v]], 1, quantile, 0.975, na.rm=T) + dat.pred$lwr <- apply(dat.out[[v]], 1, stats::quantile, 0.025, na.rm=T) + dat.pred$upr <- apply(dat.out[[v]], 1, stats::quantile, 0.975, na.rm=T) # Plotting the observed and the bias-corrected 95% CI grDevices::png(file.path(path.diagnostics, paste(ens.name, v, "day.png", sep="_")), height=6, width=6, units="in", res=220) print( ggplot2::ggplot(data=dat.pred[dat.pred$Year>=mean(dat.pred$Year)-1 & dat.pred$Year<=mean(dat.pred$Year)+1,]) + - ggplot2::geom_ribbon(ggplot2::aes(x=Date, ymin=lwr, ymax=upr, fill="corrected"), alpha=0.5) + - ggplot2::geom_line(ggplot2::aes(x=Date, y=mean, color="corrected"), size=0.5) + - ggplot2::geom_line(ggplot2::aes(x=Date, y=obs, color="original"), size=0.5) + + ggplot2::geom_ribbon(ggplot2::aes(x=.data$Date, ymin=.data$lwr, ymax=.data$upr, fill="corrected"), alpha=0.5) + + ggplot2::geom_line(ggplot2::aes(x=.data$Date, y=mean, color="corrected"), size=0.5) + + ggplot2::geom_line(ggplot2::aes(x=.data$Date, y=.data$obs, color="original"), size=0.5) + ggplot2::scale_color_manual(values=c("corrected" = "red", "original"="black")) + ggplot2::scale_fill_manual(values=c("corrected" = "red", "original"="black")) + ggplot2::guides(fill=F) + @@ -1053,14 +1051,14 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU grDevices::png(file.path(path.diagnostics, paste(ens.name, v, "day2.png", sep="_")), height=6, width=6, units="in", res=220) print( ggplot2::ggplot(data=stack.sims[stack.sims$Year>=mean(stack.sims$Year)-2 & stack.sims$Year<=mean(stack.sims$Year)+2,]) + - ggplot2::geom_line(ggplot2::aes(x=Date, y=values, color=ind), size=0.2, alpha=0.8) + + ggplot2::geom_line(ggplot2::aes(x=.data$Date, y=values, color=ind), size=0.2, alpha=0.8) + ggplot2::ggtitle(paste0(v, " - example ensemble members (daily slice)")) + ggplot2::theme_bw() ) grDevices::dev.off() # Looking tat the annual means over the whole time series to make sure we're getting decent interannual variability - dat.yr <- aggregate(dat.pred[,c("obs", "mean", "lwr", "upr")], + dat.yr <- stats::aggregate(dat.pred[,c("obs", "mean", "lwr", "upr")], by=list(dat.pred$Year), FUN=mean) names(dat.yr)[1] <- "Year" @@ -1068,9 +1066,9 @@ debias.met.regression <- function(train.data, source.data, n.ens, vars.debias=NU grDevices::png(file.path(path.diagnostics, paste(ens.name, v, "annual.png", sep="_")), height=6, width=6, units="in", res=220) print( ggplot2::ggplot(data=dat.yr[,]) + - ggplot2::geom_ribbon(ggplot2::aes(x=Year, ymin=lwr, ymax=upr, fill="corrected"), alpha=0.5) + - ggplot2::geom_line(ggplot2::aes(x=Year, y=mean, color="corrected"), size=0.5) + - ggplot2::geom_line(ggplot2::aes(x=Year, y=obs, color="original"), size=0.5) + + ggplot2::geom_ribbon(ggplot2::aes(x=.data$Year, ymin=.data$lwr, ymax=.data$upr, fill="corrected"), alpha=0.5) + + ggplot2::geom_line(ggplot2::aes(x=.data$Year, y=mean, color="corrected"), size=0.5) + + ggplot2::geom_line(ggplot2::aes(x=.data$Year, y=.data$obs, color="original"), size=0.5) + ggplot2::scale_color_manual(values=c("corrected" = "red", "original"="black")) + ggplot2::scale_fill_manual(values=c("corrected" = "red", "original"="black")) + ggplot2::guides(fill=F) + diff --git a/modules/data.atmosphere/R/download.Fluxnet2015.R b/modules/data.atmosphere/R/download.Fluxnet2015.R index 80af35b3f2b..fb98328ee30 100644 --- a/modules/data.atmosphere/R/download.Fluxnet2015.R +++ b/modules/data.atmosphere/R/download.Fluxnet2015.R @@ -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) diff --git a/modules/data.atmosphere/R/download.NARR.R b/modules/data.atmosphere/R/download.NARR.R index 071c442ad74..32d98314883 100644 --- a/modules/data.atmosphere/R/download.NARR.R +++ b/modules/data.atmosphere/R/download.NARR.R @@ -8,7 +8,7 @@ ##' @param end_date desired end date YYYY-MM-DD ##' @param ... other inputs ##' example options(download.ftp.method="ncftpget") -##' @importFrom magrittr %>% +##' @importFrom dplyr %>% ##' ##' @examples ##' \dontrun{ diff --git a/modules/data.atmosphere/R/download.NARR_site.R b/modules/data.atmosphere/R/download.NARR_site.R index 703546874e2..9de280c15d7 100644 --- a/modules/data.atmosphere/R/download.NARR_site.R +++ b/modules/data.atmosphere/R/download.NARR_site.R @@ -233,7 +233,7 @@ get_NARR_thredds <- function(start_date, end_date, lat.in, lon.in, get_dfs$data <- foreach::`%dopar%`( foreach::foreach( url = get_dfs$url, flx = get_dfs$flx, - .packages = c("PEcAn.data.atmosphere", "magrittr"), + .packages = c("PEcAn.data.atmosphere", "dplyr"), .export = c("get_narr_url", "robustly") ), PEcAn.utils::robustly(get_narr_url)(url, xy = xy, flx = flx) diff --git a/modules/data.atmosphere/R/download.NEONmet.R b/modules/data.atmosphere/R/download.NEONmet.R index 87cdb8faaaf..f1262701790 100644 --- a/modules/data.atmosphere/R/download.NEONmet.R +++ b/modules/data.atmosphere/R/download.NEONmet.R @@ -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) diff --git a/modules/data.atmosphere/R/download.US_Wlef.R b/modules/data.atmosphere/R/download.US_Wlef.R index 4d44c880f1a..b6389c3f24a 100644 --- a/modules/data.atmosphere/R/download.US_Wlef.R +++ b/modules/data.atmosphere/R/download.US_Wlef.R @@ -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 diff --git a/modules/data.atmosphere/R/download_noaa_gefs_efi.R b/modules/data.atmosphere/R/download_noaa_gefs_efi.R index 0483b8c3331..2837800f508 100644 --- a/modules/data.atmosphere/R/download_noaa_gefs_efi.R +++ b/modules/data.atmosphere/R/download_noaa_gefs_efi.R @@ -40,10 +40,10 @@ download_NOAA_GEFS_EFI <- function(sitename, outfolder, start_date, site.lat, si noaa_data[v] <- NULL #filter for met variable - curr_var <- filter(weather, .data$variable == cf_var_names[v]) + curr_var <- dplyr::filter(weather, .data$variable == cf_var_names[v]) #remove ensemble member 31 does not cover full timeseries #this is a HACK should add a generalized method for ensemble member outlier detection - curr_var <- filter(curr_var, .data$parameter <= 30) + curr_var <- dplyr::filter(curr_var, .data$parameter <= 30) noaa_data[[v]] <- list(value = curr_var$prediction, ensembles = curr_var$parameter, forecast.date = curr_var$datetime) diff --git a/modules/data.atmosphere/R/downscaling_helper_functions.R b/modules/data.atmosphere/R/downscaling_helper_functions.R index 75734f4d6d9..85532c9a186 100644 --- a/modules/data.atmosphere/R/downscaling_helper_functions.R +++ b/modules/data.atmosphere/R/downscaling_helper_functions.R @@ -112,7 +112,7 @@ downscale_repeat_6hr_to_hrly <- function(df, varName, hr = 1){ t0 <- min(df$time) df <- df %>% - dplyr::select("time", all_of(varName)) %>% + dplyr::select("time", tidyselect::all_of(varName)) %>% #Calculate time difference dplyr::mutate(days_since_t0 = difftime(.data$time, t0, units = "days")) %>% #Shift valued back because the 6hr value represents the average over the diff --git a/modules/data.atmosphere/R/extract.nc.module.R b/modules/data.atmosphere/R/extract.nc.module.R index c6fab2ab600..a65969e13a9 100644 --- a/modules/data.atmosphere/R/extract.nc.module.R +++ b/modules/data.atmosphere/R/extract.nc.module.R @@ -1,5 +1,4 @@ ##' @export -##' @import dplyr .extract.nc.module <- function(cf.id, register, dir, met, str_ns, site, new.site, con, start_date, end_date, host, overwrite = FALSE) { PEcAn.logger::logger.info("Site Extraction") diff --git a/modules/data.atmosphere/R/extract_ERA5.R b/modules/data.atmosphere/R/extract_ERA5.R index ddbcb78d4af..4bc5a9f05a0 100644 --- a/modules/data.atmosphere/R/extract_ERA5.R +++ b/modules/data.atmosphere/R/extract_ERA5.R @@ -122,13 +122,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. diff --git a/modules/data.atmosphere/R/half_hour_downscale.R b/modules/data.atmosphere/R/half_hour_downscale.R index d2302ea247e..bb14748412a 100644 --- a/modules/data.atmosphere/R/half_hour_downscale.R +++ b/modules/data.atmosphere/R/half_hour_downscale.R @@ -229,7 +229,7 @@ downscale_ShortWave_to_half_hrly <- function(df,lat, lon, hr = 0.5){ } #ShortWave.ds <- dplyr::select(data.hrly, time, surface_downwelling_shortwave_flux_in_air) - ShortWave.ds <- data.hrly %>% select("time", "surface_downwelling_shortwave_flux_in_air") + ShortWave.ds <- data.hrly %>% dplyr::select("time", "surface_downwelling_shortwave_flux_in_air") # data.hrly$group_6hr <- NA # # group <- 0 @@ -277,7 +277,7 @@ downscale_repeat_6hr_to_half_hrly <- function(df, varName, hr = 0.5){ t0 <- min(df$time) df <- df %>% - dplyr::select("time", all_of(varName)) %>% + dplyr::select("time", tidyselect::all_of(varName)) %>% #Calculate time difference dplyr::mutate(days_since_t0 = difftime(.data$time, t0, units = "days")) %>% #Shift valued back because the 6hr value represents the average over the diff --git a/modules/data.atmosphere/R/met_functions.R b/modules/data.atmosphere/R/lightME.R similarity index 68% rename from modules/data.atmosphere/R/met_functions.R rename to modules/data.atmosphere/R/lightME.R index 3addec2de33..e7e6667bdd4 100644 --- a/modules/data.atmosphere/R/met_functions.R +++ b/modules/data.atmosphere/R/lightME.R @@ -1,34 +1,4 @@ -get.weather <- function(lat, lon, met.nc = met.nc, start.date, end.date, output.dt = 1) { - # if(!is.land(lat, lon)) stop('point is in ocean') - result <- load.cfmet(lat = lat, lon = lon, met.nc = met.nc, start.date, end.date) - downscaled.result <- cfmet.downscale.time(cfmet = result, output.dt = output.dt, lat = lat) - weather <- cruncep_dt2weather(downscaled.result) -} # get.weather - - -is.land <- function(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)) - loni <- which.min(abs(Lon - lon)) - mask <- ncdf4::ncvar_get(nc = met.nc, varid = "mask", start = c(loni, lati), count = c(1, 1)) - return(mask >= 0) -} # is.land - -get.latlonbox <- function(lati, loni, Lat = Lat, Lon = Lon) { - lat <- c(mean(Lat[lati:(lati - 1)]), mean(Lat[lati:(lati + 1)])) - lon <- c(mean(Lon[loni:(loni - 1)]), mean(Lon[loni:(loni + 1)])) - return(c(sort(lat), sort(lon))) -} # get.latlonbox - -get.cruncep <- function(lat, lon, start.date, end.date) { - result <- load.cfmet(lat, lon) - hourly.result <- cruncep_hourly(result, lat = Lat[lati]) - weather <- cruncep_dt2weather(hourly.result) - return(weather) -} # get.cruncep - ##' Simulates the light macro environment ##' ##' Simulates light macro environment based on latitude, day of the year. diff --git a/modules/data.atmosphere/R/load.cfmet.R b/modules/data.atmosphere/R/load.cfmet.R index b7b03532899..f8a9e0cee56 100644 --- a/modules/data.atmosphere/R/load.cfmet.R +++ b/modules/data.atmosphere/R/load.cfmet.R @@ -1,16 +1,16 @@ -## ensures data.table objects treated as such http://stackoverflow.com/q/24501245/513006 -.datatable.aware <- TRUE ##' Load met data from PEcAn formatted met driver ##' -##' subsets a PEcAn formatted met driver file and converts to a data.table / data.frame object -##' @title load CF met +##' subsets a PEcAn formatted met driver file and converts to a data.frame object +##' ##' @param met.nc object of class ncdf4 representing an open CF compliant, PEcAn standard netcdf file with met data -##' @param lat numeric value of latutude +##' @param lat numeric value of latitude ##' @param lon numeric value of longitude ##' @param start.date format is 'YYYY-MM-DD' ##' @param end.date format is 'YYYY-MM-DD' -##' @return data.table of met data +##' @return data frame of met data +##' @importFrom rlang .data +##' @importFrom dplyr %>% ##' @export ##' @author David LeBauer load.cfmet <- function(met.nc, lat, lon, start.date, end.date) { @@ -43,15 +43,13 @@ load.cfmet <- function(met.nc, lat, lon, start.date, end.date) { base.units <- strsplit(basetime.string, " since ")[[1]][1] ## convert to days - if (!base.units == "days") { + if (base.units != "days") { time.idx <- PEcAn.utils::ud_convert(time.idx, basetime.string, paste("days since ", base.date)) } time.idx <- PEcAn.utils::ud_convert(time.idx, "days", "seconds") - date <- as.POSIXct.numeric(time.idx, origin = base.date, tz = "UTC") + date <- as.POSIXct.numeric(round(time.idx), origin = base.date, tz = "UTC") - ## data table warns not to use POSIXlt, which is induced by round() - ## but POSIXlt moves times off by a second - suppressWarnings(all.dates <- data.table::data.table(index = seq(time.idx), date = round(date))) + all.dates <- data.frame(index = seq_along(time.idx), date = date) if (start.date + lubridate::days(1) < min(all.dates$date)) { PEcAn.logger::logger.severe("run start date", start.date, "before met data starts", min(all.dates$date)) @@ -60,15 +58,15 @@ load.cfmet <- function(met.nc, lat, lon, start.date, end.date) { PEcAn.logger::logger.severe("run end date", end.date, "after met data ends", max(all.dates$date)) } - run.dates <- all.dates[date >= start.date & date <= end.date, - list(index, - date = date, - doy = lubridate::yday(date), - year = lubridate::year(date), - month = lubridate::month(date), - day = lubridate::day(date), - hour = lubridate::hour(date) + lubridate::minute(date) / 60)] - + run.dates <- all.dates %>% + dplyr::filter(.data$date >= start.date & .data$date <= end.date) %>% + dplyr::mutate( + doy = lubridate::yday(.data$date), + year = lubridate::year(.data$date), + month = lubridate::month(.data$date), + day = lubridate::day(.data$date), + hour = lubridate::hour(.data$date) + lubridate::minute(.data$date) / 60) + results <- list() @@ -83,5 +81,5 @@ load.cfmet <- function(met.nc, lat, lon, start.date, end.date) { names(vars) <- gsub("surface_pressure", "air_pressure", variables) - return(cbind(run.dates, data.table::as.data.table(vars[!sapply(vars, is.null)]))) + return(cbind(run.dates, vars[!sapply(vars, is.null)])) } # load.cfmet diff --git a/modules/data.atmosphere/R/met.process.R b/modules/data.atmosphere/R/met.process.R index 8eef41130ae..c2e253586d9 100644 --- a/modules/data.atmosphere/R/met.process.R +++ b/modules/data.atmosphere/R/met.process.R @@ -19,6 +19,7 @@ ##' *except* raw met downloads. I.e., it corresponds to: ##' ##' list(download = FALSE, met2cf = TRUE, standardize = TRUE, met2model = TRUE) +##' @importFrom rlang .data .env ##' @export ##' @author Elizabeth Cowdery, Michael Dietze, Ankur Desai, James Simkins, Ryan Kelly met.process <- function(site, input_met, start_date, end_date, model, @@ -227,25 +228,25 @@ 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) - # Fidning the tiles. - raw.tiles <- tbl(con, "inputs") %>% - filter( - site_id == register$ParentSite, - start_date >= start_date, - end_date <= end_date, - format_id == formatid + machine.id <- PEcAn.DB::get.id(table = "machines", "hostname", PEcAn.remote::fqdn(), con) + # Finding the tiles. + raw.tiles <- dplyr::tbl(con, "inputs") %>% + dplyr::filter( + .data$site_id == register$ParentSite, + .data$start_date <= .env$start_date, + .data$end_date >= .env$end_date, + .data$format_id == formatid ) %>% - filter(grepl(met, "name")) %>% - inner_join(tbl(con, "dbfiles"), by = c('id' = 'container_id')) %>% - filter(machine_id == machine.id) %>% - collect() + dplyr::filter(grepl(met, "name")) %>% + dplyr::inner_join(dplyr::tbl(con, "dbfiles"), by = c('id' = 'container_id')) %>% + dplyr::filter(.data$machine_id == machine.id) %>% + dplyr::collect() cf.id <- raw.id <- list(input.id = raw.tiles$id.x, dbfile.id = raw.tiles$id.y) } diff --git a/modules/data.atmosphere/R/met_temporal_downscale.Gaussian_ensemble.R b/modules/data.atmosphere/R/met_temporal_downscale.Gaussian_ensemble.R index 6bb62c1bfe9..375c44e1857 100644 --- a/modules/data.atmosphere/R/met_temporal_downscale.Gaussian_ensemble.R +++ b/modules/data.atmosphere/R/met_temporal_downscale.Gaussian_ensemble.R @@ -169,7 +169,10 @@ 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 @@ -196,7 +199,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) } diff --git a/modules/data.atmosphere/R/metgapfill.NOAA_GEFS.R b/modules/data.atmosphere/R/metgapfill.NOAA_GEFS.R index 4c41614bd1d..1e573a65706 100644 --- a/modules/data.atmosphere/R/metgapfill.NOAA_GEFS.R +++ b/modules/data.atmosphere/R/metgapfill.NOAA_GEFS.R @@ -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) } } } @@ -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)) { diff --git a/modules/data.atmosphere/R/metutils.R b/modules/data.atmosphere/R/metutils.R index 4ef16906842..c5662e2db66 100644 --- a/modules/data.atmosphere/R/metutils.R +++ b/modules/data.atmosphere/R/metutils.R @@ -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 diff --git a/modules/data.atmosphere/R/tdm_lm_ensemble_sims.R b/modules/data.atmosphere/R/tdm_lm_ensemble_sims.R index b9ee80e71f6..d594666611c 100644 --- a/modules/data.atmosphere/R/tdm_lm_ensemble_sims.R +++ b/modules/data.atmosphere/R/tdm_lm_ensemble_sims.R @@ -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 & @@ -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, envir = env) + mod.save <- env$mod.save + } # Pull coefficients (betas) from our saved matrix @@ -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) } } diff --git a/modules/data.atmosphere/R/tdm_model_train.R b/modules/data.atmosphere/R/tdm_model_train.R index 685f7880590..0cff73e411d 100644 --- a/modules/data.atmosphere/R/tdm_model_train.R +++ b/modules/data.atmosphere/R/tdm_model_train.R @@ -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 - @@ -41,7 +41,7 @@ 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% @@ -49,7 +49,7 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL } 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 - @@ -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, @@ -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) @@ -120,15 +120,15 @@ 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, ]) } @@ -136,47 +136,47 @@ model.train <- function(dat.subset, v, n.beta, resids = resids, threshold = NULL 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)) } diff --git a/modules/data.atmosphere/R/tdm_subdaily_pred.R b/modules/data.atmosphere/R/tdm_subdaily_pred.R index 1b02852519c..a7bca5cb515 100644 --- a/modules/data.atmosphere/R/tdm_subdaily_pred.R +++ b/modules/data.atmosphere/R/tdm_subdaily_pred.R @@ -56,7 +56,7 @@ subdaily_pred <- function(newdata, model.predict, Rbeta, resid.err = FALSE, mode # newdata <- newdata[order(newdata$ens, newdata$hour),] } - Xp <- model.matrix(eval(model.predict$formula), m, contrasts.arg=model.predict$contr) + Xp <- stats::model.matrix(eval(model.predict$formula), m, contrasts.arg=model.predict$contr) if (resid.err == TRUE) { newdata$resid <- 99999 @@ -66,7 +66,7 @@ subdaily_pred <- function(newdata, model.predict, Rbeta, resid.err = FALSE, mode resid.m <- newdata[,model.resid$factors] resid.m[,"as.ordered(hour)"] <- resid.m$hour if(length(df.hr$hour)!= length(resid.m$hour)) resid.m <- merge(resid.m, df.hr, all=T) - Xp.res <- model.matrix(eval(model.resid$formula), resid.m, contrasts.arg=model.resid$contr) + Xp.res <- stats::model.matrix(eval(model.resid$formula), resid.m, contrasts.arg=model.resid$contr) err.resid <- Xp.res[, resid.piv] %*% t(Rbeta.resid) } # End residual error diff --git a/modules/data.atmosphere/R/temporal.downscaling.R b/modules/data.atmosphere/R/temporal.downscaling.R index 4032fe61135..21d10026dbb 100644 --- a/modules/data.atmosphere/R/temporal.downscaling.R +++ b/modules/data.atmosphere/R/temporal.downscaling.R @@ -1,15 +1,15 @@ -## ensures data.table objects treated as such http://stackoverflow.com/q/24501245/513006 -.datatable.aware <- TRUE -##' Temporal downscaling of daily or subdaily met data -##' -##' @title Downscale CF met data -##' @param cfmet data.table with CF variables generated by \code{\link{load.cfmet}} -##' @param output.dt time step (hours) for output -##' @return downscaled result -##' @export -##' @author David LeBauer -cfmet.downscale.time <- cruncep_hourly <- function(cfmet, output.dt = 1, lat = lat, ...) { +#' Temporal downscaling of daily or subdaily CF met data +#' +#' @param cfmet data frame with CF variables generated by \code{\link{load.cfmet}} +#' @param output.dt time step (hours) for output +#' @param lat latitude (for calculating solar radiation) +#' @param ... ignored +#' +#' @return downscaled result +#' @export +#' @author David LeBauer +cfmet.downscale.time <- function(cfmet, output.dt = 1, lat = lat, ...) { ## time step dt_hr <- as.numeric(round(difftime(cfmet$date[2], cfmet$date[1], units = "hours"))) @@ -28,13 +28,16 @@ cfmet.downscale.time <- cruncep_hourly <- function(cfmet, output.dt = 1, lat = l } else if (dt_hr > 6 & dt_hr < 24) { # cfmet <- cfmet[,list(air_temperature_max = max(air_temperature), air_temperature_min = # min(air_temperature), ), by = 'year,doy']) dt_hr <- 24 - PEcAn.logger::logger.error("timestep of input met data is between 6 and 24 hours.\n", "PEcAn will automatically convert this to daily data\n", + PEcAn.logger::logger.error("timestep of input met data is between 6 and 24 hours.\n", "PEcAn will automatically convert this to daily data\n", "you should confirm validity of downscaling, in particular that min / max temperatures are realistic") } if (dt_hr == 24) { if (all(c("tmax", "tmin") %in% colnames(cfmet))) { - data.table::setnames(cfmet, c("tmax", "tmin"), c("air_temperature_max", "air_temperature_min")) + nm <- colnames(cfmet) + nm[nm == "tmax"] <- "air_temperature_max" + nm[nm == "tmin"] <- "air_temperature_min" + colnames(cfmet) <- nm } downscaled.result <- cfmet.downscale.daily(dailymet = cfmet, output.dt = output.dt, lat = lat) } else if (dt_hr > 24) { @@ -49,18 +52,19 @@ cfmet.downscale.time <- cruncep_hourly <- function(cfmet, output.dt = 1, lat = l ##' ##' Uses simple spline to interpolate variables with diurnal variability, otherwise uses averaging or repeating ##' for variables with no clear diurnal pattern. For all variables except temperature, negative values are set to zero. -##' @title subdaily downscaling -##' @param subdailymet data table with climate variables queried from \code{\link{load.cfmet}} +##' +##' @param subdailymet data frame with climate variables queried from \code{\link{load.cfmet}} ##' @param output.dt output timestep. default is one hour ##' @export ##' @return weather file with subdaily met variables rescaled to output time step ##' @author David LeBauer cfmet.downscale.subdaily <- function(subdailymet, output.dt = 1) { ## converting surface_downwelling_shortwave_flux_in_air from W/m2 avg to PPFD - new.date <- subdailymet[,list(hour = 0:(23 / output.dt) / output.dt), - by = c("year", "month", "day", "doy")] + new.date <- subdailymet %>% + dplyr::group_by(.data$year, .data$month, .data$day, .data$doy) %>% + dplyr::group_modify(~data.frame(hour = 0:(23 / output.dt) / output.dt)) - new.date$date <- new.date[,list(date = lubridate::ymd_h(paste(year, month, day, hour)))] + new.date$date <- lubridate::ymd_h(paste(new.date$year, new.date$month, new.date$day, new.date$hour)) downscaled.result <- list() tint <- nrow(new.date)/ nrow(subdailymet) @@ -88,7 +92,7 @@ cfmet.downscale.subdaily <- function(subdailymet, output.dt = 1) { "surface_downwelling_shortwave_flux_in_air", "ppfd", "relative_humidity")){ if(var %in% colnames(subdailymet)){ ## convert units from subdaily to hourly - hrscale <- ifelse(var %in% c("surface_downwelling_shortwave_flux_in_air", "precipitation_flux"), + hrscale <- ifelse(var %in% c("surface_downwelling_shortwave_flux_in_air", "precipitation_flux"), output.dt, 1) f <- stats::splinefun(as.double(subdailymet$date), (subdailymet[[var]] / hrscale), method = "monoH.FC") @@ -100,85 +104,116 @@ cfmet.downscale.subdaily <- function(subdailymet, output.dt = 1) { } } - downscaled.result <- cbind(new.date, data.table::as.data.table(downscaled.result)) + downscaled.result <- cbind(new.date, as.data.frame(downscaled.result)) } # cfmet.downscale.subdaily +#' Internal helper to downscale a single row from a daily file +#' +#' @param df one row from dailymet +#' @param tseq vector of hours at which to estimate +#' @param lat latitude +#' +#' @return df with one row for each hour in `tseq` +#' +downscale_one_cfmet_day <- function(df, tseq, lat) { + if (nrow(df) != 1) { + PEcAn.logger::logger.error("df must be a one-row dataframe") + } + if (length(unique(diff(tseq))) != 1) { + PEcAn.logger::logger.error("tseq has uneven intervals") + } + + n <- length(tseq) + + light <- lightME(DOY = df$doy, t.d = tseq, lat = lat) %>% + as.data.frame() %>% + dplyr::mutate( + Itot = .data$I.dir + .data$I.diff, + resC2 = (.data$Itot - min(.data$Itot)) / max(.data$Itot)) + + rhscale <- (cos(2 * pi * (tseq - 10) / n) + 1) / 2 + + data.frame( + hour = tseq, + solarR = rep( + df$surface_downwelling_shortwave_flux_in_air * 2.07 * 10^5 / 36000, + each = n) + * light$resC2, + # Note: When dt >= 6, all times get the same T prediction + Temp = df$tmin + + (sin(2 * pi * (tseq - 10) / n) + 1) / + 2 * (df$tmax - df$tmin), + relative_humidity = df$rhmin + rhscale * (df$rhmax - df$rhmin), + # TODO: Why do we divide by n? + # isn't precipitation_flux already an intensity? + precipitation_flux = rep(df$precipitation_flux / n, each = n), + wind = rep(df$wind_speed, each = n)) %>% + dplyr::mutate( + # TODO: Computation of solarR above already multiplies by resC2. + # Is multiplying it again here really correct? + # That's how the old data.table version did it + # (once when computing `solarR` and again when computing `SolarR`), + # so keeping it until proven wrong. + downwelling_photosynthetic_photon_flux = .data$solarR * light$resC2) +} + + ##' Simple, Fast Daily to Hourly Climate Downscaling ##' ##' Based on weach family of functions but 5x faster than weachNEW, -##' and requiring metric units (temperature in celsius, windspeed in kph, -##' precip in mm, relative humidity as fraction). +##' and requiring metric units (temperature in Kelvins on input and celsius on +##' output, windspeed in kph, precip in mm, relative humidity as fraction). ##' Derived from the weachDT function in the BioCro package. -##' @title daily to subdaily downscaling -##' @param dailymet data table with climate variables -##' @param lat latitude (for calculating solar radiation) +##' +##' @param dailymet data frame with climate variables ##' @param output.dt output timestep +##' @param lat latitude (for calculating solar radiation) +#' +#' @importFrom dplyr %>% +#' @importFrom rlang .data ##' @export ##' @return weather file with subdaily timesteps ##' @author David LeBauer cfmet.downscale.daily <- function(dailymet, output.dt = 1, lat) { - + tint <- 24/output.dt tseq <- seq(from = 0, to = 23, by = output.dt) - - data.table::setkeyv(dailymet, c("year", "doy")) - + if (all(c("air_temperature_max", "air_temperature_min") %in% colnames(dailymet))) { - data.table::setnames(dailymet, c("air_temperature_max", "air_temperature_min"), c("tmax", "tmin")) + nm <- colnames(dailymet) + nm[nm == "air_temperature_max"] <- "tmax" + nm[nm == "air_temperature_min"] <- "tmin" + colnames(dailymet) <- nm } - - light <- dailymet[, lightME(DOY = doy, t.d = tseq, lat = lat), by = c("year", "doy")] - - light$Itot <- light[, list(I.dir + I.diff)] - resC2 <- light[, list(resC2 = (Itot - min(Itot))/max(Itot)), by = c("year", "doy")]$resC2 - solarR <- dailymet[, list(year, doy, solarR = rep(surface_downwelling_shortwave_flux_in_air * - 2.07 * 10^5/36000, each = tint) * resC2)] - - SolarR <- cbind(resC2, solarR)[, list(SolarR = solarR * resC2)]$SolarR - - ## Temperature - Temp <- dailymet[, list(Temp = tmin + (sin(2 * pi * (tseq - 10)/tint) + 1)/2 * (tmax - tmin), - hour = tseq), by = "year,doy"]$Temp - - ## Relative Humidity - RH <- dailymet[, list(RH = rep(relative_humidity, each = tint), hour = tseq), by = "year,doy"] - data.table::setkeyv(RH, c("year", "doy", "hour")) - - # if(!'air_pressure' %in% colnames(dailymet)) air_pressure <- - qair <- dailymet[, list(year, doy, tmin, tmax, air_pressure, air_temperature, qmin = rh2qair(rh = relative_humidity/100, - T = tmin), qmax = rh2qair(rh = relative_humidity/100, T = tmax))] - - a <- qair[, list(year, doy, tmin, tmax, air_temperature, qmin, qmax, pressure = PEcAn.utils::ud_convert(air_pressure, - "Pa", "millibar"))][, list(year, doy, rhmin = qair2rh(qmin, air_temperature, pressure), rhmax = qair2rh(qmax, - air_temperature, pressure))] - rhscale <- (cos(2 * pi * (tseq - 10)/tint) + 1)/2 - RH <- a[, list(RH = rhmin + rhscale * (rhmax - rhmin)), by = c("year", "doy")]$RH - - ## Wind Speed - - - if ("wind_speed" %in% colnames(dailymet)) { - wind_speed <- rep(dailymet$wind_speed, each = tint) - } else if (all(c("eastward_wind", "northward_wind") %in% colnames(dailymet))) { - northward_wind <- rep(dailymet$northward_wind, each = tint) - eastward_wind <- rep(dailymet$eastward_wind, each = tint) - if (!"wind_speed" %in% colnames(dailymet)) { - wind_speed <- sqrt(northward_wind^2 + eastward_wind^2) + + if (! "wind_speed" %in% colnames(dailymet)) { + if (all(c("eastward_wind", "northward_wind") %in% colnames(dailymet))) { + dailymet$wind_speed <- sqrt(dailymet$northward_wind^2 + dailymet$eastward_wind^2) + } else { + PEcAn.logger::logger.error("no wind_speed found in daily met dataset") } - } else { - PEcAn.logger::logger.error("no wind_speed found in daily met dataset") } - - ## Precipitation - precip <- rep(dailymet$precipitation_flux/tint, each = tint) - - ## Hour - time <- dailymet[, list(hour = tseq), by = c("year", "doy")] - - ans <- data.table::data.table(time, downwelling_photosynthetic_photon_flux = SolarR, air_temperature = PEcAn.utils::ud_convert(Temp, - "kelvin", "celsius"), relative_humidity = RH, wind = wind_speed, precipitation_flux = precip) - return(ans) + + return(dailymet %>% + dplyr::mutate( + qmin = rh2qair(rh = .data$relative_humidity / 100, T = .data$tmin), + qmax = rh2qair(rh = .data$relative_humidity / 100, T = .data$tmax), + pressure = PEcAn.utils::ud_convert(.data$air_pressure, "Pa", "millibar"), + rhmin = qair2rh(.data$qmin, .data$air_temperature, .data$pressure), + rhmax = qair2rh(.data$qmax, .data$air_temperature, .data$pressure)) %>% + dplyr::group_by(.data$year, .data$doy) %>% + dplyr::group_modify(~downscale_one_cfmet_day(.x, tseq, lat), .keep = TRUE) %>% + dplyr::ungroup() %>% + dplyr::mutate( + air_temperature = PEcAn.utils::ud_convert(.data$Temp, "kelvin", "celsius")) %>% + dplyr::select( + "year", "doy", "hour", + "downwelling_photosynthetic_photon_flux", + "air_temperature", + "relative_humidity", + "wind", + "precipitation_flux")) } # cfmet.downscale.daily @@ -188,16 +223,16 @@ cfmet.downscale.daily <- function(dailymet, output.dt = 1, lat) { ##' streamlining extraction of data from netCDF files ##' with CF-compliant variable names ##' -##' @title Get time series vector from netCDF file ##' @param var name of variable to extract ##' @param lati,loni latitude and longitude to extract -##' @param run.dates data.table of dates to read +##' @param run.dates data frame of dates to read ##' @param met.nc netcdf file with CF variable names +##' ##' @return numeric vector ##' @export ##' @author David Shaner LeBauer get.ncvector <- function(var, lati = lati, loni = loni, run.dates = run.dates, met.nc) { - + start.idx <- c(latitude = lati, longitude = loni, time = run.dates$index[1]) count.idx <- c(latitude = 1, longitude = 1, time = nrow(run.dates)) dim.order <- sapply(met.nc$var$air_temperature$dim, function(x) x$name) @@ -205,7 +240,7 @@ get.ncvector <- function(var, lati = lati, loni = loni, run.dates = run.dates, m ans <- ncdf4::ncvar_get(nc = met.nc, varid = var, start = start.idx[dim.order], count = count.idx[dim.order]) return(as.numeric(ans)) } # ncvar_get2 - + if (var %in% attributes(met.nc$var)$names) { ans <- ncvar_get2(var) } else if (var == "air_pressure") { diff --git a/modules/data.atmosphere/man/cfmet.downscale.daily.Rd b/modules/data.atmosphere/man/cfmet.downscale.daily.Rd index 36f9e387c55..ce37b35b8b5 100644 --- a/modules/data.atmosphere/man/cfmet.downscale.daily.Rd +++ b/modules/data.atmosphere/man/cfmet.downscale.daily.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/temporal.downscaling.R \name{cfmet.downscale.daily} \alias{cfmet.downscale.daily} -\title{daily to subdaily downscaling} +\title{Simple, Fast Daily to Hourly Climate Downscaling} \usage{ cfmet.downscale.daily(dailymet, output.dt = 1, lat) } \arguments{ -\item{dailymet}{data table with climate variables} +\item{dailymet}{data frame with climate variables} \item{output.dt}{output timestep} @@ -17,12 +17,9 @@ cfmet.downscale.daily(dailymet, output.dt = 1, lat) weather file with subdaily timesteps } \description{ -Simple, Fast Daily to Hourly Climate Downscaling -} -\details{ Based on weach family of functions but 5x faster than weachNEW, -and requiring metric units (temperature in celsius, windspeed in kph, -precip in mm, relative humidity as fraction). +and requiring metric units (temperature in Kelvins on input and celsius on + output, windspeed in kph, precip in mm, relative humidity as fraction). Derived from the weachDT function in the BioCro package. } \author{ diff --git a/modules/data.atmosphere/man/cfmet.downscale.subdaily.Rd b/modules/data.atmosphere/man/cfmet.downscale.subdaily.Rd index 018699d8e17..4401bfe4ac0 100644 --- a/modules/data.atmosphere/man/cfmet.downscale.subdaily.Rd +++ b/modules/data.atmosphere/man/cfmet.downscale.subdaily.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/temporal.downscaling.R \name{cfmet.downscale.subdaily} \alias{cfmet.downscale.subdaily} -\title{subdaily downscaling} +\title{Subdaily to hourly (or less) downscaling} \usage{ cfmet.downscale.subdaily(subdailymet, output.dt = 1) } \arguments{ -\item{subdailymet}{data table with climate variables queried from \code{\link{load.cfmet}}} +\item{subdailymet}{data frame with climate variables queried from \code{\link{load.cfmet}}} \item{output.dt}{output timestep. default is one hour} } @@ -15,9 +15,6 @@ cfmet.downscale.subdaily(subdailymet, output.dt = 1) weather file with subdaily met variables rescaled to output time step } \description{ -Subdaily to hourly (or less) downscaling -} -\details{ Uses simple spline to interpolate variables with diurnal variability, otherwise uses averaging or repeating for variables with no clear diurnal pattern. For all variables except temperature, negative values are set to zero. } diff --git a/modules/data.atmosphere/man/cfmet.downscale.time.Rd b/modules/data.atmosphere/man/cfmet.downscale.time.Rd index 5e574f2899b..81c9edbaf7d 100644 --- a/modules/data.atmosphere/man/cfmet.downscale.time.Rd +++ b/modules/data.atmosphere/man/cfmet.downscale.time.Rd @@ -2,20 +2,24 @@ % Please edit documentation in R/temporal.downscaling.R \name{cfmet.downscale.time} \alias{cfmet.downscale.time} -\title{Downscale CF met data} +\title{Temporal downscaling of daily or subdaily CF met data} \usage{ cfmet.downscale.time(cfmet, output.dt = 1, lat = lat, ...) } \arguments{ -\item{cfmet}{data.table with CF variables generated by \code{\link{load.cfmet}}} +\item{cfmet}{data frame with CF variables generated by \code{\link{load.cfmet}}} \item{output.dt}{time step (hours) for output} + +\item{lat}{latitude (for calculating solar radiation)} + +\item{...}{ignored} } \value{ downscaled result } \description{ -Temporal downscaling of daily or subdaily met data +Temporal downscaling of daily or subdaily CF met data } \author{ David LeBauer diff --git a/modules/data.atmosphere/man/download.Fluxnet2015.Rd b/modules/data.atmosphere/man/download.Fluxnet2015.Rd index 3b249de1c7b..80f1e105dfc 100644 --- a/modules/data.atmosphere/man/download.Fluxnet2015.Rd +++ b/modules/data.atmosphere/man/download.Fluxnet2015.Rd @@ -17,7 +17,7 @@ download.Fluxnet2015( } \arguments{ \item{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}} \item{outfolder}{location on disk where outputs will be stored} diff --git a/modules/data.atmosphere/man/download.NEONmet.Rd b/modules/data.atmosphere/man/download.NEONmet.Rd index 7f6d63d153d..bf9708a9b7c 100644 --- a/modules/data.atmosphere/man/download.NEONmet.Rd +++ b/modules/data.atmosphere/man/download.NEONmet.Rd @@ -16,7 +16,7 @@ download.NEONmet( } \arguments{ \item{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}} \item{outfolder}{location on disk where outputs will be stored} diff --git a/modules/data.atmosphere/man/downscale_one_cfmet_day.Rd b/modules/data.atmosphere/man/downscale_one_cfmet_day.Rd new file mode 100644 index 00000000000..99798bdb53e --- /dev/null +++ b/modules/data.atmosphere/man/downscale_one_cfmet_day.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/temporal.downscaling.R +\name{downscale_one_cfmet_day} +\alias{downscale_one_cfmet_day} +\title{Internal helper to downscale a single row from a daily file} +\usage{ +downscale_one_cfmet_day(df, tseq, lat) +} +\arguments{ +\item{df}{one row from dailymet} + +\item{tseq}{vector of hours at which to estimate} + +\item{lat}{latitude} +} +\value{ +df with one row for each hour in `tseq` +} +\description{ +Internal helper to downscale a single row from a daily file +} diff --git a/modules/data.atmosphere/man/get.ncvector.Rd b/modules/data.atmosphere/man/get.ncvector.Rd index 0005102cc60..22614d0f89b 100644 --- a/modules/data.atmosphere/man/get.ncvector.Rd +++ b/modules/data.atmosphere/man/get.ncvector.Rd @@ -11,7 +11,7 @@ get.ncvector(var, lati = lati, loni = loni, run.dates = run.dates, met.nc) \item{lati, loni}{latitude and longitude to extract} -\item{run.dates}{data.table of dates to read} +\item{run.dates}{data frame of dates to read} \item{met.nc}{netcdf file with CF variable names} } @@ -19,9 +19,6 @@ get.ncvector(var, lati = lati, loni = loni, run.dates = run.dates, met.nc) numeric vector } \description{ -Get time series vector from netCDF file -} -\details{ internal convenience function for streamlining extraction of data from netCDF files with CF-compliant variable names diff --git a/modules/data.atmosphere/man/lightME.Rd b/modules/data.atmosphere/man/lightME.Rd index 8564a82d8ec..8bddbd3eeab 100644 --- a/modules/data.atmosphere/man/lightME.Rd +++ b/modules/data.atmosphere/man/lightME.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/met_functions.R +% Please edit documentation in R/lightME.R \name{lightME} \alias{lightME} \title{Simulates the light macro environment} diff --git a/modules/data.atmosphere/man/load.cfmet.Rd b/modules/data.atmosphere/man/load.cfmet.Rd index 04d4fc72e3a..f58e26beda5 100644 --- a/modules/data.atmosphere/man/load.cfmet.Rd +++ b/modules/data.atmosphere/man/load.cfmet.Rd @@ -2,14 +2,14 @@ % Please edit documentation in R/load.cfmet.R \name{load.cfmet} \alias{load.cfmet} -\title{load CF met} +\title{Load met data from PEcAn formatted met driver} \usage{ load.cfmet(met.nc, lat, lon, start.date, end.date) } \arguments{ \item{met.nc}{object of class ncdf4 representing an open CF compliant, PEcAn standard netcdf file with met data} -\item{lat}{numeric value of latutude} +\item{lat}{numeric value of latitude} \item{lon}{numeric value of longitude} @@ -18,13 +18,10 @@ load.cfmet(met.nc, lat, lon, start.date, end.date) \item{end.date}{format is 'YYYY-MM-DD'} } \value{ -data.table of met data +data frame of met data } \description{ -Load met data from PEcAn formatted met driver -} -\details{ -subsets a PEcAn formatted met driver file and converts to a data.table / data.frame object +subsets a PEcAn formatted met driver file and converts to a data.frame object } \author{ David LeBauer diff --git a/modules/data.atmosphere/man/qair2rh.Rd b/modules/data.atmosphere/man/qair2rh.Rd index 6fc85af4d85..c98cea9714e 100644 --- a/modules/data.atmosphere/man/qair2rh.Rd +++ b/modules/data.atmosphere/man/qair2rh.Rd @@ -23,7 +23,7 @@ Convert specific humidity to relative humidity 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} } \author{ David LeBauer diff --git a/modules/data.atmosphere/tests/Rcheck_reference.log b/modules/data.atmosphere/tests/Rcheck_reference.log index be30baa7ac2..855436fb1b1 100644 --- a/modules/data.atmosphere/tests/Rcheck_reference.log +++ b/modules/data.atmosphere/tests/Rcheck_reference.log @@ -1,88 +1,15 @@ -* using log directory ‘/home/tanishq010/pecan/modules/PEcAn.data.atmosphere.Rcheck’ -* using R version 4.2.1 (2022-06-23) +* using log directory ‘/tmp/Rtmp56JiPq/PEcAn.data.atmosphere.Rcheck’ +* using R version 4.1.3 (2022-03-10) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 -* using options ‘--no-tests --no-manual --as-cran’ +* using options ‘--no-manual --as-cran’ * checking for file ‘PEcAn.data.atmosphere/DESCRIPTION’ ... OK * checking extension type ... Package * this is package ‘PEcAn.data.atmosphere’ version ‘1.7.2’ * package encoding: UTF-8 -* checking CRAN incoming feasibility ... WARNING -Maintainer: ‘David LeBauer ’ - -New submission - -License components with restrictions and base license permitting such: - BSD_3_clause + file LICENSE -File 'LICENSE': - University of Illinois/NCSA Open Source License - - Copyright (c) 2012, University of Illinois, NCSA. All rights reserved. - - Permission is hereby granted, free of charge, to any person obtaining - a copy of this software and associated documentation files (the - "Software"), to deal with the Software without restriction, including - without limitation the rights to use, copy, modify, merge, publish, - distribute, sublicense, and/or sell copies of the Software, and to - permit persons to whom the Software is furnished to do so, subject to - the following conditions: - - - Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimers. - - Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimers in the - documentation and/or other materials provided with the distribution. - - Neither the names of University of Illinois, NCSA, nor the names - of its contributors may be used to endorse or promote products - derived from this Software without specific prior written permission. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR - ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE SOFTWARE. - -Unknown, possibly misspelled, fields in DESCRIPTION: - ‘Remotes’ - -Strong dependencies not in mainstream repositories: - nneo, PEcAn.DB, PEcAn.logger, PEcAn.remote, PEcAn.utils - -Found the following (possibly) invalid URLs: - URL: http://fluxnet.fluxdata.org//sites/site-list-and-pages/ (moved to https://fluxnet.org/sites/site-list-and-pages/) - From: man/download.Fluxnet2015.Rd - Status: 200 - Message: OK - URL: http://www.eol.ucar.edu/projects/ceop/dm/documents/refdata_report/eqns.html (moved to https://archive.eol.ucar.edu/projects/ceop/dm/documents/refdata_report/eqns.html) - From: man/qair2rh.Rd - Status: 200 - Message: OK - URL: http://www.fluxdata.org/DataInfo/Dataset%20Doc%20Lib/SynthDataSummary.aspx - From: man/download.FluxnetLaThuile.Rd - Status: Error - Message: libcurl error code 28: - Connection timed out after 60000 milliseconds - URL: http://www.neonscience.org/science-design/field-sites/list (moved to https://www.neonscience.org/science-design/field-sites/list) - From: man/download.NEONmet.Rd - Status: 404 - Message: Not Found - URL: https://pecanproject.github.io/pecan//modules/data.atmosphere/inst/web/index.html - From: README.md - Status: 404 - Message: Not Found - -The Date field is over a month old. * checking package namespace information ... OK -* checking package dependencies ... WARNING -Imports includes 40 non-default packages. -Importing from so many packages makes the package vulnerable to any of -them becoming unavailable. Move as many as possible to Suggests and -use conditionally. - * checking package dependencies ... NOTE -Imports includes 40 non-default packages. +Imports includes 38 non-default packages. Importing from so many packages makes the package vulnerable to any of them becoming unavailable. Move as many as possible to Suggests and use conditionally. @@ -115,136 +42,24 @@ Use \uxxxx escapes for other characters. * checking whether the namespace can be loaded with stated dependencies ... OK * checking whether the namespace can be unloaded cleanly ... OK * checking loading without being on the library search path ... OK -* checking use of S3 registration ... OK -* checking dependencies in R code ... WARNING -'library' or 'require' calls not declared from: - ‘MASS’ ‘mgcv’ -'library' or 'require' calls in package code: - ‘MASS’ ‘mgcv’ - Please use :: or requireNamespace() instead. - See section 'Suggested packages' in the 'Writing R Extensions' manual. -Package in Depends field not imported from: ‘methods’ - These packages need to be imported from (in the NAMESPACE file) - for when this namespace is loaded but not attached. +* checking dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... NOTE -cfmet.downscale.daily: no visible binding for global variable ‘doy’ -cfmet.downscale.daily: no visible binding for global variable ‘I.dir’ -cfmet.downscale.daily: no visible binding for global variable ‘I.diff’ -cfmet.downscale.daily: no visible binding for global variable ‘Itot’ -cfmet.downscale.daily: no visible binding for global variable ‘year’ -cfmet.downscale.daily: no visible binding for global variable - ‘surface_downwelling_shortwave_flux_in_air’ -cfmet.downscale.daily: no visible binding for global variable ‘tmin’ -cfmet.downscale.daily: no visible binding for global variable ‘tmax’ -cfmet.downscale.daily: no visible binding for global variable - ‘relative_humidity’ -cfmet.downscale.daily: no visible binding for global variable - ‘air_pressure’ -cfmet.downscale.daily: no visible binding for global variable - ‘air_temperature’ -cfmet.downscale.daily: no visible binding for global variable ‘qmin’ -cfmet.downscale.daily: no visible binding for global variable ‘qmax’ -cfmet.downscale.daily: no visible binding for global variable - ‘pressure’ -cfmet.downscale.daily: no visible binding for global variable ‘rhmin’ -cfmet.downscale.daily: no visible binding for global variable ‘rhmax’ -cfmet.downscale.subdaily: no visible binding for global variable ‘year’ -cfmet.downscale.subdaily: no visible binding for global variable - ‘month’ -cfmet.downscale.subdaily: no visible binding for global variable ‘day’ -cfmet.downscale.subdaily: no visible binding for global variable ‘hour’ -debias.met.regression: no visible global function definition for ‘sd’ -debias.met.regression: no visible global function definition for - ‘aggregate’ -debias.met.regression: no visible global function definition for - ‘predict’ -debias.met.regression: no visible global function definition for - ‘resid’ -debias.met.regression: no visible global function definition for ‘lm’ -debias.met.regression: no visible global function definition for ‘coef’ -debias.met.regression: no visible global function definition for ‘vcov’ -debias.met.regression: no visible global function definition for - ‘terms’ -debias.met.regression: no visible global function definition for - ‘model.frame’ -debias.met.regression: no visible global function definition for - ‘model.matrix’ -debias.met.regression : : no visible global function - definition for ‘sd’ -debias.met.regression: no visible global function definition for - ‘rnorm’ -debias.met.regression: no visible global function definition for - ‘quantile’ -debias.met.regression: no visible binding for global variable - ‘quantile’ -debias.met.regression: no visible binding for global variable ‘Date’ -debias.met.regression: no visible binding for global variable ‘lwr’ -debias.met.regression: no visible binding for global variable ‘upr’ -debias.met.regression: no visible binding for global variable ‘obs’ debias.met.regression: no visible binding for global variable ‘values’ -debias.met.regression: no visible binding for global variable ‘Year’ -download.US_Wlef: no visible global function definition for - ‘read.table’ extract.local.CMIP5: no visible binding for global variable ‘GCM’ -extract.nc.ERA5 : : no visible global function definition - for ‘setNames’ -extract.nc.ERA5: no visible global function definition for ‘setNames’ extract.nc.ERA5 : : no visible binding for global variable ‘.’ -get.cruncep: no visible binding for global variable ‘Lat’ -get.cruncep: no visible binding for global variable ‘lati’ -get.cruncep: no visible global function definition for - ‘cruncep_dt2weather’ -get.weather: no visible global function definition for - ‘cruncep_dt2weather’ get_NARR_thredds: no visible binding for global variable ‘latitude’ get_NARR_thredds: no visible binding for global variable ‘longitude’ -is.land: no visible binding for global variable ‘met.nc’ -lm_ensemble_sims: no visible global function definition for ‘quantile’ -lm_ensemble_sims: no visible binding for global variable ‘mod.save’ -lm_ensemble_sims: no visible global function definition for ‘sd’ -load.cfmet: no visible binding for global variable ‘index’ -met.process: no visible global function definition for ‘get.id’ -met.process: no visible binding for global variable ‘site_id’ -met.process: no visible binding for global variable ‘format_id’ -met.process: no visible binding for global variable ‘machine_id’ met2CF.FACE: no visible binding for global variable ‘x’ -met_temporal_downscale.Gaussian_ensemble: no visible global function - definition for ‘sd’ -met_temporal_downscale.Gaussian_ensemble: no visible global function - definition for ‘rnorm’ met_temporal_downscale.Gaussian_ensemble: no visible binding for global variable ‘temp_max’ met_temporal_downscale.Gaussian_ensemble: no visible binding for global variable ‘temp_min’ -metgapfill.NOAA_GEFS: no visible global function definition for - ‘na.omit’ -metgapfill.NOAA_GEFS: no visible global function definition for ‘lm’ -metgapfill.NOAA_GEFS: no visible global function definition for - ‘predict’ -model.train: no visible global function definition for ‘lm’ -model.train: no visible global function definition for ‘coef’ -model.train: no visible global function definition for ‘vcov’ -model.train: no visible global function definition for ‘resid’ -subdaily_pred: no visible global function definition for ‘model.matrix’ Undefined global functions or variables: - . Date GCM I.diff I.dir Itot Lat Year aggregate air_pressure - air_temperature coef cruncep_dt2weather day doy format_id get.id hour - index lati latitude lm longitude lwr machine_id met.nc mod.save - model.frame model.matrix month na.omit obs predict pressure qmax qmin - quantile read.table relative_humidity resid rhmax rhmin rnorm sd - setNames site_id surface_downwelling_shortwave_flux_in_air temp_max - temp_min terms tmax tmin upr values vcov x year -Consider adding - importFrom("datasets", "pressure") - importFrom("stats", "aggregate", "coef", "lm", "model.frame", - "model.matrix", "na.omit", "predict", "quantile", "resid", - "rnorm", "sd", "setNames", "terms", "vcov") - importFrom("utils", "read.table") -to your NAMESPACE file. + . GCM latitude longitude temp_max temp_min values x * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -262,9 +77,6 @@ See chapter ‘Writing R documentation files’ in the ‘Writing R Extensions’ manual. * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... WARNING -Undocumented arguments in documentation object 'cfmet.downscale.time' - ‘lat’ ‘...’ - Undocumented arguments in documentation object 'closest_xy' ‘slat’ ‘slon’ ‘infolder’ ‘infile’ @@ -292,24 +104,15 @@ Undocumented arguments in documentation object 'download.NARR_site' Undocumented arguments in documentation object 'download.NEONmet' ‘...’ -Undocumented arguments in documentation object 'extract.nc.ERA5' - ‘...’ - Undocumented arguments in documentation object 'extract.nc' ‘...’ -Undocumented arguments in documentation object 'get.ncvector' - ‘var’ ‘lati’ ‘loni’ ‘run.dates’ - Undocumented arguments in documentation object 'lm_ensemble_sims' ‘lags.list’ Undocumented arguments in documentation object 'met.process' ‘browndog’ -Undocumented arguments in documentation object 'met.process.stage' - ‘input.id’ ‘raw.id’ - Undocumented arguments in documentation object 'met2CF.ALMA' ‘verbose’ @@ -319,23 +122,12 @@ Undocumented arguments in documentation object 'met2CF.Ameriflux' Undocumented arguments in documentation object 'met2CF.AmerifluxLBL' ‘...’ -Undocumented arguments in documentation object 'met2CF.FACE' - ‘in.path’ ‘in.prefix’ ‘outfolder’ ‘start_date’ ‘end_date’ ‘input.id’ - ‘site’ ‘format’ ‘...’ - -Undocumented arguments in documentation object 'met2CF.NARR' - ‘in.path’ ‘in.prefix’ ‘outfolder’ ‘...’ - Undocumented arguments in documentation object 'met2CF.PalEON' ‘lat’ ‘lon’ ‘verbose’ ‘...’ Undocumented arguments in documentation object 'met2CF.PalEONregional' ‘verbose’ ‘...’ -Undocumented arguments in documentation object 'met2CF.csv' - ‘in.path’ ‘in.prefix’ ‘outfolder’ ‘start_date’ ‘end_date’ ‘lat’ ‘lon’ - ‘overwrite’ ‘...’ - Undocumented arguments in documentation object 'metgapfill.NOAA_GEFS' ‘...’ @@ -354,12 +146,6 @@ Undocumented arguments in documentation object 'permute.nc' Undocumented arguments in documentation object 'predict_subdaily_met' ‘...’ -Undocumented arguments in documentation object 'site.lst' - ‘site.id’ ‘con’ - -Undocumented arguments in documentation object 'site_from_tag' - ‘sitename’ ‘tag’ - Undocumented arguments in documentation object 'temporal.downscale.functions' ‘...’ @@ -375,9 +161,6 @@ Argument items with no description in Rd object 'gen.subdaily.models': Argument items with no description in Rd object 'merge_met_variable': ‘start_date’ ‘end_date’ ‘...’ -Argument items with no description in Rd object 'met_temporal_downscale.Gaussian_ensemble': - ‘in.path’ ‘in.prefix’ - Argument items with no description in Rd object 'split_wind': ‘start_date’ ‘end_date’ @@ -390,7 +173,7 @@ Argument items with no description in Rd object 'split_wind': by using R CMD build --resave-data old_size new_size compress cruncep_landmask.RData 39Kb 9Kb xz - narr_cruncep_ebifarm.RData 790Kb 595Kb xz + narr_cruncep_ebifarm.RData 790Kb 597Kb xz * checking files in ‘vignettes’ ... WARNING Files in the 'vignettes' directory but no files in 'inst/doc': ‘ameriflux_demo.Rmd’, ‘cfmet_downscaling.Rmd’, @@ -398,9 +181,10 @@ Files in the 'vignettes' directory but no files in 'inst/doc': Package has no Sweave vignette sources and no VignetteBuilder field. * checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK -* checking tests ... SKIPPED +* checking tests ... OK + Running ‘test.NARR.R’ + Running ‘testthat.R’ * checking for non-standard things in the check directory ... OK * checking for detritus in the temp directory ... OK * DONE - -Status: 10 WARNINGs, 2 NOTEs +Status: 7 WARNINGs, 2 NOTEs diff --git a/modules/data.atmosphere/tests/testthat/test.cf-downscaling.R b/modules/data.atmosphere/tests/testthat/test.cf-downscaling.R index 2b514f84b6a..f2dce0296a3 100644 --- a/modules/data.atmosphere/tests/testthat/test.cf-downscaling.R +++ b/modules/data.atmosphere/tests/testthat/test.cf-downscaling.R @@ -11,18 +11,18 @@ test_that( "actual values will need to be revised if (when) algorithms change"), { b <- cfmet.downscale.time(cfmet = daily.cf, lat = 40) - expect_equal(b[,unique(year)], 1951) - expect_equal(b[,range(doy)], c(2,151)) - expect_equal(b[,unique(hour)], 0:23) - expect_equal(b[,round(range(downwelling_photosynthetic_photon_flux))], c(0, 2061)) - expect_equal(b[,round(range(air_temperature))], c(-22, 31)) - # expect_equal(b[,round(range(relative_humidity))], c(0.30569194491299, 1)) - expect_equal(b[,signif(range(precipitation_flux), 3)], c(0, 1.67e-05)) - expect_equal(b[,signif(range(wind), 2)], c(0.066, 6.60)) + expect_equal(unique(b$year), 1951) + expect_equal(range(b$doy), c(2,151)) + expect_equal(unique(b$hour), 0:23) + expect_equal(round(range(b$downwelling_photosynthetic_photon_flux)), c(0, 2061)) + expect_equal(round(range(b$air_temperature)), c(-22, 31)) + # expect_equal(round(range(b$relative_humidity)), c(0.30569194491299, 1)) + expect_equal(signif(range(b$precipitation_flux), 3), c(0, 1.67e-05)) + expect_equal(signif(range(b$wind), 2), c(0.066, 6.60)) }) test_that("downscaling with timestep", { - df <- data.table::data.table( + df <- data.frame( year = 2020, doy = 100, air_temperature_min = 293.15, air_temperature_max = 303.15, air_temperature = 298.15, surface_downwelling_shortwave_flux_in_air = 1000, @@ -43,13 +43,31 @@ test_that("downscaling with timestep", { purrr::walk(~{ expect_equal(mean(.$air_temperature), (df$air_temperature - 273.15)) # input is K, output is C expect_equal(sum(.$precipitation_flux), df$precipitation_flux) - expect_true(all(.$air_pressure == df$air_pressure)) + expect_true(all(.$wind == df$wind_speed)) }) }) +test_that("output for a given day not affected by adjacent days", { + df <- data.frame( + year = 2020, doy = 100:101, + air_temperature_min = 293.15 + c(0, 10), + air_temperature_max = 303.15 + c(0, 10), + air_temperature = 298.15 + c(0, 10), + surface_downwelling_shortwave_flux_in_air = c(1000, 2000), + air_pressure = 1030, + wind_speed = 0, + relative_humidity = 0.5, + precipitation_flux = c(0, 2 / (60 * 60))) + + # print(cfmet.downscale.daily(df[2,], 6, 40)) + # print(cfmet.downscale.daily(df, 6, 40)) + expect_equal(cfmet.downscale.daily(df[1,], 6, 40), cfmet.downscale.daily(df, 6, 40)[1:4,]) + expect_equal(cfmet.downscale.daily(df[2,], 6, 40), cfmet.downscale.daily(df, 6, 40)[5:8,]) +}) + test_that("get.ncvector works",{ - run.dates <- data.table::data.table(index = 1:2, date = c(lubridate::ymd("1951-01-01 UTC"), lubridate::ymd("1951-01-02 UTC"))) + run.dates <- data.frame(index = 1:2, date = c(lubridate::ymd("1951-01-01 UTC"), lubridate::ymd("1951-01-02 UTC"))) res <- get.ncvector("air_temperature", lati = 1, loni = 1, run.dates, met.nc = daily.nc) expect_type(res, "double") expect_equal(length(res), nrow(run.dates)) diff --git a/modules/data.atmosphere/tests/testthat/test.download.NARR.R b/modules/data.atmosphere/tests/testthat/test.download.NARR.R index a4547e144d6..cb4bd4121fa 100644 --- a/modules/data.atmosphere/tests/testthat/test.download.NARR.R +++ b/modules/data.atmosphere/tests/testthat/test.download.NARR.R @@ -6,10 +6,6 @@ ntime <- as.numeric(difftime(end_date, start_date) + 1) * 24 / 3 + 1 lat.in <- 43.3724 lon.in <- -89.9071 -outfolder <- tempfile() -dir.create(outfolder) -teardown(unlink(outfolder, recursive = TRUE, force = TRUE)) - test_url <- generate_narr_url(as.POSIXct("2000-01-01"), TRUE)[["url"]] test_nc <- tryCatch(ncdf4::nc_open(test_url), error = function(e) { skip("Unable to reach NARR server") @@ -21,15 +17,18 @@ test_that("NARR download works as expected", { # Please run locally to test! skip_on_ci() skip_if_offline() - r <- download.NARR_site(outfolder, start_date, end_date, lat.in, lon.in, - progress = TRUE, parallel = TRUE, ncores = 2) - - expect_equal(nrow(r), 1) - expect_true(file.exists(r$file[1])) - nc <- ncdf4::nc_open(r$file) - temp <- ncdf4::ncvar_get(nc, "air_temperature") - precip <- ncdf4::ncvar_get(nc, "precipitation_flux") - expect_true(all(!is.na(temp)), all(temp > 0), length(temp) == ntime) - expect_true(all(!is.na(precip)), length(precip) == ntime) - ncdf4::nc_close(nc) + withr::with_tempdir({ + r <- download.NARR_site( + outfolder = getwd(), + start_date, end_date, lat.in, lon.in, + progress = TRUE, parallel = TRUE, ncores = 2) + expect_equal(nrow(r), 1) + expect_true(file.exists(r$file[1])) + nc <- ncdf4::nc_open(r$file) + temp <- ncdf4::ncvar_get(nc, "air_temperature") + precip <- ncdf4::ncvar_get(nc, "precipitation_flux") + expect_true(all(!is.na(temp)), all(temp > 0), length(temp) == ntime) + expect_true(all(!is.na(precip)), length(precip) == ntime) + ncdf4::nc_close(nc) + }) }) diff --git a/modules/data.atmosphere/tests/testthat/test.load.cfmet.R b/modules/data.atmosphere/tests/testthat/test.load.cfmet.R index 42281fcb903..b0df59ff85b 100644 --- a/modules/data.atmosphere/tests/testthat/test.load.cfmet.R +++ b/modules/data.atmosphere/tests/testthat/test.load.cfmet.R @@ -14,7 +14,6 @@ on.exit(ncdf4::nc_close(subdaily.nc), add = TRUE) test_that("data extracted from test pecan-cf met files is valid",{ expect_is(daily.cf, "data.frame") - expect_is(daily.cf, "data.table") expect_is(daily.cf$date, "POSIXct") expect_is(daily.cf$date, "POSIXt")