update SA functions

USDAForestService · May 29, 2024 · 0b660bb · 0b660bb
1 parent 162ce6f
commit 0b660bb
Show file tree

Hide file tree

Showing 7 changed files with 75 additions and 59 deletions.
diff --git a/R/DBgetPlots.R b/R/DBgetPlots.R
@@ -3464,10 +3464,14 @@ DBgetPlots <- function (states = NULL,
                     message(seed.qry)
                     return(NULL) })
         }
+        if (is.null(seedx)) {
+          message("no seedling data for ", stabbr)
+          message(seed.qry)
+        }
         if (!is.null(seedx) && nrow(seedx) != 0 && length(ssvars) > 0) {
-	      if (!"seed" %in% names(dbqueries[[state]])) {
+	        if (!"seed" %in% names(dbqueries[[state]])) {
             dbqueries[[state]]$seed <- seed.qry
-	      }
+	        }
 
           seedx <- setDT(seedx)
           seedx[, PLT_CN := as.character(PLT_CN)]

diff --git a/R/ISAinternal.R b/R/ISAinternal.R
@@ -18,7 +18,7 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
 
   ## global parameters
   stepcnt <- 1
-  maxbndxlst <-{}
+  maxbndxlst=maxbndx_intdlst <-{}
   int.pct=helperbndx.tmp <- NULL
 
   mbndlst <- list()
@@ -177,6 +177,11 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
 #        maxbndxlst <- maxbnd.gtthres
 #      } else if (multiSAdoms) {
       if (multiSAdoms) {
+        ## Create list of new maxbnd(s)
+        maxbndx_intdlst <- lapply(maxbndxlst, function(mbnd, maxbndx.intd) {
+          maxbndx.intd[maxbndx.intd[[maxbnd.unique]] %in% mbnd, ]
+        }, maxbndx.intd)
+
         if (byeach) {
           mbndlst <- maxbnd_max[[maxbnd.unique]]
           sbndlst <- lapply(maxbnd_max[[smallbnd.unique]], 
@@ -195,7 +200,7 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
               smallbndx[smallbndx[[smallbnd.unique]] %in% sbnd.att,]
           }, maxbnd_max, smallbndx, smallbnd.unique)
           names(sbndlst) <- mbndlst
-        
+
 
           ## Appends small areas from maxbnds less than threshold to maxbnds greater 
           ## than threshold based on closest centroids
@@ -245,10 +250,12 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
  #         }, maxbnd_max, smallbndx, smallbnd.unique)
 
           sbndlst <- list(smallbndx)
+          maxbndx_intdlst <- list(maxbndx.intd)
 
         } else {
           mbndlst <- maxbndxlst[1]
           sbndlst <- list(smallbndx)
+          maxbndx_intdlst <- list(maxbndx.intd)
         }
       }
     } else {
@@ -332,6 +339,7 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
   #################################################################
   helperbndxlst <- list()
   smallbndxlst <- list()
+  largebndxdlst <- list()
   SAbndlst <- list()
   SAdomsnmlst <- vector("character", length(sbndlst))
 
@@ -376,7 +384,7 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
     helperbndx <- helperbndx[, helperbnd.unique]
 
     j <- 1
-    ## Loop thru maxbndlst
+    ## Loop thru maxbndxlst
     while (nbrdom < nbrdom.minx && j <= ifelse(length(maxbndxlst) > 0, length(mbnd), 1)) {
       if (length(mbnd) > 0) {
         message("\nadding ", maxbnd.unique, ": ", mbnd[j])
@@ -457,6 +465,8 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
 
       ## Select largebnd(s) that intersect more than threshold
       largebnd_select <- largebndx.intd[largebndx.intd[[largebnd.unique]] %in% largebnd.gtthres,]
+      largebndxdlst[[mbnd]] <- largebnd_select
+
 
       ############################################################
       ## Display and save image for largebnd_intersect
@@ -511,33 +521,20 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
 
       ## For largebnds with no intersecting smallbnd, get largebndx polygons closest
       ## to smallbnd centroid within maxbnd
-      if (length(largebndxlst) > 0) {
+      if (length(largebndxlst) > 1) {
         largebndx.dist <- closest_poly(sbnd.centroid,
 			     ypoly=largebndx[largebndx[[largebnd.unique]] %in% largebndxlst,],
 			     ypoly.att=largebnd.unique, returnsf=FALSE)
-           largebndxlst <- unique(c(largebnd.ltthres, largebnd.lt0, names(largebndx.dist)))
+        largebndxlst <- unique(c(largebnd.ltthres, largebnd.lt0, names(largebndx.dist)))
       }
 
       while (!end) {
         ## Get intersecting helper polygons
-        helperbndx.tmp <- sf::st_join(helperbndx,
-						sf_dissolve(largebnd_select, largebnd.unique),
-						join=sf::st_intersects, left=FALSE, largest=TRUE)
-
-        #sf_use_s2(TRUE)
-        #helperbndx.tmp2 <- s2::s2_intersects(helperbndx,
-        #                              sf_dissolve(largebnd_select, largebnd.unique))
-        #helperbndx.tmp2
-
-        # get percent overlap of helperbndx.int and y largebndx.int
-        ############################################################
-#        helperbndx.tmp$FID <- seq(1:nrow(helperbndx.tmp))
-#        helperbndx.tmppct <- suppressWarnings(tabulateIntersections(layer1=helperbndx.tmp,
-#			layer1fld="FID", layer2=largebnd_select))
-#        FIDpct <- helperbndx.tmppct$FID[helperbndx.tmppct$int.pct > largebnd.threshold]
-#        helperbndx.tmp <- helperbndx.tmp[helperbndx.tmp$FID %in% FIDpct,]
-
-        helperbndx.tmppct <- tryCatch(suppressWarnings(tabulateIntersections(layer1=helperbndx.tmp,
+        helperbndx.tmp <- suppressWarnings(selectByIntersects(helperbndx, 
+                              sf_dissolve(largebnd_select, largebnd.unique), .1))
+
+        helperbndx.tmppct <- tryCatch(
+            suppressWarnings(tabulateIntersections(layer1=helperbndx.tmp,
 			            layer1fld=helperbnd.unique, layer2=largebnd_select)),
      	 			error=function(e) {
 					message("helperbnd intersection...")
@@ -687,7 +684,7 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
     if (polyunion && !bayes) {
       ## Remove columns in helperbndx.tmp with same names as in smallbnd attributes
       helperbndx.tmp <- helperbndx.tmp[, names(helperbndx.tmp)[!names(helperbndx.tmp) %in%
-			c(smallbnd.unique, "AOI", smallbnd.domain)]]
+		 	              c(smallbnd.unique, "AOI", smallbnd.domain)]]
 
       ## Change name of helperbnd.unique values if the same as smallbnd.unique values
       if (any(helperbndx.tmp[[helperbnd.unique]] %in% smallbndx[[smallbnd.unique]])) {
@@ -709,16 +706,17 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
     ## Add AOI to SAdoms
     if (bayes) {
       SAdoms$AOI <- 1
+      if (length(maxbnd.gtthres) > 1) {
+        ## Subset maxbndxd to only maxbnd.unique that intersects with smallbnd
+        maxgtthres2.intd <- maxbndxd[maxbndxd[[maxbnd.unique]] %in% 
+                                      maxbnd.gtthres[!maxbnd.gtthres %in% mbndlst],]
+        helperbndx2.tmp <- suppressWarnings(selectByIntersects(helperbndx, maxgtthres2.intd, 1))
+        helperbndx2.tmp <- suppressWarnings(selectByIntersects(helperbndx2.tmp, sbnd, .1)) 
+        helperbndx2.tmp$AOI <- 1
+        SAdoms <- rbind(SAdoms, helperbndx2.tmp)
+      }
       setnames(SAdoms, helperbnd.unique, "DOMAIN")
-
-      helperbndx.test <- helperbndx.tmp[helperbndx.tmp$PROVINCE == 'M331',]
-      test <- sf::st_intersection(sbnd[, smallbnd.unique], 
-                                  helperbndx.test[, helperbnd.unique])
 
-
-      sbnd <- sf::st_cast(sf::st_intersection(sbnd[, smallbnd.unique], 
-                      helperbndx.tmp[, helperbnd.unique]))
-      setnames(sbnd, helperbnd.unique, "DOMAIN")
     } else {
       ## Add 0 to non-AOI
       SAdoms[is.na(SAdoms$AOI), "AOI"] <- 0
@@ -743,9 +741,16 @@ helper.select <- function(smallbndx, smallbnd.unique, smallbnd.domain=NULL,
 
   names(SAdomslst) <- SAdomsnmlst
   names(smallbndxlst) <- SAdomsnmlst
-
+  #names(SAdomslst) <- mbndlst
+  #names(smallbndxlst) <- mbndlst
+  names(helperbndxlst) <- SAdomsnmlst
+  names(largebndxdlst) <- SAdomsnmlst
+  if (!is.null(maxbndx_intdlst)) {
+    names(maxbndx_intdlst) <- SAdomsnmlst
+  }
+
   returnlst <- list(SAdomslst=SAdomslst, helperbndxlst=helperbndxlst,
-				smallbndxlst=smallbndxlst)
+				smallbndxlst=smallbndxlst, largebndxlst=largebndxdlst, maxbndxlst=maxbndx_intdlst)
 
   if (!is.null(maxbndx)) {
     returnlst$maxbndx_intersect <- maxbndx_intersect[maxbndx.pct[[maxbnd.unique]] %in% maxbndxlst, ]

diff --git a/R/modSAarea.R b/R/modSAarea.R
@@ -787,7 +787,8 @@ modSAarea <- function(SApopdatlst = NULL,
 			    showsteps=showsteps, savesteps=savesteps,
 			    stepfolder=stepfolder, prior=prior, 
 			    modelselect=modelselect, multest=multest,
-			    SApackage=SApackage, SAmethod=SAmethod),
+			    SApackage=SApackage, SAmethod=SAmethod, bayes=bayes,
+			    save4testing=FALSE, vars2keep=vars2keep),
      	      error=function(e) {
 			      message("error with estimates of ", response, "...")
 			      message(e, "\n")
@@ -882,7 +883,8 @@ modSAarea <- function(SApopdatlst = NULL,
 				showsteps=showsteps, savesteps=savesteps, 
 				stepfolder=stepfolder, prior=prior, 
 				modelselect=modelselect, multest=multest,
-				SApackage=SApackage, SAmethod=SAmethod), 
+				SApackage=SApackage, SAmethod=SAmethod, bayes=bayes,
+				vars2keep=vars2keep), 
 			 error=function(e) {
 			 message("error with estimates of ", response, "...")
 			 message(e, "\n")

diff --git a/R/modSApop.R b/R/modSApop.R
@@ -584,9 +584,9 @@ modSApop <- function(popType = "VOL",
   if (!"AOI" %in% names(pltassgnx)) {
     pltassgnx$AOI <- 1
   }
-  
+
   if (!is.null(pvars2keep)) {
-    pvars2keep <- pvars2keep[pvars2keep %in% names(pltx)]
+    pvars2keep <- pvars2keep[pvars2keep %in% names(pltx) & !pvars2keep %in% names(pltassgnx)]
     if (length(pvars2keep) > 0) {
       pltassgnx <- merge(pltassgnx, pltx[, c(puniqueid, pvars2keep), with=FALSE], 
                        by.x=pltassgnid, by.y=puniqueid)

diff --git a/R/modSAtree.R b/R/modSAtree.R
@@ -855,7 +855,7 @@ modSAtree <- function(SApopdatlst = NULL,
                stepfolder=stepfolder, prior=prior, 
                modelselect=modelselect, multest=multest,
                SApackage=SApackage, SAmethod=SAmethod, bayes=bayes, # TODO: pass bayes_opts
-               save4testing=TRUE, vars2keep=vars2keep),
+               save4testing=FALSE, vars2keep=vars2keep),
         error=function(e) {
           message("error with estimates of ", response, "...")
           message(e, "\n")
@@ -906,11 +906,11 @@ modSAtree <- function(SApopdatlst = NULL,
       ## Merge SAdom attributes to dunit_totest
       if (addSAdomsdf) {
         pdomdat <- merge(setDT(SAdomsdf)[, 
-			unique(c("DOMAIN", "AOI", SAdomvars)), with=FALSE], 
-			pdomdat, by=c("DOMAIN", "AOI"))
+			       unique(c("DOMAIN", "AOI", SAdomvars)), with=FALSE], 
+			       pdomdat, by=c("DOMAIN", "AOI"))
         dunitlut <- merge(setDT(SAdomsdf)[, 
-			unique(c("DOMAIN", "AOI", SAdomvars)), with=FALSE], 
-			dunitlut, by=c("DOMAIN", "AOI"))
+			       unique(c("DOMAIN", "AOI", SAdomvars)), with=FALSE], 
+			       dunitlut, by=c("DOMAIN", "AOI"))
       }
       pdomdatlst[[SApopdatnm]] <- pdomdat
       dunitlutlst[[SApopdatnm]] <- dunitlut
@@ -936,7 +936,7 @@ modSAtree <- function(SApopdatlst = NULL,
 #largebnd.val=largebnd.vals
 #domain=rowcolinfo$rowvar
 #largebnd.unique="LARGEBND"
- 
+
       dunit_rowestlst <- 
 		tryCatch(
 			lapply(largebnd.vals, SAest.large, 
@@ -948,7 +948,8 @@ modSAtree <- function(SApopdatlst = NULL,
 				showsteps=showsteps, savesteps=savesteps,
 				stepfolder=stepfolder, prior=prior, 
 				modelselect=modelselect, multest=multest, 
-				SApackage=SApackage, SAmethod=SAmethod),
+				SApackage=SApackage, SAmethod=SAmethod, bayes=bayes,
+				vars2keep=vars2keep),
      	 	error=function(e) {
 			message("error with estimates of ", response, " by ", rowvar, "...")
 			message(e, "\n")

diff --git a/R/spGetAuxiliary.R b/R/spGetAuxiliary.R
@@ -393,18 +393,22 @@ spGetAuxiliary <- function(xyplt = NULL,
 
   ## Check continuous rasters
   ###################################################################
-  rastlst.contfn <- tryCatch(
+  if (!is.null(rastlst.cont)) {
+    rastlst.contfn <- tryCatch(
               getrastlst(rastlst.cont, rastfolder, quiet=TRUE, gui=gui),
      	 	            error=function(e) {
 			              message(e, "\n")
 			              return("stop") })
+  }
   if (!is.null(rastlst.contfn)) {
     if (length(rastlst.contfn) == 1) {
       if (rastlst.contfn == "stop") {
-        stop()
+        stop("invalid rastlst.contfn: \n",
+             toString(sapply(rastlst.cont, normalizePath)))
       }
     }
   }
+
   if (!is.null(rastlst.contfn)) {
     band.cont <- sapply(rastlst.contfn, function(x) rasterInfo(x)$nbands)
     nlayers.cont <- sum(band.cont)
@@ -460,18 +464,21 @@ spGetAuxiliary <- function(xyplt = NULL,
 
   ## Check categorical rasters
   ###################################################################
-  rastlst.catfn <- tryCatch(
+  if (!is.null(rastlst.cat)) {
+    rastlst.catfn <- tryCatch(
              getrastlst(rastlst.cat, rastfolder, quiet=TRUE, gui=gui),
      	 	           error=function(e) {
 			               message(e, "\n")
 			             return("stop") })
+  }
   if (is.null(rastlst.contfn) && is.null(rastlst.catfn)) {
     message("both rastlst.cont and rastlst.cat are NULL")
   }
   if (!is.null(rastlst.catfn)) {
     if (length(rastlst.catfn) == 1) {
       if (rastlst.catfn == "stop") {
-        stop()
+        stop("invalid rastlst.catfn: \n",
+             toString(sapply(rastlst.cat, normalizePath)))
       }
     }
     band.cat <- sapply(rastlst.catfn, function(x) rasterInfo(x)$nbands)
@@ -622,7 +629,6 @@ spGetAuxiliary <- function(xyplt = NULL,
     polyvarlst <- c("UNITVAR", vars2keep)
   }
 
-
   #############################################################################
   ## 2) Set up outputs - unitzonal, prednames, inputdf, zonalnames
   #############################################################################

diff --git a/R/spGetSAdoms.R b/R/spGetSAdoms.R
@@ -527,9 +527,7 @@ spGetSAdoms <- function(smallbnd,
   largebndx <- pcheck.spatial(layer=largebnd, dsn=largebnd_dsn, 
 		caption="large boundary")
   if (!all(sf::st_is_valid(largebndx))) {
-    largebndx <- sf::st_make_valid(largebndx, 
-                                 geos_method = 'valid_structure', 
-                                 geos_keep_collapsed = FALSE)
+    largebndx <- sf::st_make_valid(largebndx)
   }
 
   ## Check largebndx
@@ -573,10 +571,8 @@ spGetSAdoms <- function(smallbnd,
   #############################################################################
   helperbndx <- pcheck.spatial(layer=helperbnd, dsn=helperbnd_dsn, 
 		caption="helper boundary")
-  if (!all(sf::st_is_valid(smallbndx))) {
-    helperbndx <- sf::st_make_valid(helperbndx, 
-                                  geos_method = 'valid_structure', 
-                                  geos_keep_collapsed = FALSE)
+  if (!all(sf::st_is_valid(helperbndx))) {
+    helperbndx <- sf::st_make_valid(helperbndx)
   }
 
   if (is.null(largebndx)) {
@@ -693,6 +689,8 @@ spGetSAdoms <- function(smallbnd,
     SAdomslst <- autoselectlst$SAdomslst
     helperbndxlst <- autoselectlst$helperbndxlst
     smallbndxlst <- autoselectlst$smallbndxlst
+    largebndxlst <- autoselectlst$largebndxlst
+    maxbndxlst <- autoselectlst$maxbndxlst
 
   } else {