Merge pull request #127 from ncss-tech/fixjindex

remove implicit conversion to SPDF with unit length j-index #125

R/SoilProfileCollection-methods.R

@@ -858,31 +858,38 @@ setMethod("[", signature=c("SoilProfileCollection", i="ANY", j="ANY"),
     # valid spatial data is now tested via validSpatialData()
     # also need to test that there is only 1 horizon/slice per location
   	# only produces a SPDF when j index is present
-    if(validSpatialData(x) & length(p.ids) == nrow(h) & !missing(j)) {
-      # combine with coordinates
-      message('result is a SpatialPointsDataFrame object')
-      # note that we are filtering based on 'i' - an index of selected profiles
-			
-      # since the order of our slices and coordinates are the same
-      # it is safe to use 'match.ID=FALSE'
-      # this gets around a potential problem when dimnames(x)[[1]] aren't consecutive 
-      # values-- often the case when subsetting has been performed
-
-      ## TODO: there should always be something in @site
-      # if site data, join hz+site
-      if(nrow(s) > 0) {
-      	return(SpatialPointsDataFrame(as(x, 'SpatialPoints')[i, ], data=join(h, s, by=idname(x)), match.ID=FALSE))
-      }
-      ## TODO: can this ever happen?
-      # no site data
-      else {
-      	return(SpatialPointsDataFrame(as(x, 'SpatialPoints')[i, ], data=h, match.ID=FALSE))	
-      }
-    }
+    # if(validSpatialData(x) & length(p.ids) == nrow(h) & !missing(j)) {
+    #   # combine with coordinates
+    #   message('result is a SpatialPointsDataFrame object')
+    #   # note that we are filtering based on 'i' - an index of selected profiles
+    # 
+    #   # since the order of our slices and coordinates are the same
+    #   # it is safe to use 'match.ID=FALSE'
+    #   # this gets around a potential problem when dimnames(x)[[1]] aren't consecutive 
+    #   # values-- often the case when subsetting has been performed
+    #   
+    #   ## TODO: there should always be something in @site
+    #   # if site data, join hz+site
+    #   if(nrow(s) > 0) {
+    #   	return(SpatialPointsDataFrame(as(x, 'SpatialPoints')[i, ], data=join(h, s, by=idname(x)), match.ID=FALSE))
+    #   }
+    #   ## TODO: can this ever happen?
+    #   # no site data
+    #   else {
+    #   	return(SpatialPointsDataFrame(as(x, 'SpatialPoints')[i, ], data=h, match.ID=FALSE))	
+    #   }
+    # }
 
     # in this case there may be missing coordinates, or we have more than 1 slice of hz data
-    else {
-      res <- SoilProfileCollection(idcol=idname(x), depthcols=horizonDepths(x), metadata=aqp::metadata(x), horizons=h, site=s, sp=sp, diagnostic=d, restrictions=r)
+    #else {
+      res <- SoilProfileCollection(idcol = idname(x), 
+                                   depthcols = horizonDepths(x), 
+                                   metadata = aqp::metadata(x), 
+                                   horizons = h, 
+                                   site = s, 
+                                   sp = sp, 
+                                   diagnostic = d, 
+                                   restrictions = r)
 
       # preserve one off slots that may have been customised relative to defaults 
       #  in prototype or resulting from construction of SPC 
@@ -904,7 +911,7 @@ setMethod("[", signature=c("SoilProfileCollection", i="ANY", j="ANY"),
 
 
       return(res)
-    }
+    #}
 
   # done
   }

R/glom.R

@@ -7,12 +7,12 @@
 # gloms a set of horizons for a single-profile SPC `p`
 #  the horizons are aggregated by depth using 
 #  clod.hz.ids() defined below
-glom <- function(p, z1, z2=NA, ids = FALSE, df = FALSE, truncate = FALSE) {
+glom <- function(p, z1, z2=NA, ids = FALSE, df = FALSE, truncate = FALSE, modality = "all") {
   # aka glom.by.depth; 
   if(length(p) > 1)
     stop("glom is intended for single-profile SPCs", call.=FALSE)
 
-  idx <- clod.hz.ids(p, z1, z2)
+  idx <- clod.hz.ids(p, z1, z2, modality)
 
   # short circuit to get hzIDs of intersection
   if(ids)
@@ -49,7 +49,7 @@ glom <- function(p, z1, z2=NA, ids = FALSE, df = FALSE, truncate = FALSE) {
 
 # returns unique index to all horizons occuring over the depth interval [z1, z2]. 
 # z2 is optional, in which case a single horizon with depth range containing z1 is returned
-clod.hz.ids <- function (p, z1, z2 = NA, as.list = FALSE) 
+clod.hz.ids <- function (p, z1, z2 = NA, modality = "all", as.list = FALSE) 
 {
   # access SPC slots to get important info about p
   hz <- horizons(p)
@@ -140,6 +140,13 @@ clod.hz.ids <- function (p, z1, z2 = NA, as.list = FALSE)
       }
     }
 
+    if(modality == "thickest") {
+      sub.thk <- bdep[idx.top:idx.bot] - tdep[idx.top:idx.bot]
+      max.sub.thk <- max(sub.thk)
+      first.thickest.idx <- which(sub.thk == max.sub.thk)[1]
+      idx.top <- idx.bot <- idx.top - 1 + first.thickest.idx
+    }
+
     # get the ID values out of horizon table
     idval <- hz[idx.top:idx.bot, hzid]
 
@@ -151,6 +158,11 @@ clod.hz.ids <- function (p, z1, z2 = NA, as.list = FALSE)
     return(list(hzid = hzid, hz.idx = idx.top:idx.bot, value = idval))
   }
 
+  if(modality == "thickest") {
+    first.thickest.idx <- which(bdep - tdep == max(bdep - tdep))[1]
+    idx.top <- first.thickest.idx
+  }
+
   idval <- hz[idx.top, hzid]
 
   if (!as.list) 

R/glomApply.R

@@ -5,18 +5,19 @@
 #' @param object A SoilProfileCollection
 #' @param .fun A function that returns vector with top and bottom depth (z1 and z2 arguments to \code{glom}) for a single profile `p` (as passed by \code{profileApply})
 #' @param truncate Truncate horizon top and bottom depths to z1 and z2? 
+#' @param modality Aggregation method for glom result. Default "all": all horizons; "thickest": return (shallowest) thickest horizon
 #' @param ... A set of comma-delimited R expressions that resolve to a transformation to be applied to a single profile e.g \code{glomApply(hzdept = max(hzdept) - hzdept)}
 #' @param chunk.size Chunk size parameter for \code{profileApply}
 #' @return A SoilProfileCollection.
 #' @author Andrew G. Brown.
 #' 
 #' @rdname glomApply
 #' @export glomApply
-glomApply <- function(object, .fun=NULL, truncate = FALSE, ..., chunk.size = 100) {
+glomApply <- function(object, .fun=NULL, truncate = FALSE, modality="all", ..., chunk.size = 100) {
   if(is.null(.fun) | !inherits(.fun, 'function'))
     stop("function `.fun`` to return glom boundaries for profiles is missing", call. = FALSE)
   aqp::union(profileApply(object, function(p, ...) {
     dep <- .fun(p, ...)
-    return(glom(p, dep[1], dep[2], truncate = truncate))
+    return(glom(p, dep[1], dep[2], truncate = truncate, modality = modality))
   }, simplify = FALSE, chunk.size = chunk.size))
 }

tests/testthat/test-SPC-objects.R

@@ -206,13 +206,21 @@ test_that("SPC spatial operations ", {
   sp1.spdf <- suppressMessages(as(sp1, 'SpatialPointsDataFrame'))
   expect_true(inherits(sp1.spdf, 'SpatialPointsDataFrame'))
 
+  # Unit-length j-index SPDF downgrading DEPRECATED
+
   # implicity down-grade to SPDF via hz-subsetting
-  sp1.spdf <- suppressMessages(sp1[, 1])
-  expect_true(inherits(sp1.spdf, 'SpatialPointsDataFrame'))
+  # sp1.spdf <- suppressMessages(sp1[, 1])
+  # expect_true(inherits(sp1.spdf, 'SpatialPointsDataFrame'))
+  # again, with profile indexing
+  # sp1.spdf <- suppressMessages(sp1[1, 1])
+  # expect_true(inherits(sp1.spdf, 'SpatialPointsDataFrame'))
 
+  # retain SPC object when using unit-length j index
+  sp1.spc <- suppressMessages(sp1[, 1])
+  expect_true(inherits(sp1.spc, 'SoilProfileCollection'))
   # again, with profile indexing
-  sp1.spdf <- suppressMessages(sp1[1, 1])
-  expect_true(inherits(sp1.spdf, 'SpatialPointsDataFrame'))
+  sp1.spc <- suppressMessages(sp1[1, 1])
+  expect_true(inherits(sp1.spc, 'SoilProfileCollection'))
 
 })