Merge pull request #76 from ncss-tech/horizon-id-madness

Horizon id madness

NAMESPACE

@@ -25,7 +25,9 @@ importFrom(methods,
   setMethod,
   as,
   show,
-  slot
+  slot,
+  .hasSlot,
+  slotNames
 )
 
 

NEWS

@@ -4,8 +4,12 @@
 * new functions: previewColors(), colorQuantiles(), plotColorQuantiles()
 * new function: horizonDepths()<-, edit top/bottom names after SPC init
 * new function: profile_id()<-, edit profile IDs after init; be careful!
-* rbind method for SoilProfileCollection objects [...]
+* new functions: hzID() and hzID()<-, get/set unique horizon IDs
+* new functions: hzidname() and hzidname()<-, get/set column containing unique horizon IDs
+* rbind.SoilProfileCollection() has been deprecated in favor of union(), gains new functionality: 
 * bug fixes in sanity checks for horizonNames()<-
+* !!! SoilProfileCollection internal structure has changed: 
+* 
 
 -------------------------- aqp 1.16-6 (2018-12-12) ------------------------
 * partial bug fix in test_hz_logic() related to missing top AND bottom depths, needs work: https://github.com/ncss-tech/aqp/issues/65

R/Class-SoilProfileCollection.R

@@ -1,3 +1,4 @@
+# https://github.com/ncss-tech/aqp/issues/75
 # class def for main class within aqp
 .SoilProfileCollectionValidity <- function(object) {
 
@@ -43,6 +44,7 @@ setClass(
   Class='SoilProfileCollection', 
   representation=representation(
     idcol='character', # column name containing IDs
+    hzidcol='character',
     depthcols='character', # 2 element vector with column names for hz top, bottom
     metadata='data.frame', # single-row dataframe with key-value mapping
     horizons='data.frame', # all horizons sorted by ID, top
@@ -52,6 +54,7 @@ setClass(
   ),
   prototype=prototype(
     idcol='id',
+    hzidcol='hzID',
     depthcols=c('top','bottom'),
     metadata=data.frame(stringsAsFactors=FALSE), # default units are unkown
     horizons=data.frame(stringsAsFactors=FALSE),

R/SoilProfileCollection-coercion.R

@@ -1,4 +1,37 @@
 ## Coercition methods: general
+
+# safely deconstruct as list
+setAs("SoilProfileCollection", "list", function(from) {
+
+  # get slot names from prototype
+  sn <- slotNames(from)
+
+  # test for presence of all slots
+  # copy contents over to list with same name
+  # if missing return NULL + warning
+  s.list <- lapply(sn, function(i) {
+    if(.hasSlot(from, name=i)) {
+      res <- slot(from, i)
+    } else {
+      res <- NULL
+    }
+    return(res)
+  })
+
+  # copy slot names
+  names(s.list) <- sn
+
+  # test for missing slots
+  if(any(sapply(s.list, is.null))) {
+    warning("some slots were missing, use reBuildSPC to fix", call. = FALSE)
+  }
+
+  return(s.list)
+
+}
+)
+
+
 setAs("SoilProfileCollection", "data.frame", function(from) {
 
   # horizons + site + coordinates

R/SoilProfileCollection-methods.R

@@ -67,6 +67,28 @@ setMethod("idname", "SoilProfileCollection",
     return(object@idcol)
 )
 
+## horizon ID name
+if (!isGeneric("hzidname"))
+  setGeneric("hzidname", function(object, ...) standardGeneric("hzidname"))
+
+setMethod("hzidname", "SoilProfileCollection",
+          function(object)
+            return(object@hzidcol)
+)
+
+## get horizon IDs
+if (!isGeneric("hzID"))
+  setGeneric("hzID", function(object, ...) standardGeneric("hzID"))
+
+setMethod("hzID", "SoilProfileCollection",
+          function(object) {
+            h <- horizons(object)
+            res <- h[[hzidname(object)]]
+            return(res)
+          }
+
+)
+
 
 ## distinct profile IDs
 if (!isGeneric("profile_id"))
@@ -189,103 +211,6 @@ setMethod("horizonNames", "SoilProfileCollection",
 ## overloads
 ##
 
-### This will be greatly improved with new class structure
-## concatentation
-## # https://github.com/ncss-tech/aqp/issues/71
-## TODO: concatenation of data with duplicated IDs in @site, but unique data in other @site fields, will result in corrupt SPC
-## TODO: duplicates in @sp will cause errors
-## TODO: duplicates are removed in all other slots... does this make sense?
-rbind.SoilProfileCollection <- function(...) {
-	# setup some defaults
-	options(stringsAsFactors=FALSE)
-
-	# parse dots
-	objects <- list(...)
-	names(objects) <- NULL
-
-	# short-circuits
-	if(length(objects) == 0)
-		return(NULL)
-	if(length(objects) == 1)
-		return(objects[1])
-
-	## TODO: normalize idname and horizonDepths
-	# profile_id() <- 
-	# horizonDepths() <- 
-
-
-	# combine pieces
-	# should have length of 1
-	o.idname <- unique(lapply(objects, idname))
-	o.depth.units <- unique(lapply(objects, depth_units))
-	o.hz.depths <- unique(lapply(objects, horizonDepths))
-	o.m <- unique(lapply(objects, aqp::metadata))
-	o.coords <- unique(lapply(objects, function(i) ncol(coordinates(i))))
-	o.p4s <- unique(lapply(objects, proj4string))
-
-	# should have length > 1
-	o.h <- lapply(objects, horizons)
-	o.s <- lapply(objects, site)
-	o.d <- lapply(objects, diagnostic_hz)
-	o.sp <- lapply(objects, slot, 'sp')
-
-	# sanity checks:
-	if(length(o.idname) > 1)
-		stop('inconsistent ID names', call.=FALSE)
-	if(length(o.depth.units) > 1)
-		stop('inconsistent depth units', call.=FALSE)
-	if(length(o.hz.depths) > 1)
-		stop('inconsistent depth columns', call.=FALSE)
-	if(length(o.m) > 1)
-		stop('inconsistent metadata', call.=FALSE)
-
-	# spatial data may be missing...
-	if(length(o.coords) > 1)
-		stop('inconsistent spatial data', call.=FALSE)
-	if(length(o.p4s) > 1)
-		stop('inconsistent CRS', call.=FALSE)
-
-	# generate new SPC components
-	# using plyr::rbind.fill seems to solve the problem on non-conformal DF
-	# is it safe?
-	# https://github.com/ncss-tech/aqp/issues/71
-	o.h <- unique(do.call('rbind.fill', o.h)) # horizon data
-	o.s <- unique(do.call('rbind.fill', o.s)) # site data
-	o.d <- unique(do.call('rbind.fill', o.d)) # diagnostic data, leave as-is
-
-	## 2015-12-18: removed re-ordering, was creating corrupt SPC objects
-	##             site and horizon data
-
-	# spatial points require some more effort when spatial data are missing
-	o.1.sp <- objects[[1]]@sp
-	if(ncol(coordinates(o.1.sp)) == 1) # missing spatial data
-		o.sp <- o.1.sp # copy the first filler
-
-	## 2015-12-18: added call to specific function: "sp::rbind.SpatialPoints"
-	# not missing spatial data
-	else
-		o.sp <- do.call("rbind.SpatialPoints", o.sp) # rbind properly
-
-	# make SPC and return
-	res <- SoilProfileCollection(idcol=o.idname[[1]], depthcols=o.hz.depths[[1]], metadata=o.m[[1]], horizons=o.h, site=o.s, sp=o.sp, diagnostic=o.d)
-
-# 	# one more final check:
-# 	print(profile_id(res))
-# 	print( site(res)[[idname(res)]])
-# 	print(site(res))
-
-	if(length(profile_id(res)) != length(site(res)[[idname(res)]]))
-	  stop("SPC object corruption. This shouldn't happen and will be fixed in aqp 2.0", call. = FALSE)
-	if(! all.equal(profile_id(res), site(res)[[idname(res)]]))
-	  stop("SPC object corruption. This shouldn't happen and will be fixed in aqp 2.0", call. = FALSE)
-
-	return(res)
-}
-
-
-## TODO: this doesn't work as expected ... fix in 2.0
-## overload rbind
-#setMethod("rbind", "SoilProfileCollection", .rbind.SoilProfileCollection)
 
 
 

R/SoilProfileCollection-slice-methods.R

@@ -153,6 +153,7 @@ slice.fast <- function(object, fm, top.down=TRUE, just.the.data=FALSE, strict=TR
   if(just.the.data)
     return(hd.slices)
 
+  ## TODO: WTF (AGB: loafercreek[, 2])
   # if spatial data and only a single slice: SPDF
   if(nrow(coordinates(object)) == length(object) & length(z) == 1) {
     cat('result is a SpatialPointsDataFrame object\n')
@@ -169,7 +170,15 @@ slice.fast <- function(object, fm, top.down=TRUE, just.the.data=FALSE, strict=TR
 
 
   # otherwise return an SPC, be sure to copy over the spatial data
-  depths(hd.slices) <- as.formula(paste(id, '~', top, '+', bottom))
+  # NOTE: suppressing warning due to non-unique horizon IDs, don't panic
+  suppressWarnings(depths(hd.slices) <- as.formula(paste(id, '~', top, '+', bottom)))
+
+  # reset auto-generated horizon ID so that we know it is now the slice ID
+  idx <- match(hzidname(hd.slices), horizonNames(hd.slices))
+  horizonNames(hd.slices)[idx] <- 'sliceID'
+  hzidname(hd.slices) <- 'sliceID'
+
+  # copy spatial data
   hd.slices@sp <- object@sp
 
   # if site data: return an SPC + @site

R/checkSPC.R

@@ -0,0 +1,23 @@
+
+# test for valid SPC, based on presence / absense of slots as compared to 
+# class prototype
+# likely only used between major versions of aqp where internal structure of SPC has changed
+checkSPC <- function(x) {
+
+  # get slot names from prototype
+  sn <- slotNames(x)
+
+  # test for all slots in the prototype
+  s.test <- sapply(sn, function(i) .hasSlot(x, name=i))
+
+  # a valid object will have all slots present
+  if(all(s.test)) {
+    res <- TRUE
+  } else {
+    res <- FALSE
+  }
+
+  return(res)
+}
+
+

R/rebuildSPC.R

@@ -0,0 +1,24 @@
+# repair an SPC by breaking into pieces and re-assembling
+# likely only used to fix outdated SPC objects that are missing slots
+rebuildSPC <- function(x) {
+
+  # break into pieces as list
+  x.list <- suppressWarnings(as(x, 'list'))
+
+  # seed object for new SPC
+  res <- x.list$horizons
+
+  # init SPC from pieces
+  # note: using depths<- because it will generate a horizon ID
+  fm <- as.formula(sprintf("%s ~ %s + %s", x.list$idcol, x.list$depthcols[1], x.list$depthcols[2]))
+  depths(res) <- fm
+
+  # add additional pieces
+  metadata(res) <- x.list$metadata
+  site(res) <- x.list$site
+  res@sp <- x.list$sp
+  diagnostic_hz(res) <- x.list$diagnostic
+
+  return(res)
+}
+