Skip to content

Commit

Permalink
fix notes/warnings in data.land package
Browse files Browse the repository at this point in the history
  • Loading branch information
moki1202 committed Sep 15, 2023
1 parent ce327b9 commit 7056153
Show file tree
Hide file tree
Showing 15 changed files with 133 additions and 132 deletions.
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", .)
) %>%
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", .)
) %>%
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", .)) %>%
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", .)) %>%
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", .)) %>%
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", .)) %>%
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"
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::ines(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
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

0 comments on commit 7056153

Please sign in to comment.