R/cellMetrics_speciesRaster.R

Defines functions convertNAtoEmpty cellMetrics_speciesRaster

Documented in cellMetrics_speciesRaster

##' @title Cell Metrics
##' 
##' @description Calculate various morphological and phylogenetic community
##' metrics for every cell in a \code{speciesRaster} object. 
##' 
##' @param x object of class \code{speciesRaster}
##' 
##' @param metric name of metric to use, see Details. 
##'
##' @param var If a univariate morphological metric is specified, and the 
##' 	data in \code{x} are multivariate, which trait should be used?
##' 	This can also specify which subset of columns a multivariate metric should be applied to.
##' 
##' @param nreps Number of repetitions for Foote metric distribution.
##'
##' @param verbose Intended primarily for debugging, prints progress to the console
##' 
##' @return object of class \code{speciesRaster} where the raster represents
##' 	calculations of the metric at every cell. 
##' 
##' @details 
##' 	Univariate trait metrics
##' 	\itemize{
##' 		\item{mean}
##' 		\item{median}
##' 		\item{range}
##'			\item{mean_NN_dist:} {mean nearest neighbor distance}
##'			\item{min_NN_dist:} {minimum nearest neighbor distance}
##' 		\item{variance}
##' 		\item{arithmeticWeightedMean} (see below)
##' 		\item{geometricWeightedMean} (see below)
##'			\item{phylosignal:} {Blomberg's K for phylogenetic signal}
##' 	}
##' 	Multivariate trait metrics
##'		\itemize{
##'			\item{disparity} 
##' 		\item{range}
##' 		\item{rangePCA}
##' 		\item{mean_NN_dist:} {mean nearest neighbor distance}
##'			\item{min_NN_dist:} {minimum nearest neighbor distance}
##' 		\item{phylosignal:} {Blomberg's K for phylogenetic signal, as implemented in geomorph::physignal}
##' 	}
##' 	Phylogenetic metrics
##' 	\itemize{
##'			\item{pd:} {Faith's phylogenetic diversity, including the root}
##'			\item{meanPatristic}
##'			\item{meanPatristicNN:} {mean nearest neighbor in patristic distance}
##'			\item{minPatristicNN:} {minimum nearest neighbor in patristic distance}
##'			\item{phyloDisparity:} {sum of squared deviations in patristic distance}
##'			\item{PSV:} {Phylogenetic Species Variability}
##'			\item{DR:} {non-parametric estimate of speciation rates}
##' 	}
##' 	Range-weighted metrics
##' 	\itemize{
##'			\item{weightedEndemism}
##'			\item{correctedWeightedEndemism:} {Weighted endemism standardized by species richness}
##'			\item{phyloWeightedEndemism:}
##' 	}
##'
##'		If data slot contains a pairwise matrix, \code{var} is ignored.
##'		Weighted mean options are available where, for each cell, a weighting scheme 
##' 	(inverse of species range sizes) is applied such that small-ranged species are 
##' 	up-weighted, and broadly distributed species are down-weighted. 
##' 	This can be a useful way to lessen the influence of broadly distributed species 
##' 	in the geographic mapping of trait data. 
##'
##' @examples
##' tamiasSpRas <- addPhylo_speciesRaster(tamiasSpRas, tamiasTree)
##' tamiasSpRas <- addTraits_speciesRaster(tamiasSpRas, tamiasTraits)
##'
##' # univariate morphological example
##' x <- cellMetrics_speciesRaster(tamiasSpRas, metric='mean', var='V2')
##' plot(x) 
##'
##' # multivariate morphological
##' x <- cellMetrics_speciesRaster(tamiasSpRas, metric='disparity')
##' plot(x)
##' 
##' # phylogenetic metrics
##' x <- cellMetrics_speciesRaster(tamiasSpRas, metric='meanPatristic')
##' plot(x)
##' 
##' @export


cellMetrics_speciesRaster <- function(x, metric, var = NULL, nreps = 20, verbose = FALSE) {
	
	if (!inherits(x, 'speciesRaster')) {
		stop('x must be of class speciesRaster.')
	}
	
	if (!identical(names(x), c('raster', 'speciesList', 'cellCommInd', 'geogSpecies', 'cellCount', 'data', 'phylo'))) {
		stop('speciesRaster object is missing some components.')
	}

	if (length(metric) > 1) {
		stop('You can only specify one metric.')
	}
	
	metric <- match.arg(metric, choices = c('mean', 'median', 'range', 'variance', 'arithmeticWeightedMean', 'geometricWeightedMean', 'rangePCA', 'disparity', 'mean_NN_dist', 'min_NN_dist', 'pd', 'meanPatristic', 'meanPatristicNN', 'minPatristicNN', 'phyloDisparity', 'PSV', 'DR', 'weightedEndemism', 'correctedWeightedEndemism', 'phyloWeightedEndemism', 'phylosignal'))
	
	pairwise <- FALSE
	
	# if data is pairwise matrix, then set flag appropriately
	if (inherits(x[['data']], c('matrix', 'data.frame'))) {
		if (identical(rownames(x[['data']]), colnames(x[['data']]))) {
			if (verbose) message('\t...detected pairwise distance matrix...\n') 
			var <- NULL
			pairwise <- TRUE
			# make the diagonal and lower triangle NA
			x[['data']][lower.tri(x[['data']], diag = TRUE)] <- NA
			if (all(is.na(x[['data']][upper.tri(x[['data']])]))) {
				stop('There are no values in the upper triangle of the pairwise matrix.')
			}
		}
	}
	
	# if var is defined but data are univariate, disable var. 
	if (is.vector(x[['data']]) & !is.null(var)) {
		var <- NULL
	}
		
	if (!is.null(var) & inherits(x[['data']], c('matrix', 'data.frame'))) {
		if (is.character(var)) {
			if (any(!var %in% colnames(x[['data']]))) {
				stop('var not a valid column name of the data.')
			}
		} else if (is.numeric(var)) {
			if (any(!all(var %in% 1:ncol(x[['data']])))) {
				stop('Requested data column indices are out of range.')
			}
		}
	}
	
	if (inherits(x[['data']], c('matrix', 'data.frame')) | is.vector(x[['data']])) {
		if (inherits(x[['data']], c('matrix', 'data.frame'))) {
			metricType <- 'multiVar'
			if (!is.null(var) & length(var) == 1) {
				metricType <- 'uniVar'
			}
		} else {
			metricType <- 'uniVar'
		}
	} else {
		metricType <- 'none'
	}
	
	# if a subset of data columns are requested, subset the data table
	if (!is.null(var) & inherits(x[['data']], c('matrix', 'data.frame'))) {
		if (length(var) > 1) {
			x[['data']] <- x[['data']][, var]
		} else {
			x[['data']] <- setNames(x[['data']][, var], rownames(x[['data']]))
		}
		var <- NULL
	}
	
	if (metric %in% c('mean', 'median', 'variance', 'arithmeticWeightedMean', 'geometricWeightedMean') & metricType == 'multiVar' & !pairwise) {
		stop('If a univariate metric is requested from a multivariate dataset, a column name must be provided as var.')
	}
	
	if (metric %in% c('rangePCA', 'disparity') & metricType == 'uniVar' & !pairwise) {
		stop('A multivariate metric cannot be applied to a univariate dataset.')
	}
	
	if (metric == 'geometricWeightedMean') {
		if (any(x[['data']] < 0)) {
			stop('Negative trait values cannot be supplied to geometricWeightedMean because the log of a negative number is NA.')
		}
	}
	
	# For metrics that require both trait and phylogenetic data, prune the data and tree to the shared taxon set
	if (metric == 'phylosignal') {
		if (is.vector(x[['data']])) {
			sharedTaxa <- intersect(names(x[['data']]), x[['phylo']]$tip.label)
			x[['data']] <- x[['data']][sharedTaxa]
		} else {
			sharedTaxa <- intersect(rownames(x[['data']]), x[['phylo']]$tip.label)
			x[['data']] <- x[['data']][sharedTaxa,]
		}
		x[['phylo']] <- ape::keep.tip(x[['phylo']], sharedTaxa)
	}
	
	# Prune species list according to metric of interest
	if (metric %in% c('mean', 'median', 'disparity', 'range', 'variance', 'arithmeticWeightedMean', 'geometricWeightedMean', 'rangePCA', 'mean_NN_dist', 'min_NN_dist', 'phylosignal')) {
		
		if (verbose) message('\t...dropping species that are not in trait data...\n')
		# check that there is data in speciesRaster object
		if (is.null(x[['data']])) {
			stop('speciesRaster object does not contain trait data!')
		}
		
		# prune speciesRaster object down to species shared with data
		if (is.vector(x[['data']])) {
			x[['speciesList']] <- intersectList(x[['speciesList']], names(x[['data']]))
		} else {
			x[['speciesList']] <- intersectList(x[['speciesList']], rownames(x[['data']]))
		}
		
	 } else if (metric %in% c('pd', 'meanPatristic', 'meanPatristicNN', 'minPatristicNN', 'phyloDisparity', 'phyloWeightedEndemism', 'PSV', 'DR')) {
	 	
	 	# check that there is a phylogeny in speciesRaster object
		if (is.null(x[['phylo']])) {
			stop('speciesRaster object does not contain a phylo object!')
		}
		
		if (verbose) message('\t...dropping species that are not in phylo data...\n')
	 	# prune speciesRaster object down to species shared with phylogeny
		x[['speciesList']] <- intersectList(x[['speciesList']], x[['phylo']]$tip.label)
	
	} else if (!metric %in% c('weightedEndemism', 'correctedWeightedEndemism')) {
		stop('Metric not recognized!')
	} else {
		# set empty entries to NA
		emptyEntries <- which(lengths(x[[2]]) == 0)
		x[['speciesList']][emptyEntries] <- rep(list(NA), length(emptyEntries))
	}
	
	uniqueComm <- x[['speciesList']]

	## ----------------------------------
	## MORPHOLOGY-RELATED METRICS
	
	## UNIVARIATE
	
	if (metric %in% c('mean', 'median', 'variance', 'range', 'mean_NN_dist', 'min_NN_dist') & !pairwise & metricType == 'uniVar') {
		if (verbose) message('\t...calculating univariate metric: ', metric, '...\n')
		trait <- x[['data']]
		resVal <- cellAvg(uniqueComm, trait = trait, stat = metric)
		
		# to be consistent with other metrics, we will have 'variance' return 0 for single-species cells
		if (metric == 'variance') {
			ind <- which(lengths(uniqueComm) == 1 & sapply(uniqueComm, anyNA) == FALSE)
			resVal[ind] <- 0
		}
	}
	
	if (metric %in% c('arithmeticWeightedMean', 'geometricWeightedMean') & metricType == 'uniVar') {
		if (verbose) message('\t...calculating univariate metric: ', metric, '...\n')
		trait <- x[['data']]
			
		# get inverse range sizes (1 / cell counts)
		inverseCellCount <- 1 / x[['cellCount']]
		
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(lengths(uniqueComm) > 0 & sapply(uniqueComm, anyNA) == FALSE)

		if (metric == 'arithmeticWeightedMean') {
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
				(1 / sum(inverseCellCount[y])) * sum(inverseCellCount[y] * trait[y])
			})
		}
		if (metric == 'geometricWeightedMean') {
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
				exp(sum(inverseCellCount[y] * log(trait[y])) / sum(inverseCellCount[y]))
			})			
		}	
	}
	
	if (metric %in% c('mean', 'median', 'variance') & pairwise) {
	
		# if pairwise matrix
		if (verbose) message('\t...calculating pairwise metric: ', metric, '...\n')
		
		resVal <- rep(NA, length(uniqueComm)) # set up with NA
		
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, anyNA) == FALSE)

		if (metric == 'mean') {
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) mean(unlist(x[['data']][y, y]), na.rm = TRUE))
		} else if (metric == 'median') {
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) stats::median(unlist(x[['data']][y, y]), na.rm = TRUE))
		} else if (metric == 'variance') {
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) stats::var(unlist(x[['data']][y, y]), na.rm = TRUE))
			ind <- which(lengths(uniqueComm) == 1 & sapply(uniqueComm, anyNA) == FALSE)
			resVal[ind] <- 0
		}
	}
	
	if (metric == 'phylosignal' & metricType == 'uniVar') {
		if (verbose & metricType == 'uniVar') message('\t...calculating univariate metric: ', metric, '...\n')
		resVal <- rep(NA, length(uniqueComm)) # set up with NA
		ind <- which(lengths(uniqueComm) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) geomorph::physignal(x[['data']][y], phy = ape::keep.tip(x$phylo, y), iter = 999, print.progress = FALSE)$phy.signal)	
	}

	## MULTIVARIATE
	
	if (metric == 'disparity' & metricType == 'multiVar') {
		if (verbose) message('\t...calculating multivariate metric: ', metric, '...\n')
		# sum of the diagonal of the covariance matrix
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		# resVal <- rep(NA, length(uniqueComm)) # set up with NA
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, length) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) sum(diag(cov(x[['data']][y,]))))
	}
	
	if (metric == 'range' & metricType == 'multiVar') {
		if (verbose) message('\t...calculating multivariate metric: ', metric, '...\n')
		# maximum of the distance matrix (0 if one sp)
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, length) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) max(dist(x[['data']][y, ])))
	}
	
	if (metric == 'rangePCA' & metricType == 'multiVar') {
		if (verbose) message('\t...calculating multivariate metric: ', metric, '...\n')
		pc <- prcomp(x[['data']])
		# retain 99% of the variation
		keep <- 1:which(cumsum(((pc$sdev)^2) / sum(pc$sdev^2)) >= 0.99)[1]
		if (length(keep) == 1) { keep <- 1:ncol(pc$x) }
		pc <- pc$x[,keep]
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, length) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
			sum(apply(pc[y,], 2, function(z) diff(range(z))))
		})	
	}
	
	if (metric == 'mean_NN_dist' & metricType == 'multiVar') {
		if (verbose) message('\t...calculating multivariate metric: ', metric, '...\n')
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, length) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) mean(nnDist(x[['data']][y, ], Nrep = nreps)$mean_dist))
	}

	if (metric == 'min_NN_dist' & metricType == 'multiVar') {
		if (verbose) message('\t...calculating multivariate metric: ', metric, '...\n')
		resVal <- numeric(length = length(uniqueComm)) # set up with zeros
		resVal[sapply(uniqueComm, anyNA)] <- NA
		ind <- which(sapply(uniqueComm, length) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) min(nnDist(x[['data']][y, ], Nrep = nreps)$mean_dist))
	}
	
	if (metric == 'phylosignal' & metricType == 'multiVar') {
		if (verbose & metricType == 'multiVar') message('\t...calculating multivariate metric: ', metric, '...\n')
		resVal <- rep(NA, length(uniqueComm)) # set up with NA
		ind <- which(lengths(uniqueComm) > 1)
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) geomorph::physignal(x[['data']][y, ], phy = ape::keep.tip(x$phylo, y), iter = 999, print.progress = FALSE)$phy.signal)	
	}
	
	## ----------------------------------
	## PHYLOGENY-RELATED METRICS
	
	if (metric %in% c('pd', 'meanPatristic', 'minPatristicNN', 'meanPatristicNN', 'phyloDisparity', 'PSV', 'DR')) {
		
		# calculate pairwise patristic distance
		patdist <- cophenetic(x[['phylo']])
		diag(patdist) <- NA
		
		if (metric == 'pd') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			mat <- matrix(0, nrow = length(uniqueComm), ncol = length(x$phylo$tip.label))
			colnames(mat) <- x$phylo$tip.label
			for (i in 1:length(uniqueComm)) {
				mat[i, colnames(mat) %in% uniqueComm[[i]]] <- 1
			}
			resVal <- picante::pd(mat, x$phylo, include.root = TRUE)$PD
			resVal[sapply(uniqueComm, anyNA)] <- NA
		}
		
		if (metric == 'meanPatristic') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# meanPatristic is 0 if 1 species, NA if no species
			patdist[upper.tri(patdist, diag = TRUE)] <- NA
			resVal <- pbapply::pbsapply(uniqueComm, function(y) mean(unlist(patdist[y, y]), na.rm = TRUE))
			resVal[sapply(uniqueComm, length) == 1] <- 0
			resVal[sapply(uniqueComm, anyNA)] <- NA
		}
		
		if (metric == 'meanPatristicNN') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# the mean of the minimum patristic distance for each species present
			resVal <- numeric(length = length(uniqueComm)) # set up with zeros
			resVal[sapply(uniqueComm, anyNA)] <- NA
			ind <- which(sapply(uniqueComm, length) > 1)
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
				return(mean(apply(patdist[y, y], MARGIN = 1, min, na.rm = TRUE)))
			})
		}

		if (metric == 'minPatristicNN') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# the minimum of the minimum patristic distance for each species present
			resVal <- numeric(length = length(uniqueComm)) # set up with zeros
			resVal[sapply(uniqueComm, anyNA)] <- NA
			ind <- which(sapply(uniqueComm, length) > 1)
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
				return(min(apply(patdist[y, y], MARGIN = 1, min, na.rm = TRUE)))
			})
		}

		if (metric == 'phyloDisparity') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# the sum of the squared deviations from the mean
			# value of 0 if 1 species, NA if no species
			
			resVal <- numeric(length = length(uniqueComm)) # set up with zeros
			resVal[sapply(uniqueComm, anyNA)] <- NA
			ind <- which(sapply(uniqueComm, length) > 1)

			patdist[upper.tri(patdist, diag = TRUE)] <- NA
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) sum((patdist[y, y] - mean(patdist[y,y], na.rm = TRUE)) ^ 2, na.rm = TRUE) / choose(length(y), 2))
		}		
		
		if (metric == 'PSV') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# Phylogenetic Species Variability is based on the variance-covariance matrix
			# measure of phylogenetic diversity, ranges from 0-1, not confounded by species richness
			
			vmat <- ape::vcv.phylo(x[['phylo']], corr = TRUE)
			
			resVal <- numeric(length = length(uniqueComm)) # set up with zeros
			resVal[sapply(uniqueComm, anyNA)] <- NA
			ind <- which(sapply(uniqueComm, length) > 1)
			
			resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) (length(y) * sum(diag(vmat[y,y])) - sum(vmat[y,y]))/(length(y) * (length(y) - 1)))
		}
		
		if (metric == 'DR') {
			if (verbose) message('\t...calculating phylo metric: ', metric, '...\n')
			# DR statistic for speciation rates
			
			tipDR <- DRstat(x[['phylo']])			
			resVal <- pbapply::pbsapply(uniqueComm, function(y) mean(tipDR[y]))
		}
	}
	
	## ----------------------------------
	## RANGE-WEIGHTED METRICS
	
	if (metric %in% c('weightedEndemism', 'correctedWeightedEndemism')) {
		if (verbose) message('\t...calculating weighted endemism metric...\n')
		
		# get inverse range sizes (1 / cell counts)
		inverseCellCount <- 1 / x[['cellCount']]

		resVal <- rep(NA, length(uniqueComm)) # set up with NA
		ind <- which(sapply(uniqueComm, anyNA) == FALSE)		
		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) sum(inverseCellCount[y]))
	
		# if corrected metric, we will standardize by species richness
		if (metric == 'correctedWeightedEndemism') {
			resVal[ind] <- resVal[ind] / lengths(uniqueComm[ind])
		}
	}
	
	if (metric == 'phyloWeightedEndemism') {
		if (verbose) message('\t...calculating phylogenetic weighted endemism...\n')
		spEdges <- getRootToTipEdges(x[['phylo']])
		if (!'edgeArea' %in% names(x)) {
			if (verbose) message('\t...calculating branch-specific range sizes...\n')

			# phylo branch ranges must be based on full list of cell communities, therefore we need to expand the speciesList
			fullSpList <- expandSpeciesCellList(x)
			x[['edgeArea']] <- do.call(cbind, phyloBranchRanges(x[['phylo']], convertNAtoEmpty(fullSpList), spEdges))
		}
		tipIndVec <- sapply(x[['phylo']]$tip.label, function(z) which(x[['phylo']]$tip.label == z))
		resVal <- rep(NA, length(uniqueComm)) # set up with NA
		ind <- which(sapply(uniqueComm, anyNA) == FALSE)

		resVal[ind] <- pbapply::pbsapply(uniqueComm[ind], function(y) {
			commEdges <- unique(unlist(spEdges[tipIndVec[y]])) + 1
			sub <- x[['edgeArea']][commEdges,]
			if (inherits(sub, 'numeric')) {
				sub <- matrix(sub, nrow = 1)
			}
			sum(sub[,1] / sub[,2])
		})
	}
	
 	
	
	## ----------------------------------
	## REMAP RESULTS TO RELEVANT CELLS AND ASSIGN TO RASTERS
	
	resVal[is.nan(resVal)] <- NA
	cellVec <- numeric(length = length(x[['cellCommInd']]))
	for (i in 1:length(resVal)) {
		cellVec[x[['cellCommInd']] == i] <- resVal[i]
	}
	resRas <- raster::raster(x[[1]])
	raster::values(resRas) <- cellVec
	names(resRas) <- metric
	
	# prepare output object
	obj <- x
	obj[[1]] <- resRas
	obj[['geogSpecies']] <- sort(unique(unlist(x[['geogSpecies']])))
	
	class(obj) <- 'speciesRaster'
	return(obj)	
}




# # function to calculate phylo branch specific range areas

# phyloBranchRanges <- function(x) {

	# if (!'speciesRaster' %in% class(x)) {
		# stop('x must be of class speciesRaster.')
	# }
	
	# # check that there is a phylogeny in speciesRaster object
	# if (is.null(x[['phylo']])) {
		# stop('speciesRaster object does not contain a phylo object!')
	# }
	
	# phylo <- x[['phylo']]
	
	# allLeaves <- phangorn::Descendants(phylo, phylo$edge[,2], type = 'tips')
	# allLeaves <- lapply(allLeaves, function(y) phylo$tip.label[y])
	
	# branchTable <- matrix(nrow = length(phylo$edge.length), ncol = 2)
	# colnames(branchTable) <- c('branchLength', 'branchArea')
	# branchTable[,1] <- phylo$edge.length
	
	# for (i in 1:length(phylo$edge.length)) {
		
		# inCell <- sapply(x[['speciesList']], function(y) any(allLeaves[[i]] %in% y))
		# branchTable[i, 2] <- sum(inCell)	
	# }

	# x[['edgeArea']] <- branchTable
	# return(x)	
# }



# getCommEdges2 <- function(phylo, comm, indList) {
	
	# tipInd <- sapply(comm, function(z) which(phylo$tip.label == z))
	# edges <- unique(unlist(indList[tipInd]))
	# edges <- edges + 1
	# # edgeColors <- ifelse(1:length(phylo$edge.length) %in% edges, 'red','black')
	# # edgeWidths <- ifelse(1:length(phylo$edge.length) %in% edges, 2, 1)
	# # plot(phylo, tip.color=ifelse(phylo$tip.label %in% comm, 'blue','gray'), root.edge=TRUE, edge.color = edgeColors, edge.width = edgeWidths)
	# # nodelabels(node=nodes, frame='circle', bg='red', cex=0.1)
	
	# return(edges)
# }




# # library(raster)
# library(maptools)
# library(Rcpp)
# # convert polygon ranges to raster
# ranges <- rasterStackFromPolyList(tamiasPolyList, resolution = 20000)
# x <- createSpeciesRaster(ranges = ranges)
# x <- addPhylo_speciesRaster(x, tamiasTree)


convertNAtoEmpty <- function(spCellList) {
	ind <- which(sapply(spCellList, anyNA) == TRUE)
	spCellList[ind] <- rep(list('empty'), length(ind))
	return(spCellList)
}



# # phylo2 <- x[['phylo']]
# class(phylo2) <- 'list'

# sourceCpp('~/Desktop/testing.cpp')

# tipEdges <- getRootToTipEdges(phylo2)

# q <- phyloBranchRanges(phylo2, spCellList, tipEdges)
# q <- do.call(cbind, q)
ptitle/speciesRaster documentation built on June 11, 2022, 1:35 a.m.