Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

fix notes/warnings in data.land package #3219

Merged
merged 22 commits into from
Feb 14, 2024
Merged
Show file tree
Hide file tree
Changes from 20 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion docker/depends/pecan_package_dependencies.csv
Original file line number Diff line number Diff line change
Expand Up @@ -478,7 +478,6 @@
"RPostgres","*","base/db","Suggests",FALSE
"RPostgreSQL","*","base/db","Suggests",FALSE
"RPostgreSQL","*","models/biocro","Suggests",FALSE
"RPostgreSQL","*","modules/data.land","Suggests",FALSE
"Rpreles","*","models/preles","Suggests",FALSE
"RSQLite","*","base/db","Suggests",FALSE
"sf","*","modules/assim.sequential","Suggests",FALSE
Expand Down
1 change: 0 additions & 1 deletion modules/data.land/DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,6 @@ Suggests:
PEcAn.settings,
redland,
raster,
RPostgreSQL,
testthat (>= 1.0.2)
License: BSD_3_clause + file LICENSE
Copyright: Authors
Expand Down
28 changes: 14 additions & 14 deletions modules/data.land/R/IC_BADM_Utilities.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ Read.IC.info.BADM <-function(lat, long){
#Reading in the DB
#
U.S.SB <-
read.csv(system.file("data","BADM.csv", package = "PEcAn.data.land"))
utils::read.csv(system.file("data","BADM.csv", package = "PEcAn.data.land"))


Regions <- EPA_ecoregion_finder(lat, long)
Expand All @@ -31,8 +31,8 @@ Read.IC.info.BADM <-function(lat, long){
#L2
biomass.df <- U.S.SB %>%
dplyr::filter(
NA_L2CODE == Code_Level,
VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)
.data$NA_L2CODE == Code_Level,
.data$VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)
infotroph marked this conversation as resolved.
Show resolved Hide resolved
) %>%
dplyr::select("SITE_ID", "GROUP_ID", "VARIABLE_GROUP", "VARIABLE", "DATAVALUE")

Expand All @@ -43,8 +43,8 @@ Read.IC.info.BADM <-function(lat, long){

biomass.df <- U.S.SB %>%
dplyr::filter(
NA_L1CODE == Code_Level,
VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)
.data$NA_L1CODE == Code_Level,
.data$VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)
infotroph marked this conversation as resolved.
Show resolved Hide resolved
) %>%
dplyr::select("SITE_ID", "GROUP_ID", "VARIABLE_GROUP", "VARIABLE", "DATAVALUE")
}
Expand All @@ -54,7 +54,7 @@ Read.IC.info.BADM <-function(lat, long){
if (nrow(biomass.df) < 3) {
Code_Level <- "ALL"
biomass.df <- U.S.SB %>%
dplyr::filter(VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)) %>%
dplyr::filter(.data$VARIABLE %>% grepl("ROOT_|AG_BIOMASS|SOIL_STOCK|SOIL_CHEM", .)) %>%
infotroph marked this conversation as resolved.
Show resolved Hide resolved
dplyr::select("SITE_ID", "GROUP_ID", "VARIABLE_GROUP", "VARIABLE", "DATAVALUE")
}

Expand Down Expand Up @@ -94,8 +94,8 @@ Read.IC.info.BADM <-function(lat, long){
type <- type[-which(type == "*_BIOMASS")]
#----------------- Unit conversion
unit.in <- Gdf %>%
dplyr::filter(VARIABLE %>% grepl("UNIT", .)) %>%
dplyr::pull(DATAVALUE) %>%
dplyr::filter(.data$VARIABLE %>% grepl("UNIT", .)) %>%
infotroph marked this conversation as resolved.
Show resolved Hide resolved
dplyr::pull(.data$DATAVALUE) %>%
as.character()


Expand All @@ -113,8 +113,8 @@ Read.IC.info.BADM <-function(lat, long){
unit.ready <- "kg/m^2"

Date.in <- Gdf %>%
dplyr::filter(VARIABLE %>% grepl("DATE", .)) %>%
dplyr::pull(DATAVALUE) %>%
dplyr::filter(.data$VARIABLE %>% grepl("DATE", .)) %>%
infotroph marked this conversation as resolved.
Show resolved Hide resolved
dplyr::pull(.data$DATAVALUE) %>%
as.Date()

if (length(Date.in) == 0)
Expand All @@ -124,8 +124,8 @@ Read.IC.info.BADM <-function(lat, long){
# if it's biomass
if (type == "*_BIOMASS") {
Oregan.in <- Gdf %>%
dplyr::filter(VARIABLE %>% grepl("ORGAN", .)) %>%
dplyr::pull(DATAVALUE)
dplyr::filter(.data$VARIABLE %>% grepl("ORGAN", .)) %>%
infotroph marked this conversation as resolved.
Show resolved Hide resolved
dplyr::pull(.data$DATAVALUE)


PlantWoodIni <-
Expand All @@ -134,8 +134,8 @@ Read.IC.info.BADM <-function(lat, long){

} else if (type == "*SOIL") {
val <- Gdf %>%
dplyr::filter(VARIABLE %>% grepl("SOIL_STOCK_C_ORG", .)) %>% #"SOIL_STOCK_C_ORG"
dplyr::pull(DATAVALUE) %>%
dplyr::filter(.data$VARIABLE %>% grepl("SOIL_STOCK_C_ORG", .)) %>% #"SOIL_STOCK_C_ORG"
infotroph marked this conversation as resolved.
Show resolved Hide resolved
dplyr::pull(.data$DATAVALUE) %>%
as.numeric()

if (length(val) > 0)
Expand Down
14 changes: 7 additions & 7 deletions modules/data.land/R/InventoryGrowthFusion.R
Original file line number Diff line number Diff line change
Expand Up @@ -404,8 +404,8 @@ model{

## JAGS initial conditions
init <- list()
if(is.mcmc.list(restart)){
init <- mcmc.list2init(restart)
if(coda::is.mcmc.list(restart)){
init <- PEcAn.utils::mcmc.list2init(restart)
nchain <- length(init)
} else {
nchain <- 3
Expand Down Expand Up @@ -438,7 +438,7 @@ model{
}

PEcAn.logger::logger.info("RUN MCMC")
load.module("dic")
rjags::load.module("dic")
for(k in avail.chunks){

## determine whether to sample states
Expand All @@ -457,17 +457,17 @@ model{
save(jags.out,file=ofile)

## update restart
if(!is.null(restart) & ((is.logical(restart) && restart) || is.mcmc.list(restart))){
if(!is.null(restart) & ((is.logical(restart) && restart) || coda::is.mcmc.list(restart))){
ofile <- paste("IGF",model,"RESTART.RData",sep=".")
jags.final <- coda.samples(model = j.model, variable.names = c("x",out.variables), n.iter = 1)
jags.final <- rjags::coda.samples(model = j.model, variable.names = c("x",out.variables), n.iter = 1)
k_restart = k + 1 ## finished k, so would restart at k+1
save(jags.final,k_restart,file=ofile)
}

## check for convergence and break from loop early
D <- as.mcmc.list(lapply(jags.out,function(x){x[,'deviance']}))
D <- coda::as.mcmc.list(lapply(jags.out,function(x){x[,'deviance']}))
gbr <- coda::gelman.diag(D)$psrf[1,1]
trend <- mean(sapply(D,function(x){coef(lm(x~seq_len(n.chunk)))[2]}))
trend <- mean(sapply(D,function(x){stats::coef(stats::lm(x~seq_len(n.chunk)))[2]}))
if(gbr < 1.005 & abs(trend) < 0.5) break
}

Expand Down
76 changes: 38 additions & 38 deletions modules/data.land/R/InventoryGrowthFusionDiagnostics.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,13 @@ InventoryGrowthFusionDiagnostics <- function(jags.out, combined=NULL) {
out <- as.matrix(jags.out)
x.cols <- which(substr(colnames(out), 1, 1) == "x")
if(length(x.cols) > 0){
ci <- apply(out[, x.cols], 2, quantile, c(0.025, 0.5, 0.975))
ci <- apply(out[, x.cols], 2, stats::quantile, c(0.025, 0.5, 0.975))
ci.names <- parse.MatrixNames(colnames(ci), numeric = TRUE)

### DBH par(mfrow=c(3,2))
if(length(x.cols) > 0){
layout(matrix(1:8, 4, 2, byrow = TRUE))
ci <- apply(out[, x.cols], 2, quantile, c(0.025, 0.5, 0.975))
graphics::layout(matrix(1:8, 4, 2, byrow = TRUE))
ci <- apply(out[, x.cols], 2, stats::quantile, c(0.025, 0.5, 0.975))
ci.names <- parse.MatrixNames(colnames(ci), numeric = TRUE)

smp <- sample.int(data$ni, min(8, data$ni))
Expand All @@ -26,29 +26,29 @@ InventoryGrowthFusionDiagnostics <- function(jags.out, combined=NULL) {
plot(data$time, ci[2, sel], type = "n",
ylim = range(rng), ylab = "DBH (cm)", main = i)
PEcAn.visualization::ciEnvelope(data$time, ci[1, sel], ci[3, sel], col = "lightBlue")
points(data$time, data$z[i, ], pch = "+", cex = 1.5)
graphics::points(data$time, data$z[i, ], pch = "+", cex = 1.5)
# lines(data$time,z0[i,],lty=2)

## growth
sel <- which(ci.names$row == i)
inc.mcmc <- apply(out[, x.cols[sel]], 1, diff)
inc.ci <- apply(inc.mcmc, 1, quantile, c(0.025, 0.5, 0.975)) * 5
inc.ci <- apply(inc.mcmc, 1, stats::quantile, c(0.025, 0.5, 0.975)) * 5
# inc.names = parse.MatrixNames(colnames(ci),numeric=TRUE)

plot(data$time[-1], inc.ci[2, ], type = "n",
ylim = range(inc.ci, na.rm = TRUE), ylab = "Ring Increment (mm)")
PEcAn.visualization::ciEnvelope(data$time[-1], inc.ci[1, ], inc.ci[3, ], col = "lightBlue")
points(data$time, data$y[i, ] * 5, pch = "+", cex = 1.5, type = "b", lty = 2)
graphics::points(data$time, data$y[i, ] * 5, pch = "+", cex = 1.5, type = "b", lty = 2)
}
}
}

if (FALSE) {
## check a DBH
plot(out[, which(colnames(out) == "x[3,31]")])
abline(h = z[3, 31], col = 2, lwd = 2)
hist(out[, which(colnames(out) == "x[3,31]")])
abline(v = z[3, 31], col = 2, lwd = 2)
graphics::abline(h = z[3, 31], col = 2, lwd = 2)
graphics::hist(out[, which(colnames(out) == "x[3,31]")])
graphics::abline(v = z[3, 31], col = 2, lwd = 2)
}

## process model
Expand All @@ -59,105 +59,105 @@ InventoryGrowthFusionDiagnostics <- function(jags.out, combined=NULL) {
grep("alpha",colnames(out)),
grep("deviance",colnames(out)))]

par(mfrow = c(1, 1))
graphics::par(mfrow = c(1, 1))
for (i in vars) {
hist(out[, i], main = colnames(out)[i])
abline(v=0,lwd=3)
graphics::hist(out[, i], main = colnames(out)[i])
graphics::abline(v=0,lwd=3)
}
if (length(vars) > 1 && length(vars) < 10) {
pairs(out[, vars])
graphics::pairs(out[, vars])
}

if("deviance" %in% colnames(out)){
hist(out[,"deviance"])
graphics::hist(out[,"deviance"])
vars <- c(vars,which(colnames(out)=="deviance"))
}


## rebuild coda for just vars
var.out <- as.mcmc.list(lapply(jags.out,function(x){ x[,vars]}))
var.out <- coda::as.mcmc.list(lapply(jags.out,function(x){ x[,vars]}))

## convergence
gelman.diag(var.out)
coda::gelman.diag(var.out)

#### Diagnostic plots
plot(var.out)

if("deviance" %in% colnames(out)){
hist(out[,"deviance"])
graphics::hist(out[,"deviance"])
vars <- c(vars,which(colnames(out)=="deviance"))
}

## rebuild coda for just vars
var.out <- as.mcmc.list(lapply(jags.out,function(x){ x[,vars]}))
var.out <- coda::as.mcmc.list(lapply(jags.out,function(x){ x[,vars]}))

## convergence
gelman.diag(var.out)
coda::gelman.diag(var.out)

#### Diagnostic plots
plot(var.out)


## Standard Deviations layout(matrix(c(1,2,3,3),2,2,byrow=TRUE))
par(mfrow = c(2, 3))
graphics::par(mfrow = c(2, 3))
prec <- out[, grep("tau", colnames(out))]
for (i in seq_along(colnames(prec))) {
hist(1 / sqrt(prec[, i]), main = colnames(prec)[i])
graphics::hist(1 / sqrt(prec[, i]), main = colnames(prec)[i])
}
cor(prec)
stats::cor(prec)
# pairs(prec)

### alpha
par(mfrow = c(1, 1))
graphics::par(mfrow = c(1, 1))
alpha.cols <- grep("alpha", colnames(out))
if (length(alpha.cols) > 0) {
alpha.ord <- 1:length(alpha.cols)
ci.alpha <- apply(out[, alpha.cols], 2, quantile, c(0.025, 0.5, 0.975))
ci.alpha <- apply(out[, alpha.cols], 2, stats::quantile, c(0.025, 0.5, 0.975))
plot(alpha.ord, ci.alpha[2, ], type = "n",
ylim = range(ci.alpha, na.rm = TRUE), ylab = "Random Effects")
PEcAn.visualization::ciEnvelope(alpha.ord, ci.alpha[1, ], ci.alpha[3, ], col = "lightBlue")
lines(alpha.ord, ci.alpha[2, ], lty = 1, lwd = 2)
abline(h = 0, lty = 2)
graphics::lines(alpha.ord, ci.alpha[2, ], lty = 1, lwd = 2)
graphics::abline(h = 0, lty = 2)
}

par(mfrow = c(1, 1))
graphics::par(mfrow = c(1, 1))
### alpha
alpha.cols <- grep("alpha", colnames(out))
if (length(alpha.cols) > 0) {
alpha.ord <- 1:length(alpha.cols)
ci.alpha <- apply(out[, alpha.cols], 2, quantile, c(0.025, 0.5, 0.975))
ci.alpha <- apply(out[, alpha.cols], 2, stats::quantile, c(0.025, 0.5, 0.975))
plot(alpha.ord, ci.alpha[2, ], type = "n",
ylim = range(ci.alpha, na.rm = TRUE), ylab = "Random Effects")
PEcAn.visualization::ciEnvelope(alpha.ord, ci.alpha[1, ], ci.alpha[3, ], col = "lightBlue")
lines(alpha.ord, ci.alpha[2, ], lty = 1, lwd = 2)
abline(h = 0, lty = 2)
graphics::lines(alpha.ord, ci.alpha[2, ], lty = 1, lwd = 2)
graphics::abline(h = 0, lty = 2)
}

### YEAR
year.cols <- grep("year", colnames(out))
if (length(year.cols > 0)) {
ci.yr <- apply(out[, year.cols], 2, quantile, c(0.025, 0.5, 0.975))
ci.yr <- apply(out[, year.cols], 2, stats::quantile, c(0.025, 0.5, 0.975))
plot(data$time, ci.yr[2, ], type = "n",
ylim = range(ci.yr, na.rm = TRUE), ylab = "Year Effect")
PEcAn.visualization::ciEnvelope(data$time, ci.yr[1, ], ci.yr[3, ], col = "lightBlue")
lines(data$time, ci.yr[2, ], lty = 1, lwd = 2)
abline(h = 0, lty = 2)
graphics::lines(data$time, ci.yr[2, ], lty = 1, lwd = 2)
graphics::abline(h = 0, lty = 2)
}

### INDIV
ind.cols <- which(substr(colnames(out), 1, 3) == "ind")
if (length(ind.cols) > 0 & !is.null(combined)) {
boxplot(out[, ind.cols], horizontal = TRUE, outline = FALSE, col = as.factor(combined$PLOT))
abline(v = 0, lty = 2)
graphics::boxplot(out[, ind.cols], horizontal = TRUE, outline = FALSE, col = as.factor(combined$PLOT))
graphics::abline(v = 0, lty = 2)
tapply(apply(out[, ind.cols], 2, mean), combined$PLOT, mean)
table(combined$PLOT)

spp <- combined$SPP
# boxplot(out[order(spp),ind.cols],horizontal=TRUE,outline=FALSE,col=spp[order(spp)])
boxplot(out[, ind.cols], horizontal = TRUE, outline = FALSE, col = spp)
abline(v = 0, lty = 2)
graphics::boxplot(out[, ind.cols], horizontal = TRUE, outline = FALSE, col = spp)
graphics::abline(v = 0, lty = 2)
spp.code <- levels(spp)[table(spp) > 0]
legend("bottomright", legend = rev(spp.code), col = rev(which(table(spp) > 0)), lwd = 4)
graphics::legend("bottomright", legend = rev(spp.code), col = rev(which(table(spp) > 0)), lwd = 4)
tapply(apply(out[, ind.cols], 2, mean), combined$SPP, mean)
}
} # InventoryGrowthFusionDiagnostics
Expand Down
1 change: 0 additions & 1 deletion modules/data.land/R/Soilgrids_SoilC_prep.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@
#' @return A data frame containing AGB median and sd for each site and each time step.
#' @export
#'
#' @examples
#' @author Dongchen Zhang
#' @importFrom magrittr %>%
Soilgrids_SoilC_prep <- function(site_info, start_date, end_date, time_points,
Expand Down
6 changes: 2 additions & 4 deletions modules/data.land/R/download_NEON_soilmoisture.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
##' @name download_NEON_soilmoist
##' @description:
##' Download NEON Soil Water Content and Soil Salinity data by date and site name
##'
##' @param site four letter NEON site code name(s). If no site is specified, it will download all of them (chr) (e.g "BART" or c("SRER", "KONA", "BART"))
Expand All @@ -17,8 +15,8 @@
##'
##' @author Juliette Bateman
##'
##' @example
##' \run{
##' @examples
##' \dontrun{
##' test <- download_NEON_soilmoisture(
##' site = c("SRER", "BART", "KONA"),
##' avg = 30,
Expand Down
6 changes: 3 additions & 3 deletions modules/data.land/R/extract_FIA.R
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ extract_FIA <- function(lon, lat, start_date, end_date, gridres = 0.075, dbparms
veg_info <- list()

fia.con <- PEcAn.DB::db.open(dbparms$fia)
on.exit(db.close(fia.con), add = TRUE)
on.exit(PEcAn.DB::db.close(fia.con), add = TRUE)

lonmin <- lon - gridres
lonmax <- lon + gridres
Expand Down Expand Up @@ -49,7 +49,7 @@ extract_FIA <- function(lon, lat, start_date, end_date, gridres = 0.075, dbparms
" AND p.lat <= ", latmax, " AND p.measyear >= ", min.year,
" AND p.measyear <= ", max.year, " GROUP BY p.cn")

plot.info <- db.query(query, con = fia.con)
plot.info <- PEcAn.DB::db.query(query, con = fia.con)
if (nrow(plot.info) == 0) {
PEcAn.logger::logger.severe("No plot data found on FIA.")
}
Expand Down Expand Up @@ -97,7 +97,7 @@ extract_FIA <- function(lon, lat, start_date, end_date, gridres = 0.075, dbparms
" and p.lon < ", lonmax,
" and p.lat >= ", latmin,
" and p.lat < ", latmax)
tree.info <- db.query(query, con = fia.con)
tree.info <- PEcAn.DB::db.query(query, con = fia.con)
names(tree.info) <- tolower(names(tree.info))

if (nrow(tree.info) == 0) {
Expand Down
Loading
Loading