R/plot.R

In AlexCast/lisst: Functions for manipulation of the Laser In-Situ Scattering and Transmissometry (LISST) data

#' Auxiliary plotting functions:

.xaxn <- function(x, by) {
	switch(by,
		sample = list("Sample Number", 1:nrow(x)),
		depth  = list(make_unit_label("'Depth'", x$Depth), x$Depth),
		time   = list("Time", ltime(x))
	)
}

.yaxn <- function(x) {
	lmodl <- attr(x, "lmodl")
	ity   <- attr(x, "lproc")$ity
	switch(attr(x, "type"),
		raw = ,
		cor = ,
		cal = list("'Bins'", set_units(1:lmodl$nring, 1)),
		vsf = list("'Scattering Angle'", lmodl$wang[, 3]),
		vol = ,
		pnc = ,
		csa = list("'Size'", lmodl$binr[[ity]][, 3])
	)
}

.zaxn <- function(x) {
	switch(attr(x, "type"),
		raw = ,
		cor = "'Digital Number'",
		cal = "'Laser Power'",
		vsf = "'Volume Scattering Function'",
		vol = "'Volume Concentration'",
		pnc = "'Number Concentration'",
		csa = "'Cross sectional Concentration'"
	)
}

.map2color <- function(x, pal, limits){
    if(missing(limits)) limits <- range(x, na.rm = T)
    pal[findInterval(x, seq(limits[1], limits[2], length.out = length(pal)+1), all.inside = TRUE)]
}



#' Hovmöller diagram for lisst objects
#'
#' Create Hovmöller diagram for lisst objects based on sample number, 
#' measurement time or depth.
#'
#' @param x      A lisst object.
#' @param by     Ordinate axis dimension. Must be one of 'sample', 'time' or 
#'               'depth'.
#' @param norm   Logic. Should the magnitude values de normalized? See details.
#' @param legend Logic. Should a legend be added to the plot?
#' @param total  Logic. Should an upper pannel with the total magnitude be 
#'               added?
#' @param col    A color pallete to be used. If not specified a default pallete 
#'               will be used.
#' @param yu     The units for the abscissa. See details.
#' @param zu     The units for the magnitude values. See details.
#' @param brks   Passed to lstat function. See ?lstat for details.
#'
#' @details
#' The abscissa of the Hovmöller diagram will show the ring number if data is 
#' 'raw' or 'cor' and the median angle or size of the bin, as appropriate for 
#' the lisst object type.
#' 
#' The z scale is provided in log10 for appropriate color mapping. However, 
#' LISST SOP will return 0 volume concentration if concentration is below 0.001
#' ppm. In that case, zeros as treated as half this minimum value, 0.0005 ppm 
#' for ploting purposes only.
#'
#' If the lisst object is not regularly spaced in the chosen dimension, it will 
#' be binned for raster representation. The nbins parameter will control the 
#' number of bins, with all bins having the same interval size. The sample 
#' dimension will always be regular.
#'
#' Normalization is performed by the sum of all detector values. It is provided 
#' to aid visual separation from changes in concentration and distribution. If 
#' TRUE a top pannel is added to the plot, with a line plot indicating the 
#' magnitude of the sum of detector values.
#'
#' If legend is set to TRUE a bottom pannel will be added to the plot with a 
#' magnitude scale. The axis label and units are affected by the zlab and zu 
#' parameters.
#'
#' Note that unit conversion is performed 'automatcaly' using the units package. 
#' If there is no standard conversion between the default and the request units
#' the function will exit with an error. See \code{link{units}} for details.
#'
#' @examples
#' flp <- system.file("extdata", "DN_27_rs.asc", package = "lisst")
#' model <- "100CX" 
#' lop <- read_lisst(flp, model = model)
#' lhov(lop, by = 'sample')
#' lhov(lop, by = 'sample', norm = FALSE)
#' lhov(lop, by = 'sample', norm = FALSE, legend = FALSE)
#'
#' @export

lhov <- function(x, by = 'sample', norm = TRUE, legend = TRUE, total = TRUE, col, 
        xlab, ylab, zlab, yu, zu, brks) {

	lmodl  <- attr(x, "lmodl")
	lty    <- attr(x, "type")

	# Get axis, labels and set units:
	xaxn <- .xaxn(x, by)
	if(missing(xlab)) xlab <- xaxn[[1]]
	yaxn <- .yaxn(x)
	if(!missing(yu)) units(yaxn[[2]]) <- yu
	if(missing(ylab)) ylab <- units::make_unit_label(yaxn[[1]], yaxn[[2]])
	zaxn <- .zaxn(x)
	if(missing(zlab)) {
		if(norm) zlab <- paste("Normalized", parse(text = zaxn))
		else zlab <- units::make_unit_label(zaxn, x[, 1])
	}

	# Check ordinates:
	if(by == 'time' && any(diff(xaxn[[2]]) < 0))
		stop('Time index must be monotonicaly increasing', call. = F)
	if(by == 'depth') {
		if(diff(range(diff(drop_units(xaxn[[2]])))) != 0) {
			warning('Bining irregular data', call. = F)
                        if(missing(brks)) {
				hx <- hist(drop_units(xaxn[[2]]), plot = F)
				bx <- hx$breaks
				bx <- c(bx[1], rep(bx[-c(1, length(bx))], each = 2), bx[length(bx)])
				dim(bx) <- c(2, length(bx) / 2)
				brks <- list()
				for(i in 1:ncol(bx)) brks[[i]] <- paste(bx[, i], collapse = "|")
			} else {
                                hx <- list()
				hx$mids <- apply(sapply(sapply(brks, strsplit, split = "|", fixed = T), as.numeric), 2, mean)
                        }
			x <- lstat(x, brks = brks, 'mean')
			units(hx$mids) <- units(xaxn[[2]])
			xaxn[[2]] <- hx$mids
		}
	}

	# Deal with possible zeros in 'vol' and 'pnc'. 
	# The LISST SOP will give a 0 if concentration is below 0.001 ppm.
#	if(lty == 'vol' | lty == 'pnc') {
#		x <- lget(x, 'vol')
#		for(i in 1:lmodl$nring)
#			x[, i][which(drop_units(x[, i]) == 0)] <- 0.0005
#		x <- lget(x, lty)				
#	}

	xm     <- do.call(rbind, x[, 1:lmodl$nring])
	units(xm) <- units(x[, 1])
	if(!missing(zu)) xm <- set_units(xm, zu)

	# Set device division:
	if(legend && total) {
		layout(matrix(c(1, 2, 3), ncol = 1), heights = c(1, 4, 1))
	} else if(legend) {
		layout(matrix(c(1, 2), ncol = 1), heights = c(4, 1))
	} else if(total) {
		layout(matrix(c(1, 2), ncol = 1), heights = c(1, 4))
	}

	if(missing(col)) {
		col <- colorRampPalette(c("#011F4B", "#03396C", "#005B96", "#6497B1", "#B3CDE0"))(256)
	}

	xs <- xm[1, ]
	for(i in 1:ncol(xm)) xs[i] <- sum(xm[, i], na.rm = T)
	xs[drop_units(xs) == 0] <- NA
	if(norm) {
		xm <- xm / rep(xs, each = nrow(xm))
	}

	if(total) {
		par(mar = c(0, 5.5, 0.5, 1), las = 1)
		nylab <- units::make_unit_label('Total', x[, 1])
		plot(xaxn[[2]], drop_units(xs), axes = F, xlab = "", type = "l", xaxs = 'i', 
			ylab = '', log = 'y')
		mtext(nylab, side = 2, line = 3.5, las = 0, cex = 0.7)
		axis(2, crt = 90)
	}

	par(mar = c(5, 5.5, 2, 1))
	plot(NA, xlim = c(1, ncol(xm) + 1), ylim = c(1, lmodl$nring + 1), bty = "l", xaxs = 'i', 
		yaxs = 'i', xaxt = 'n', yaxt = 'n', xlab = xlab, ylab = '')
	mtext(ylab, side = 2, line = 3.5, las = 0, cex = 0.7)
	id <- c(seq(1, lmodl$nring, by = 4), lmodl$nring)
	axis(2, at = id+0.5, labels = round(yaxn[[2]][id], 4))
	if(by == 'time') {
		axis(1, at = axTicks(1) + 0.5, labels = format(xaxn[[2]][axTicks(1)], "%Y-%m-%d\n%H:%M:%S"), padj = 0.5)
	} else {
		axis(1, at = axTicks(1) + 0.5, labels = xaxn[[2]][axTicks(1)])
	}
	xm   <- log10(drop_units(xm))
	xm[abs(xm) == Inf] <- NA
	xrst <- matrix(.map2color(xm, col), ncol = ncol(xm))[nrow(xm):1, ]
	xrst[is.na(xrst)] <- '#000000' 
	rasterImage(xrst, xleft = 1, ybottom = 1, xright = ncol(xm) + 1, ytop = lmodl$nring + 1, 
		interpolate = F)

	if(legend) {
		par(mar = c(5, 8, 1, 4))
		plot(NA, xlim = range(xm, na.rm = T), ylim = c(0, 1), xaxs = 'i', yaxs = 'i', xaxt = 'n', 
			yaxt = 'n', xlab = zlab, ylab = "")
		digits <- as.numeric(sapply(strsplit(formatC(10^axTicks(1), format = 'e'), "e"), 
			'[[', 2))
		digits[digits > 0] <- 0
		axis(1, at = axTicks(1), labels = round(10^axTicks(1), abs(digits)))
		xrst <- t(as.raster(.map2color(seq(min(xm, na.rm = T), max(xm, na.rm = T), length.out = length(col)), col)))
		rasterImage(xrst, xleft = min(xm, na.rm = T), ybottom = 0, xright = max(xm, na.rm = T), ytop = 1, interpolate = T)
	}
}

#' Plot Bins
#'
#' @param x    A lisst object.
#' @param bins A list of the rings to be summed. See details and examples.
#' @param by   Ordinate. One of 'sample', 'time' or 'depth'
#' @param col  A vector of colors.
#' @param lty  A vector of line types.
#'
#' @export

plotBins <- function(x, bins, by, col, lty) {
	stopifnot(is.lisst(x))
	if(!is(bins, 'list')) stop('bins must be a list', call. = FALSE)

	lmodl <- attr(x, 'lmodl')
	xm    <- as.matrix(x[, 1:lmodl$nring])
	xm[xm == 0] <- 0.005
	units(xm) <- units(x[, 1])
	x[, 1:lmodl$nring] <- as.data.frame(xm)

	xsm  <- list()
	for(i in 1:length(bins)) xsm[[i]] <- lrings(x, bins[[i]])
	xsm  <- do.call(rbind, xsm)
	xsm[xsm == 0] <- NA
	ylab <- units::make_unit_label(.zaxn(x), x[, 1])	
	ylim <- range(xsm, na.rm = T)
	xaxn <- .xaxn(x, by)

	if(missing(col)) {
		col <- rainbow(length(bins), start = 0.1)
	} else if(length(col) != length(bins)) {
		col <- rep(col, length(bins))[1:length(bins)] 
	}
	if(missing(lty)) {
		lty <- rep(1, length(bins))
	} else if(length(lty) != length(bins)) {
		lty <- rep(lty, length(bins))[1:length(bins)] 
	}

	plot(xaxn[[2]], xsm[1, ], col = col[1], xlab = xaxn[[1]], ylab = ylab,
		ylim = ylim, log = 'y', type = 'l', lty = lty[1])
	for(i in 1:nrow(xsm)) lines(xaxn[[2]], xsm[i, ], col = col[i], lty = lty[i])
	
	invisible(xsm)
}



#l200ps1 <- lrings(l200p, bins = 1:12)
#l200ps2 <- lrings(l200p, bins = 25:36)

#ylim <- range(c(l200ps1, l200ps2), na.rm = T)
#par(mar = c(5, 5, 1, 1))
#plot(ltime(l200p), drop_units(l200ps1), type = 'l', xlab = 'Time', ylab = ylab, ylim = ylim, col = 'blue', log = 'y')
#lines(ltime(l200p), l200ps2, col = 'red')
#legend('topleft', c('1 to 10 µm', '66 to 500 µm'), lty = 1, col = c("blue", 'red'), bty = 'n')

#}




# if xu is not defined the default unit will be used.

plot.lisst <- function(lo, xu, type = "rings", ...) {

	linst <- attr(lo, "linst")
	lmodl <- attr(lo, "lmodl")
	typ   <- attr(lo, "type")

	if(missing(type))
		type <- attr(lo, "type")

	if(type == 'rings') {
		
	}


	if(type == "vsf" || type == "pf") {
		if(attr(lo, "type") != "vsf")
			stop("vsf or pf plot requires a vsf type lisst object", call. = FALSE)
		if(!missing(xu) && !(xu == "degree" || xu == "radian"))
			stop("xu for vsf plots must be 'degree' or 'rad'", call. = FALSE)
		x <- lmodl$wang[, 3]
		if(!missing(xu)) units(x) <- xu
		id <- which(lo[, 1] == set_units(0, 1/m/sr))
		y <- set_units(as.matrix(lo[, 1:lmodl$nring]), 1/m/sr)
		if(length(id) > 0) y <- y[-id, ]
		if(type == "pf") {
			cp <- lo[, "Beam attenuation"]
			aw670 <- units::set_units(0.439, 1/m)
			bw670 <- units::set_units(0.0005808404, 1/m)
			cp <- cp - ((aw670 + bw670) * as.numeric(lmodl$pl))
			if(length(id) > 0) cp <- cp[-id, ]
			y <- y / cp
			plot.vsf(x, y, pf = TRUE, ...)
		} else {
			plot.vsf(x, y, pf = FALSE, ...)
		}
	}
	if(type == 'vol' || type == 'pnc' || type == 'psd') {
		x <- lmodl$binr[[attr(lo, "lproc")$ity]]
		if(!missing(xu)) units(x) <- xu
		id <- which(as.numeric(lo[, 1]) == 0)
		if((type == 'vol' && typ == 'pnc') || type == 'psd') lo <- lgetvol(lo)
		if(type == 'pnc' && typ == 'vol') lo <- lgetpnc(lo)
		y <- as.matrix(lo[, 1:lmodl$nring, drop = FALSE])
		units(y) <- units(lo[, 1])
		y2 <- NULL
		if(type == 'psd' || type == 'pnc') {
			lo <- lgetpnc(lo)
			y2 <- as.matrix(lo[, 1:lmodl$nring, drop = FALSE])
			units(y2) <- units(lo[, 1])
		}
		if(length(id) > 0) {
			y <- y[-id, ]
			if(!is.null(y2)) y2 <- y2[-id, ]
		}
		plot.psd(x, y, y2, type, ...)

	}
}

plot.psd <- function(x, y, y2, type, ...) {
	par(mar = c(5, 5, 4, 5))
	yn <- "'Particle Volume Concentration'"
	ylab <- units::make_unit_label(yn, y)
	if(type == 'pnc') {
		yn <- "'Particle Number Concentration'"
		ylab <- units::make_unit_label(yn, y2)
	}
	xlab <- units::make_unit_label("'Particle size'", x)
	ylim <- range(y, na.rm = T)
	if(type == 'pnc') ylim <- range(y2, na.rm = T)
	plot(NA, xlim = range(c(x[1, 1], x[, 2])), ylim = ylim, log = "xy", 
	xlab = xlab, ylab = ylab, bty = "l", yaxs = "i", ...)
	if(type == 'psd' || type == 'vol') {
		for(i in 1:ncol(y)) {
			rect(xleft = x[i, 1], ybottom = ylim[1], xright = x[i, 2], 
			ytop = y[1, i])
		}
	} else {
		for(i in 1:nrow(y2)) lines(x[, 3], y2[i, ], lty = 1)
	}
	if(type == 'psd') {
		par(new = T)
		plot(NA, xlim = range(c(x[1, 1], x[, 2])), ylim = range(y2), 
		log = "xy", xlab = "", ylab = "", bty = "n", main = "", xaxt = "n", 
		yaxt = "n")
		bor <- par("usr")
		lines(c(10^bor[2], 10^bor[2]), c(10^bor[3], 10^bor[4]), col = "red")
		axis(4, col = "red", col.axis = "red")
		lines(x[, 3], y2[1, ], lty = 1, col = "red")
		yn <- "'Particle Number Concentration'"
		ylab <- units::make_unit_label(yn, y2)
		mtext(ylab, side = 4, line = 3, col = "red")
	}
}

plot.vsf <- function(x, y, pf = FALSE, ...) {
	par(mar = c(5, 5, 3, 2))
	ylim <- range(y, na.rm = T)
	yn <- "'Particle Volume Scattering Function'"
	if(pf) yn <- "'Particle Phase Function'"
	ylab <- units::make_unit_label(yn, y)
	xlab <- units::make_unit_label("'In water scattering angle'", x) 
	plot(x, y[1, ], type = "l", log = "xy", ylim = ylim, xlab = xlab, ylab = ylab, ...)
	if(nrow(y) > 1)
		for(i in 1:nrow(y))
			lines(x, y[i, ])
}