Cran 1.0.2 (#99)

* Fix Atlas Issue Ref #94 checks if using Atlas and whether the shape is invalid if so it sets the precision * Close #96 - expanded warning message * Update NEWS * Bump package version
DistanceDevelopment · Nov 17, 2023 · 154a9db · 154a9db
1 parent 8a4b590
commit 154a9db
Show file tree

Hide file tree

Showing 11 changed files with 256 additions and 189 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -11,7 +11,7 @@ Suggests: knitr,
 VignetteBuilder: knitr
 Type: Package
 Title: Distance Sampling Survey Design
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(
   person("Laura", "Marshall", email = "[email protected]", role = c("aut", "cre")),
   person("Rexstad", "Eric", email = "[email protected]", role = "ctb"))
@@ -59,6 +59,7 @@ Collate:
     'generate.systematic.points.R'
     'get.intersection.points.R'
     'line.coords.as.dataframe.R'
+    'mat.mult.R'
     'point.coords.as.dataframe.R'
     'run.coverage.R'
     'write.transects.R'
diff --git a/NAMESPACE b/NAMESPACE
@@ -45,11 +45,15 @@ importFrom(methods,new)
 importFrom(methods,show)
 importFrom(methods,validObject)
 importFrom(sf,as_Spatial)
+importFrom(sf,st_is_valid)
+importFrom(sf,st_make_valid)
+importFrom(sf,st_set_precision)
 importFrom(sf,st_zm)
 importFrom(stats,median)
 importFrom(stats,na.omit)
 importFrom(stats,rbinom)
 importFrom(stats,runif)
 importFrom(stats,sd)
+importFrom(utils,sessionInfo)
 importFrom(utils,write.csv)
 importFrom(utils,write.table)
diff --git a/NEWS b/NEWS
@@ -1,3 +1,14 @@
+dssd 1.0.2
+----------
+
+Bug Fixes
+
+ * Fixed Atlas bug by setting precision in Atlas environment when invalid shapes are detected.
+
+Enhancements
+
+ * Enhanced message regarding strata names when they are provided but do not match the number of strata. 
+
 dssd 1.0.1
 ----------
 

diff --git a/R/Class.Constructors.R b/R/Class.Constructors.R
@@ -99,9 +99,9 @@ make.region <- function(region.name = "region",
   }else if(strata.count != length(strata.name)){
     if(length(sf.shape) <= 26){
       strata.name <- LETTERS[1:length(sf.shape[[sf.column]])]
-      warning("Automatically naming strata as no strata names provided. Assigned strata names: ", paste(strata.name, collapse = ", "), call. = F, immediate. = T)
+      warning("Automatically naming strata as no (or incorrect number of) strata names provided. Assigned strata names: ", paste(strata.name, collapse = ", "), call. = F, immediate. = T)
     }else{
-      stop("Too many strata (>26) for strata names to be assigned default names, please provide strata names.", call. = FALSE)
+      stop("Too many strata (>26) for strata names to be assigned default names, please provide the correct number of strata names.", call. = FALSE)
     }
   }
   # Check the format of the shape

diff --git a/R/Line.Transect.Design.R b/R/Line.Transect.Design.R
diff --git a/R/calc.region.width.R b/R/calc.region.width.R
@@ -19,6 +19,15 @@ calc.region.width <- function(design, strata.id = NULL){
     theta <- ifelse(rot.angle.rad == 0, 0, 2*pi-rot.angle.rad)
     rot.mat <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2, byrow = FALSE)
     rot.strata <- strata*rot.mat
+    # if we are using atlas and the shape is not valid
+    if(grepl("atlas", sessionInfo()$BLAS) && is.na(sf::st_is_valid(rot.strata))){
+      # turn it into and sfc shape
+      tmp <- sf::st_sfc(rot.strata)
+      # make valid with setting the precision
+      tmp <- sf::st_make_valid(sf::st_set_precision(tmp,1e8))
+      # extract shape again
+      rot.strata <- tmp[[1]]
+    }
     #Find the width of the region
     bbox <- sf::st_bbox(rot.strata)
     width <- width + (bbox$xmax - bbox$xmin)

diff --git a/R/generate.eqspace.zigzags.R b/R/generate.eqspace.zigzags.R
@@ -11,7 +11,7 @@ generate.eqspace.zigzags <- function(design, strata.id, samplers, line.length, s
   rot.angle.rad <- rot.angle.rad/180*pi
   theta <- ifelse(rot.angle.rad == 0, 0, 2*pi-rot.angle.rad)
   rot.mat <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2, byrow = FALSE)
-  rot.strata <- strata*rot.mat
+  rot.strata <- mat.mult(strata, rot.mat)
   #Buffer strata for plus sampling?
   if(design@edge.protocol[strata.id] == "plus"){
     rot.strata <- sf::st_buffer(rot.strata, design@truncation)
@@ -207,7 +207,6 @@ generate.eqspace.zigzags <- function(design, strata.id, samplers, line.length, s
   #Rotate back again
   reverse.theta <- rot.angle.rad
   rot.mat.rev <- matrix(c(cos(reverse.theta), sin(reverse.theta), -sin(reverse.theta), cos(reverse.theta)), ncol = 2, byrow = FALSE)
-  mat.mult <- function(x,y){return(x*y)}
   lines.unrotated <- lapply(to.keep, mat.mult, y=rot.mat.rev)
   transects <- lines.unrotated
   #Also rotate covered region

diff --git a/R/generate.parallel.lines.R b/R/generate.parallel.lines.R
@@ -10,7 +10,7 @@ generate.parallel.lines <- function(design, strata.id, samplers, line.length, sp
   rot.angle.rad <- design@design.angle[strata.id]/180*pi
   theta <- ifelse(rot.angle.rad == 0, 0, 2*pi-rot.angle.rad)
   rot.mat <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2, byrow = FALSE)
-  rot.strata <- strata*rot.mat
+  rot.strata <- mat.mult(strata, rot.mat)
   #Buffer strata for plus sampling?
   if(design@edge.protocol[strata.id] == "plus"){
     rot.strata <- st_buffer(rot.strata, design@truncation)
@@ -136,7 +136,6 @@ generate.parallel.lines <- function(design, strata.id, samplers, line.length, sp
   #Rotate back again
   reverse.theta <- rot.angle.rad
   rot.mat.rev <- matrix(c(cos(reverse.theta), sin(reverse.theta), -sin(reverse.theta), cos(reverse.theta)), ncol = 2, byrow = FALSE)
-  mat.mult <- function(x,y){return(x*y)}
   lines.unrotated <- lapply(to.keep, mat.mult, y=rot.mat.rev)
   transects <- lines.unrotated
   #Also rotate covered region

diff --git a/R/generate.segmented.grid.R b/R/generate.segmented.grid.R
@@ -12,7 +12,7 @@ generate.segmented.grid <- function(design, strata.id, samplers, line.length, sp
   rot.angle.rad <- rot.angle/180*pi
   theta <- ifelse(rot.angle.rad == 0, 0, 2*pi-rot.angle.rad)
   rot.mat <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2, byrow = FALSE)
-  rot.strata <- strata*rot.mat
+  rot.strata <- mat.mult(strata, rot.mat)
   #Buffer strata for plus sampling?
   if(design@edge.protocol[strata.id] == "plus"){
     rot.strata <- st_buffer(rot.strata, design@truncation)
@@ -163,7 +163,6 @@ generate.segmented.grid <- function(design, strata.id, samplers, line.length, sp
   #Rotate back again
   reverse.theta <- rot.angle.rad
   rot.mat.rev <- matrix(c(cos(reverse.theta), sin(reverse.theta), -sin(reverse.theta), cos(reverse.theta)), ncol = 2, byrow = FALSE)
-  mat.mult <- function(x,y){return(x*y)}
   lines.unrotated <- lapply(to.keep, mat.mult, y=rot.mat.rev)
   transects <- lines.unrotated
   #Also rotate covered region

diff --git a/R/generate.systematic.points.R b/R/generate.systematic.points.R
@@ -1,3 +1,4 @@
+#' @importFrom utils sessionInfo
 generate.systematic.points <- function(design, strata.id, spacing, samplers, coverage.grid = FALSE, calc.cov.area = TRUE, clip.to.strata = TRUE, quiet = FALSE){
   #Generates systematic points
   region <- design@region
@@ -9,7 +10,7 @@ generate.systematic.points <- function(design, strata.id, spacing, samplers, cov
   rot.angle.rad <- design@design.angle[strata.id]/180*pi
   theta <- ifelse(rot.angle.rad == 0, 0, 2*pi-rot.angle.rad)
   rot.mat <- matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2, byrow = FALSE)
-  rot.strata <- strata*rot.mat
+  rot.strata <- mat.mult(strata, rot.mat)
   #Buffer strata for plus sampling?
   if(design@edge.protocol[strata.id] == "plus"){
     rot.strata <- sf::st_buffer(rot.strata, design@truncation)
@@ -59,7 +60,6 @@ generate.systematic.points <- function(design, strata.id, spacing, samplers, cov
     #Rotate back again
     reverse.theta <- rot.angle.rad
     rot.mat.rev <- matrix(c(cos(reverse.theta), sin(reverse.theta), -sin(reverse.theta), cos(reverse.theta)), ncol = 2, byrow = FALSE)
-    mat.mult <- function(x,y){return(x*y)}
     points.unrotated <- lapply(points.inside, mat.mult, y=rot.mat.rev)
     transects <- points.unrotated
   }

diff --git a/R/mat.mult.R b/R/mat.mult.R
@@ -0,0 +1,37 @@
+#' @importFrom utils sessionInfo
+#' @importFrom sf st_make_valid st_set_precision st_is_valid
+mat.mult <- function(x,y){
+  # Internal function for performing matrix multiplication with atlas checks for validity
+  # reverse rotation
+  unrotate <- x*y
+  # if we are using atlas and the shape is not valid
+  if(grepl("atlas", sessionInfo()$BLAS) && is.na(sf::st_is_valid(unrotate))){
+    # validity flag
+    is.valid <- FALSE
+    # counter
+    index <- 1
+    # Precision values to try
+    precision <- c(1e10,1e9,1e8,1e7,1e6,1e5,1e4)
+    # turn it into and sfc shape
+    tmp <- sf::st_sfc(unrotate)
+    # Keep trying to fix across a range of precision values
+    while(!is.valid && index <= length(precision)){
+      tmp2 <- try(sf::st_make_valid(sf::st_set_precision(tmp, precision[index])), silent = TRUE)
+      if(inherits(tmp2, "try-error")){
+        # If the problem is not fixed try a lower precision
+        index <- index + 1
+      }else{
+        # Otherwise it is fixed continue
+        is.valid <- TRUE
+        tmp <- tmp2
+      }
+    }
+    # What if none of the precisions succeed?
+    if(!is.valid){
+      stop("Problem with matrix multiplication due to BLAS/ATLAS setup. dssd has tried to resolve this issue by reducing the precision to 1e4 (see ?sf::st_precision) but differences persist in what should be identical coordinates.", call. = FALSE)
+    }
+    # extract shape again
+    unrotate <- tmp[[1]]
+  }
+  return(unrotate)
+}