Nothing
R/plot.R
In argoFloats: Analysis of Oceanographic Argo Floats
# vim:textwidth=128:expandtab:shiftwidth=4:softtabstop=4
colDefaults <- list(core="7", bgc="#05f076", deep="6")
#' Plot an argoFloats Object
#'
#' The action depends on the `type` of the object, and
#' this is set up by the function that created the object;
#' see \dQuote{Details}. These are basic plot styles, with
#' somewhat limited scope for customization. Since the data with
#' [argoFloats-class] objects are easy to extract, users should
#' not find it difficult to create their own plots to meet a
#' particular aesthetic. See \dQuote{Examples} and Kelley et al. (2021)
#' for more plotting examples.
#'
#' The various plot types are as follows.
#'
#' * For `which="map"`, a map of profile locations is created if subtype
#' is equal to cycles, or a rectangle highlighting the trajectory of a
#' float ID is created when subtype is equal to trajectories. This
#' only works if the `type` is `"index"` (meaning that `x` was created
#' by [getIndex()] or a subset of such an object, created with
#' [subset,argoFloats-method()]), or `argos` (meaning that
#' `x` was created with [readProfiles()]. The plot range is
#' auto-selected. If the `ocedata` package
#' is available, then its `coastlineWorldFine` dataset is used to draw
#' a coastline (which will be visible only if the plot region
#' is large enough); otherwise, if the `oce` package is available, then its
#' `coastlineWorld` dataset is used.
#' The `bathymetry` argument controls whether (and how) to draw a map underlay
#' that shows water depth. There are three possible values for `bathymetry`:
#' 1. `FALSE` (the default, via [argoDefaultBathymetry()]), meaning not to draw bathymetry;
#' 2. `TRUE`, meaning to draw bathymetry as an
#' image, using data downloaded with [oce::download.topo()].
#' Example 4 illustrates this, and also shows how to adjust
#' the margins after the plot, in case there is a need to add
#' extra graphical elements using [points()], [lines()], etc.
#' 3. A list with items controlling both the bathymetry data and its
#' representation in the plot (see Example 4). Those items are:
#'
#' 1. `source`, a mandatory value that is either
#' (a) the string `"auto"` (the default) to use
#' [oce::download.topo()] to download the data
#' or (b) a value returned by [oce::read.topo()].
#' 2. `contour`, an optional logical value (with `FALSE` as the default)
#' indicating whether to represent
#' bathymetry with contours (with depths of 100m, 200m, 500m shown,
#' along with 1km, 2km up to 10km), as opposed to an image;
#' 3. `colormap`, ignored if `contour` is `TRUE`,
#' an optional value that is either the string `"auto"` (the default)
#' for a form of GEBCO colors computed with
#' [oce::oceColorsGebco()], or a value computed with [oce::colormap()]
#' applied to the bathymetry data; and
#' 4. `palette`, ignored if `contour` is `TRUE`,
#' an optional logical value (with `TRUE` as the default)
#' indicating whether to draw a depth-color
#' palette to the right of the plot.
#' Note that the default value for `bathymetry` is set by a call
#' to [argoDefaultBathymetry()], and that this second function can only handle possibilities
#' 1 and 2 above.
#'
#' * For `which="profile"`, a profile plot is created, showing the variation of
#' some quantity with pressure or potential density anomaly, as specified by the
#' `profileControl` argument.
#'
#' * For `which="QC"`, two time-series panels are shown, with
#' time being that recorded in the individual profile in the dataset.
#' An additional argument named `parameter` must be given, to name the
#' quantity of interest. The function only works if `x` is an
#' [`argoFloats-class`] object created with [readProfiles()].
#' The top panel shows the percent of data flagged with codes
#' 1 (meaning good data), 2 (probably good), 5 (changed)
#' or 8 (estimated), as a function of time (lower axis) and
#' (if all cycles are from a single Argo float)
#' cycle number (upper axis, with smaller font).
#' Thus, low values on the top panel reveal
#' profiles that are questionable. Note that if all of data at a given time
#' have flag 0, meaning not assessed, then a quality of 0 is plotted at that
#' time. The bottom panel shows the mean value of the parameter in question
#' regardless of the flag value.
#'
#' * For `which="summary"`, one or more time-series panels are shown
#' in a vertical stack. If there is only one ID in `x`, then the cycle
#' values are indicated along the top axis of the top panel. The choice
#' of panels is set by the `summaryControl` argument.
#'
#' * For `which="TS"`, a TS plot is created, by calling [plotArgoTS()]
#' with the specified `x`, `xlim`, `ylim`, `xlab`, `ylab`,
#' `type`, `cex`, `col`, `pch`, `bg`, `eos`, and `TSControl`,
#' along with `debug-1`. See the help for [plotArgoTS()]
#' for how these parameters are interpreted.
#'
#' @param x an [`argoFloats-class`] object.
#'
#' @param which a character value indicating the type of plot. The possible
#' choices are `"map"`, `"profile"`, `"QC"`, `"summary"` and `"TS"`;
#' see \dQuote{Details}.
#'
#' @param bathymetry an argument used only if `which="map"`, to control
#' whether (and how) to indicate water depth. Note that the default
#' was `TRUE` until 2021-12-02, but was changed to `FALSE` on that
#' date, to avoid a bathymetry download, which can be a slow operations.
#' See \dQuote{Details} for details,
#' and Example 4 for a hint on compensating for the margin
#' adjustment done if an image is used to show bathymetry.
#'
#' @param geographical flag indicating the style of axes
#' for the `which="map"` case, but only if no projection is called
#' for in the `mapControl` argument. With `geographical=0` (which
#' is the default), the axis ticks are labeled with signed
#' longitudes and latitudes, measured in degrees. The signs are dropped
#' with `geographical=1`.
#' In the `geographical=4` case, signs are also dropped, but hemispheres
#' are indicated by writing `S`, `N`, `W` or `E` after axis tick labels,
#' except at the equator and prime meridian.
#' Note that this scheme mimics that used by [oce::plot,coastline-method()].
#'
#' @param xlim,ylim numerical values, each a two-element vector, that
#' set the `x` and `y` limits of plot axes, as for [plot.default()] and
#' other conventional plotting functions.
#'
#' @param xlab a character value indicating the name for the horizontal axis, or
#' `NULL`, which indicates that this function should choose an appropriate name
#' depending on the value of `which`. Note that `xlab` is not obeyed if
#' `which="TS"`, because altering that label can be confusing to the user.
#'
#' @param ylab as `xlab`, but for the vertical axis.
#'
#' @param type a character value that controls the line type, with `"p"` for
#' unconnected points, `"l"` for line segments between undrawn points, etc.;
#' see the docs for [par()], If `type` not specified, it defaults to `"p"`.
#'
#' @param cex a character expansion factor for plot symbols, or `NULL`, to get an
#' value that depends on the value of `which`.
#'
#' @param col the colour to be used for plot symbols, or `NULL`, to get an value
#' that depends on the value of `which` (see \dQuote{Details}). If `which="TS"`,
#' then the `TSControl` argument takes precedence over `col`.
#'
#' @param bg the colour to be used for plot symbol interior, for `pch`
#' values that distinguish between the interior of the symbol and the
#' border, e.g. for `pch=21`.
#'
#' @param pch an integer or character value indicating the type of plot symbol,
#' or `NULL`, to get a value that depends on the value of `which`.
#' (See [par()] for more on specifying `pch`.)
#'
#' @param eos a character value indicating the equation of state to use
#' if `which="TS"`. This must be `"gsw"` (the default) or `"unesco"`;
#' see [oce::plotTS()].
#'
#' @param mapControl a list that permits particular control of the `which="map"`
#' case. If provided, it may contain elements named `bathymetry` (which
#' has the same effect as the parameter `bathymetry`), `colLand`
#' (which indicates the colour of the land), and `projection` (which
#' may be `FALSE`, meaning to plot longitude and latitude on rectilinear axes,
#' `TRUE`, meaning to plot with [oce::mapPlot()], using Mollweide projection that
#' is suitable mainly for world-scale views, or a character value that will be
#' supplied to [oce::mapPlot()]. If a projection is used, then the positions
#' of the Argo floats are plotted with [oce::mapPoints()], rather than with
#' [points()], and if the user wishes to locate points with mouse clicks,
#' then [oce::mapLocator()] must be used instead of [locator()]. If `bathymetry`
#' is not contained in `mapControl`, it defaults to `FALSE`, and if `projection`
#' is not supplied, it defaults to `FALSE`. Note that `mapControl` takes
#' precedence over the `bathymetry` argument, if both are provided.
#' Also note that, at present, bathymetry cannot be shown with map projections.
#'
#' @param profileControl a list that permits control of the `which="profile"`
#' case. If provided, it may contain elements named `parameter` (a character value
#' naming the quantity to plot on the x axis), `ytype` (a character value equal to
#' either `"pressure"` or `"sigma0"`) and `connect` (a logical value indicating
#' whether to skip across `NA` values if the `type` argument is `"l"`, `"o"`,
#' or `"b"`).
#' If `profileControl` is not provided, it defaults to
#' `list(parameter="temperature", ytype="pressure", connect=FALSE)`. Alternatively,
#' if `profileControl` is missing any of the three elements, then they are
#' given defaults as in the previous sentence.
#'
#' @param QCControl a list that permits control of the `which="QC"`
#' case. If provided, it may contain an element named `parameter`, a character
#' value naming the quantity for which the quality-control information is
#' to be plotted, and an element named `dataStateIndicator`, a logical
#' value controlling whether to add a panel showing this quantity
#' (see Reference Table 6 of reference 1 to learn more
#' about what is encoded in `dataStateIndicator`).
#' If not provided, `QCControl` defaults to
#' `list(parameter="temperature",dataStateIndicator=FALSE)`.
#'
#' @param summaryControl a list that permits control of the `which="summary"`.
#' If provided, it should contain an element named `items`, a character vector
#' naming the items to be shown. The possible entries in this vector are
#' `"dataStateIndicator"` (see Reference Table 6 of reference 1,
#' for more information on this quantity), `"length"` (the number of levels
#' in the profile), `"deepest"` (the highest pressure recorded),
#' `"longitude"` and `"latitude"`. If `summaryControl` is not provided,
#' all of these will be shown. If all the elements of `x` have the same
#' `ID`, then the top panel will have ticks on its top axis, indicating
#' the `cycle`.
#'
#' @param TSControl a list that permits control of the `which="TS"`
#' case, and is ignored for the other cases. If `TSControl` is
#' supplied, then other arguments that normally control plots
#' are ignored. `TSControl` may have elements named
#' `drawByCycle`, `cex`, `col`, `pch`, `lty, `lwd`
#' and `type`. The first of these is a logical value indicating
#' whether to plot individual cycles with different aesthetic
#' characteristics, and the other elements describe those characteristics
#' in the usual way. By default, `type` is `"p"`, and the point
#' characteristics are given by the other parameters, which are passed to
#' [rep()] to have one value per cycle. For example,
#' `TSControl=list(drawByCycle=TRUE,col=1:2,type="l")` indicates
#' to draw lines for the individual cycles, with the first being
#' black, the second red, the third black, etc.
#' **FIXME: THIS IS NOT CODED YET.**
#'
#' @param debug an integer specifying the level of debugging.
#'
#' @param \dots extra arguments passed to the plot calls that are made
#' within this function.
#'
#' @examples
#' # Example 1: map profiles in index
#' library(argoFloats)
#' data(index)
#' plot(index)
#'
#' # Example 2: as Example 1, but narrow the margins and highlight floats
#' # within a circular region of diameter 100 km.
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0))
#' plot(index)
#' lon <- index[["longitude"]]
#' lat <- index[["latitude"]]
#' near <- oce::geodDist(lon, lat, -77.06, 26.54) < 100
#' R <- subset(index, near)
#' points(R[["longitude"]], R[["latitude"]],
#' cex=0.6, pch=20, col="red")
#' par(oldpar)
#'
#' # Example 3: TS of a built-in data file.
#' f <- system.file("extdata", "SR2902204_131.nc", package="argoFloats")
#' a <- readProfiles(f)
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(3.3, 3.3, 1, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
#' plot(a, which="TS")
#' par(oldpar)
#'\donttest{
#' # Example 4: map with bathymetry. Note that par(mar) is adjusted
#' # for the bathymetry palette, so it must be adjusted again after
#' # the plot(), in order to add new plot elements.
#' # This example specifies a coarse bathymetry dataset that is provided
#' # by the 'oce' package. In typical applications, the user will use
#' # a finer-scale dataset, either by using bathymetry=TRUE (which
#' # downloads a file appropriate for the plot view), or by using
#' # an already-downloaded file.
#' data(topoWorld, package="oce")
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 2), mgp=c(2, 0.7, 0)) # narrow margins for a map
#' plot(index, bathymetry=list(source=topoWorld))
#' # For bathymetry plots that use images, plot() temporarily
#' # adds 2.75 to par("mar")[4] so the same must be done, in order
#' # to correctly position additional points (shown as black rings).
#' par(mar=par("mar") + c(0, 0, 0, 2.75))
#' points(index[["longitude"]], index[["latitude"]],
#' cex=0.6, pch=20, col="red")
#' par(oldpar)
#'
#' # Example 5. Simple contour version, using coarse dataset (ok on basin-scale).
#' # Hint: use oce::download.topo() to download high-resolution datasets to
#' # use instead of topoWorld.
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 1))
#' data(topoWorld, package="oce")
#' plot(index, bathymetry=list(source=topoWorld, contour=TRUE))
#' par(oldpar)
##
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 5B. World view with Mollweide projection (Canada Day, 2020)
## jul1 <- subset(getIndex(), time=list(from="2020-09-01", to="2020-09-02"))
## par(mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0)) # narrow margins for a map
## plot(jul1, which="map", mapControl=list(projection=TRUE), bathymetry=FALSE,
## pch=20, col=4, cex=0.75)
#'
#' # Example 6. Customized map, sidestepping plot,argoFloats-method().
#' lon <- topoWorld[["longitude"]]
#' lat <- topoWorld[["latitude"]]
#' asp <- 1/cos(pi/180*mean(lat))
#' # Limit plot region to float region.
#' xlim <- range(index[["longitude"]])
#' ylim <- range(index[["latitude"]])
#' # Colourize 1km, 2km, etc, isobaths.
#' contour(x=lon, y=lat, z=topoWorld[["z"]], xlab="", ylab="",
#' xlim=xlim, ylim=ylim, asp=asp,
#' col=1:6, lwd=2, levels=-1000*1:6, drawlabels=FALSE)
#' # Show land
#' data(coastlineWorldFine, package="ocedata")
#' polygon(coastlineWorldFine[["longitude"]],
#' coastlineWorldFine[["latitude"]], col="lightgray")
#' # Indicate float positions.
#' points(index[["longitude"]], index[["latitude"]], pch=20)
#'
## # Example 7: TS plot for a particular Argo
## a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
## par(mar=c(3.3, 3.3, 1, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
## plot(a, which="TS")
##
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 7: Temperature QC plot for 25 cycles of a float in the Arabian Sea
## ais <- getIndex(filename="synthetic")
## sub <- subset(subset(ais, ID='2902123'), 50:75)
## profiles <- getProfiles(sub)
## argos <- readProfiles(profiles)
## par(mar=c(3.3, 3.3, 2, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
## plot(argos, which="QC") # defaults to temperature
## plot(argos, which="QC", QCControl=list(parameter="salinity"))
## plot(argos, which="QC", QCControl=list(parameter="salinity",dataStateIndicator=TRUE))
#'
#' # Example 7: Temperature profile of the 131st cycle of float with ID 2902204
#' a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
#' oldpar <- par(no.readonly=TRUE)
#' par(mfrow=c(1, 1))
#' par(mgp=c(2, 0.7, 0)) # mimic the oce::plotProfile() default
#' par(mar=c(1,3.5,3.5,2)) # mimic the oce::plotProfile() default
#' plot(a, which="profile")
#' par(oldpar)
#'
#' # Example 8: As Example 7, but showing temperature dependence on potential density anomaly.
#' a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
#' oldpar <- par(no.readonly=TRUE)
#' par(mgp=c(2, 0.7, 0)) # mimic the oce::plotProfile() default
#' par(mar=c(1,3.5,3.5,2)) # mimic the oce::plotProfile() default
#' plot(a, which="profile", profileControl=list(parameter="temperature", ytype="sigma0"))
#' par(oldpar)
#'
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 10: Summary plot
## par(mar=c(2, 3.3, 2, 1))
## a <- readProfiles(getProfiles(subset(getIndex(), ID=1901584)))
## plot(a, which="summary")
#'}
#'
#' @references
#' 1. Carval, Thierry, Bob Keeley, Yasushi Takatsuki, Takashi Yoshida, Stephen Loch,
#' Claudia Schmid, and Roger Goldsmith. Argo User's Manual V3.3. Ifremer, 2019.
#' `doi:10.13155/29825`
#'
#' 2. Kelley, D. E., Harbin, J., & Richards, C. (2021). argoFloats: An R package for analyzing
#' Argo data. Frontiers in Marine Science, (8), 636922.
#' \doi{10.3389/fmars.2021.635922}
#'
#' @importFrom graphics abline axis box contour par plot.window points polygon rect text
#'
#' @importFrom grDevices extendrange gray rgb
#' @importFrom utils data
#' @importFrom oce as.ctd colormap drawPalette imagep oceColorsGebco oce.plot.ts plotTS
#'
#' @export
#'
#' @aliases plot,argoFloats-method
#'
#' @return None (invisible NULL).
#'
#' @author Dan Kelley and Jaimie Harbin
setMethod(f="plot",
signature=signature("argoFloats"),
definition=function(x, # plot,argoFloats-method
which="map",
bathymetry=argoDefaultBathymetry(),
geographical=0,
xlim=NULL, ylim=NULL,
xlab=NULL, ylab=NULL,
type="p", cex=par("cex"), col=par("fg"), pch=par("pch"), bg=par("bg"),
eos="gsw",
mapControl=NULL,
profileControl=NULL,
QCControl=NULL,
summaryControl=NULL,
TSControl=NULL,
debug=0,
...)
{
if (!requireNamespace("oce", quietly=TRUE))
stop("must install.packages(\"oce\") for plot() to work")
debug <- if (debug > 3) 3 else max(0, floor(debug + 0.5))
argoFloatsDebug(debug, "plot(x, which=\"", which, "\") {\n", sep="", unindent=1, style="bold")
dots <- list(...)
if (!inherits(x, "argoFloats"))
stop("In plot() : method is only for objects of class \"argoFloats\"", call.=FALSE)
if (length(which) != 1)
stop("\"which\" must contain only one item")
if (!geographical %in% c(0, 1, 4))
stop("In plot() : \"geographical\" must be 0, 1, or 4", call.=FALSE)
istraj <- identical(x@metadata$subtype, "trajectories")
if (which == "map") {
data("coastlineWorld", package="oce", envir=environment())
argoFloatsDebug(debug, "map plot\n", sep="")
if (!istraj) {
longitude <- x[["longitude", debug=debug]]
latitude <- x[["latitude", debug=debug]]
}
if (istraj) {
lat1 <- as.numeric(x[["latitude_min"]])
lat2 <- as.numeric(x[["latitude_max"]])
lon1 <- as.numeric(x[["longitude_min"]])
lon2 <- as.numeric(x[["longitude_max"]])
latitude <- c(lat1, lat2)
longitude <- c(lon1, lon2)
}
n <- x[["length"]]
## Find type of each cycle, for colour-coding
cycleType <- rep("core", n)
cycleType[("849" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[("862" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[("864" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[!is.na(x@data$index$parameters)] <- "Bgc"
colLand <- mapControl$colLand
if (is.null(mapControl))
mapControl <- list(bathymetry=bathymetry, projection=FALSE, colLand="lightgray")
if (!"bathymetry" %in% names(mapControl))
mapControl$bathymetry <- TRUE
if (!"projection" %in% names(mapControl))
mapControl$projection <- FALSE
if (is.logical(mapControl$projection)) {
mapControl$projection <- if (mapControl$projection) "+proj=moll" else "none"
}
if (substr(mapControl$projection, 1, 6) != "+proj=" && mapControl$projection != "none")
stop("In plot,argoFloats-method(): mapControl$projection must start with \"+proj=\"", call.=FALSE)
if (mapControl$projection != "none") {
data("coastlineWorld", package="oce", envir=environment())
coastlineWorld <- get("coastlineWorld")
oce::mapPlot(coastlineWorld, col= "lightgray", projection=mapControl$projection, drawBox=FALSE)
if (!istraj) {
oce::mapPoints(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
col=if (is.null(col)) "white" else col,
pch=if (is.null(pch)) 21 else pch,
bg=if (is.null(bg)) "red" else bg,
type=if (is.null(type)) "p" else type,
...)
} else if (istraj) {
rect(lon1,lat1, lon2,lat2)
}
## warning("In plot,argoFloats-method(): projected maps do not (yet) show bathymetry", call.=FALSE)
return(invisible(NULL))
}
if (is.null(xlim))
if (!istraj) {
xlim <- extendrange(longitude)
} else if (istraj) {
xlim <- extendrange(range(lon1,lon2))
}
if (is.null(ylim))
if (!istraj) {
ylim <- extendrange(latitude)
} else if (istraj) {
ylim <- extendrange(range(lat1, lat2))
}
## Draw empty plot box, with axes, to set par("usr") for later use with bathymetry.
xlab <- if (is.null(xlab)) "" else xlab
ylab <- if (is.null(ylab)) "" else ylab
## Decode bathymetry
if (is.logical(mapControl$bathymetry)) {
drawBathymetry <- mapControl$bathymetry
bathymetry <- list(source="auto", contour=FALSE, colormap="auto", palette=TRUE)
} else if (is.list(mapControl$bathymetry)) {
drawBathymetry <- TRUE
if (!("source" %in% names(mapControl$bathymetry)))
stop("In plot() : \"bathymetry\" is a list, it must contain \"source\", at least", call.=FALSE)
if ("keep" %in% names(bathymetry))
warning("bathymetry$keep (an old argument) is no longer used\n")
if (is.null(bathymetry$contour))
bathymetry$contour <- FALSE
if (is.null(bathymetry$colormap))
bathymetry$colormap <- "auto"
if (is.null(bathymetry$palette))
bathymetry$palette <- TRUE
} else {
stop("In plot() : \"bathymetry\" must be either a logical or a list value", call.=FALSE)
}
if (!is.logical(bathymetry$palette))
stop("In plot() : \"bathymetry$palette\" must be a logical value", call.=FALSE)
argoFloatsDebug(debug, "drawBathymetry calculated to be ", drawBathymetry, "\n", sep="")
asp <- 1 / cos(pi/180*mean(range(latitude, na.rm=TRUE)))
argoFloatsDebug(debug, "asp=", asp, "\n", sep="")
if (drawBathymetry) {
argoFloatsDebug(debug, "handling bathymetry\n", sep="")
if (!requireNamespace("oce", quietly=TRUE))
stop("must install.packages(\"oce\") to plot with bathymetry")
## Handle bathymetry file downloading (or the use of a supplied value)
bathy <- NULL
if (is.character(bathymetry$source) && bathymetry$source == "auto") {
argoFloatsDebug(debug, "must either download bathymetry or use existing file\n", sep="")
if (!bathymetry$contour && bathymetry$palette) {
omar <- par("mar")
par(mar=omar + c(0, 0, 0, 2.75))
on.exit(par(mar=omar))
argoFloatsDebug(debug, "temporarily set par(mar=c(", paste(par("mar"), collapse=", "), ")) to allow for the palette\n", sep="")
}
argoFloatsDebug(debug, "using plot.window() to determine area for bathymetry download, with\n",
" extendrange(longitude)=c(", paste(extendrange(longitude), collapse=","), ")\n",
" extendrange(latitude)=c(", paste(extendrange(latitude), collapse=","), ")\n",
" xlim=c(", paste(xlim, collapse=","), ")\n",
" ylim=c(", paste(ylim, collapse=","), ")\n",
" asp=", asp, "\n", sep="")
plot.window(xlim=xlim, ylim=ylim,
xaxs="i", yaxs="i",
asp=asp,
xlab=xlab, ylab=ylab)
usr <- par("usr")
argoFloatsDebug(debug, "after using plot.window(), usr=c(", paste(round(usr, 4), collapse=", "), ")\n", sep="")
latitudeSpan <- usr[4] - usr[3]
longitudeSpan <- usr[2] - usr[1]
Dlon <- 0.05 * longitudeSpan
Dlat <- 0.05 * latitudeSpan
argoFloatsDebug(debug, "Dlon=", Dlon, " (5% of longitude span, which is ", longitudeSpan, ")\n", sep="")
argoFloatsDebug(debug, "Dlat=", Dlat, " (5% of latitude span, which is ", latitudeSpan, ")\n", sep="")
resolution <- as.integer(round(1 + 60 * latitudeSpan / 400))
if (resolution < 1)
stop("calculated resolution (=", resolution, ") cannot be under 1 minute")
argoFloatsDebug(debug, "resolution=", resolution, "\n", sep="")
## Notice the use of 0.1 and 1 digit in lon and lat, which means
## we are enlarging the field by 10km in both lon and lat.
west <- round(usr[1]-Dlon, 1)
east <- round(usr[2]+Dlon, 1)
south <- round(usr[3]-Dlat, 1)
north <- round(usr[4]+Dlat, 1)
argoFloatsDebug(debug, "will get topo in domain west=", west, ", east=", east,
", south=", south, ", north=", north, "\n")
topo <- try(oce::read.topo(oce::download.topo(west=west,
east=east,
south=south,
north=north,
resolution=resolution,
debug=debug)),
silent=FALSE)
if (inherits(topo, "try-error")) {
warning("could not download bathymetry from NOAA server\n")
drawBathymetry <- FALSE
} else {
bathy <- -topo[["z"]]
dimBathy <- dim(bathy)
argoFloatsDebug(debug, "grid size", paste(dimBathy, collapse="x"), "\n")
bathy <- as.integer(bathy)
dim(bathy) <- dimBathy
rownames(bathy) <- topo[["longitude"]]
colnames(bathy) <- topo[["latitude"]]
class(bathy) <- "bathy"
}
} else if (inherits(bathymetry$source, "bathy")) {
argoFloatsDebug(debug, "using supplied bathymetry$source\n", sep="")
bathy <- bathymetry$source
} else if (inherits(bathymetry$source, "oce") && inherits(bathymetry$source, "topo")) {
argoFloatsDebug(debug, "using oce-style topo object (converted to NOAA bathy)\n", sep="")
dim <- dim(bathymetry$source[["z"]])
## Note the negative sign, to get depth.
bathy <- matrix(-as.integer(bathymetry$source[["z"]]), nrow=dim[1], ncol=dim[2])
rownames(bathy) <- bathymetry$source[["longitude"]]
colnames(bathy) <- bathymetry$source[["latitude"]]
class(bathy) <- "bathy"
} else {
stop("cannot determine bathymetry data source")
}
argoFloatsDebug(debug, "drawBathymetry = ", drawBathymetry, " (after trying to get bathymetry)\n", sep="")
## Handle colormap
if (drawBathymetry && !inherits(bathy, "try-error") && !bathymetry$contour && is.character(bathymetry$colormap) && length(bathymetry$colormap) == 1 && bathymetry$colormap == "auto") {
argoFloatsDebug(debug, "bathy class: \"", class(bathy), "\"\n", sep="")
deepest <- max(bathy, na.rm=TRUE)
argoFloatsDebug(debug, "bathy range: ", paste(range(bathy, na.rm=TRUE), collapse=" to "), "\n")
argoFloatsDebug(debug, "setting a default colormap for ", deepest, "m to 0m depth\n", sep="")
bathymetry$colormap <- oce::colormap(zlim=c(0, deepest),
col=function(n) rev(oce::oceColorsGebco(n)))
}
## Handle optional drawing of the palette. Note space for one line of text is removed
## from the LHS of the palette and added to the RHS. This is because we know that
## there is no axis to the left of the palette, so we do not need space for one.
## We recover the stolen space by putting it back at the RHS, where it can be
## useful, especially if a map plot has another plot to its right.
if (drawBathymetry && !inherits(bathy, "try-error") && !bathymetry$contour && bathymetry$palette) {
argoFloatsDebug(debug, "drawing a bathymetry palette\n")
omai <- par("mai") # save this, since oce::drawPalette() may change it
oce::drawPalette(colormap=bathymetry$colormap, cex=if (is.null(cex)) 1 else par("cex.axis"))
on.exit(par(mai=omai))
} else {
argoFloatsDebug(debug, "not drawing a bathymetry palette, as instructed (or failed bathymetry download)\n")
}
}
argoFloatsDebug(debug, "initial xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
## Extend to world edge, if we are close and scale is large
if (diff(xlim) > 340) {
if (xlim[1] < -170) xlim[1] <- -180
if (xlim[2] > 170) xlim[2] <- 180
}
if (diff(ylim) > 160) {
if (ylim[1] < -70) ylim[1] <- -90
if (ylim[2] > 70) ylim[2] <- 90
}
argoFloatsDebug(debug, "after to-world-edge, xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
## Trim to world domain
xlim <- ifelse(xlim < -180, -180, ifelse(xlim > 180, 180, xlim))
ylim <- ifelse(ylim < -80, -90, ifelse(ylim > 80, 90, ylim))
argoFloatsDebug(debug, "after domain-trim, xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
argoFloatsDebug(debug, "about to plot, with\n",
" extendrange(longitude)=c(", paste(extendrange(longitude), collapse=","), ")\n",
" extendrange(latitude)=c(", paste(extendrange(latitude), collapse=","), ")\n",
" xlim=c(", paste(xlim, collapse=","), ")\n",
" ylim=c(", paste(ylim, collapse=","), ")\n",
" asp=", asp, "\n", sep="")
plot(xlim, ylim,
xaxs="i", yaxs="i",
asp=asp,
xlab=xlab, ylab=ylab, type="n", axes=FALSE)
## Draw box and axes. Note that we cannot just use axis() because we
## do not want to show e.g. lat beyond 90N.
usrTrimmed <- par("usr")
usrTrimmed[1] <- max(-180, usrTrimmed[1])
usrTrimmed[2] <- min(180, usrTrimmed[2])
usrTrimmed[3] <- max(-90, usrTrimmed[3])
usrTrimmed[4] <- min(90, usrTrimmed[4])
rect(usrTrimmed[1], usrTrimmed[3], usrTrimmed[2], usrTrimmed[4], lwd=1)
## Construct axes
xat <- pretty(usrTrimmed[1:2])
if (diff(range(xat)) > 300)
xat <- seq(-180, 180, 45)
yat <- pretty(usrTrimmed[3:4])
if (diff(range(yat)) > 150)
yat <- seq(-90, 90, 45)
if (geographical == 0) {
xlabels <- xat
ylabels <- yat
} else if (geographical == 1) {
xlabels <- abs(xat)
ylabels <- abs(yat)
} else if (geographical == 4) {
xlabels <- paste(abs(xat), ifelse(xat < 0, "W", ifelse(xat > 0, "E", "")), sep="")
ylabels <- paste(abs(yat), ifelse(yat < 0, "S", ifelse(yat > 0, "N", "")), sep="")
} else {
stop("In plot,argoFloats-method() : programming error: \"geographical\" must be 0, 1, or 4", call.=FALSE)
}
axis(1, at=xat, labels=xlabels, pos=usrTrimmed[3])
axis(2, at=yat, labels=ylabels, pos=usrTrimmed[1])
axis(3, at=xat, labels=FALSE, pos=usrTrimmed[4])
axis(4, at=yat, labels=FALSE, pos=usrTrimmed[2])
## argoFloatsMapAxes(geographical=geographical)
if (drawBathymetry) {
if (bathymetry$contour) {
argoFloatsDebug(debug, "indicating bathymetry with contours\n")
contour(as.numeric(rownames(bathy)),
as.numeric(colnames(bathy)),
bathy,
levels=c(100, 200, 500, seq(1e3, 10e3, 1e3)),
labels=c("100m", "200m", "500m", paste(1:10, "km", sep="")),
labcex=0.9, add=TRUE)
} else {
argoFloatsDebug(debug, "indicating bathymetry with an image\n")
oce::imagep(as.numeric(rownames(bathy)),
as.numeric(colnames(bathy)),
bathy, colormap=bathymetry$colormap, add=TRUE)
}
rect(usrTrimmed[1], usrTrimmed[3], usrTrimmed[2], usrTrimmed[4], lwd=1)
}
if (is.null(col)) {
argoFloatsDebug(debug, "defaulting col to indicate cycle type\n")
col <- rep("", n)
col[cycleType == "core"] <- colDefaults$core
col[cycleType == "Bgc"] <- colDefaults$bgc
col[cycleType == "deep"] <- colDefaults$deep
}
## Select coastline. Unlike in oce::plot,coastline-method, we base our choice
## on just the distance spanned in the north-south direction.
ocedataIsInstalled <- requireNamespace("ocedata", quietly=TRUE)
if (ocedataIsInstalled) {
usr <- par("usr")
mapSpan <- diff(range(latitude, na.rm=TRUE)) * 111 # km
##? l <- oce::geodDist(pinlon(usr[1]), pinlat(usr[3]), pinlon(usr[1]), pinlat(usr[4]))
##? r <- oce::geodDist(pinlon(usr[2]), pinlat(usr[3]), pinlon(usr[2]), pinlat(usr[4]))
##? b <- oce::geodDist(pinlon(usr[1]), pinlat(usr[1]), pinlon(usr[2]), pinlat(usr[1]))
##? t <- oce::geodDist(pinlon(usr[1]), pinlat(usr[4]), pinlon(usr[2]), pinlat(usr[4]))
##? mapSpan <- max(l, r, b, t) # largest length [km]
C <- 2 * 3.14 * 6.4e3 # circumferance of earth [km]
argoFloatsDebug(debug, "mapSpan=", mapSpan, ", C=", C, "\n")
if (mapSpan < 500) {
argoFloatsDebug(debug, "using coastlineWorldFine from ocedata package\n")
data("coastlineWorldFine", package="ocedata", envir=environment())
coastlineWorldFine <- get("coastlineWorldFine")
polygon(coastlineWorldFine[["longitude"]], coastlineWorldFine[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
} else if (mapSpan < C / 4) {
argoFloatsDebug(debug, "using coastlineWorldMedium from ocedata package\n")
data("coastlineWorldMedium", package="ocedata", envir=environment())
coastlineWorldMedium <- get("coastlineWorldMedium")
polygon(coastlineWorldMedium[["longitude"]], coastlineWorldMedium[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
} else {
argoFloatsDebug(debug, "using coastlineWorld from oce package, since the span is large\n")
data("coastlineWorld", package="oce", envir=environment())
coastlineWorld <- get("coastlineWorld")
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
}
if (is.null(pch))
pch <- 21
type <- if (is.null(type)) "p" else type
if (length(pch) == 1 && pch == 21 && !istraj) {
if (type != "l") {
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= type,
...)
} else if (type == "l") {
index1 <- subset(x, ID=unique(x[["ID"]][1]), silent=TRUE)
longitude <- index1[["longitude"]]
latitude <- index1[["latitude"]]
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= "l",
...)
for (i in unique(x[["ID"]])[2:length(unique(x[["ID"]]))]) {
index2 <- subset(x, ID=i, silent=TRUE)
points(unlist(index2[["longitude"]]), unlist(index2[["latitude"]]),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= "l",
...)
}
}
} else if (!istraj) {
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
col=col,
type=if (is.null(type)) "p" else type,
...)
} else if (istraj) {
rect(lon1,lat1, lon2,lat2)
}
} else {
argoFloatsDebug(debug, "using coastlineWorld from oce package, since the ocedata package is not installed\n")
data("coastlineWorld", package="ocedata", envir=environment())
coastlineWorld <- get("coastlineWorld")
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col=col)
}
} else if (which == "summary" && !istraj) {
argoFloatsDebug(debug, "summary plot\n", sep="")
if (x[["length"]] < 2)
stop("In plot,argoFloats-method() : cannot draw a summary plot with only one float cycle", call.=FALSE)
if (!all(names(dots) %in% c("summaryControl")))
stop(names(dots)," is not a valid argument for this plot type. Try summaryControl instead.")
if (is.null(summaryControl))
summaryControl <- list(items=c("dataStateIndicator", "longitude", "latitude", "length", "deepest"))
if (!"items" %in% names(summaryControl))
stop("summaryControl must be a list containing a character vector named \"items\"")
if (!(summaryControl$items) %in% c("dataStateIndicator", "longitude", "latitude", "length", "deepest"))
stop("The names of the items must be dataStateIndicator, longitude, latitude, length, or deepest, not ", paste(summaryControl$items, collapse=" "))
items <- summaryControl$items
nitems <- length(items)
if (nitems) {
omfrow <- par("mfrow")
omar <- par("mar")
omgp <- par("mgp")
par(mar=c(3,3,2,1.5), mgp=c(1.75, 0.5, 0))
if (nitems < 3)
par(mfrow=c(nitems, 1))
else if (nitems == 4)
par(mfrow=c(2, 2))
else
par(mfrow=c(3, 2))
on.exit(par(mfrow=omfrow, mar=omar, mgp=omgp))
time <- as.POSIXct(unlist(x[["time"]]), origin="1970-01-01", tz="UTC")
o <- order(time)
for (iitem in seq_len(nitems)) {
argoFloatsDebug(debug, "items[", iitem, "]=\"", items[iitem], "\"\n", sep="")
if (items[iitem] == "dataStateIndicator") {
y <- unlist(x[["dataStateIndicator"]])
argoFloatsDebug(debug, oce::vectorShow(y))
if (length(y)) {
u <- sort(unique(y))
yy <- seq_along(u)
oce::oce.plot.ts(range(time), range(yy), ylab="Data State Ind.",
drawTimeRange=FALSE, type="n", axes=FALSE,
xaxs="i", ylim=range(0.8, (max(yy)+0.2)),
mar=par("mar"), mgp=par("mgp"))
abline(h=seq_along(u), col="gray")
points(time[o], factor(y)[o])
oce::oce.axis.POSIXct(side=1, drawTimeRange=FALSE)
box()
axis(side=2, at=yy, labels=u)
}
} else if (items[iitem] == "length") {
y <- sapply(x[["argos"]], function(a) length(a[["pressure"]]))
oce::oce.plot.ts(time[o], y[o], ylab="Length", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "longitude") {
y <- sapply(x[["longitude"]], function(a) a[1])
oce::oce.plot.ts(time[o], y[o], ylab="Longitude", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "latitude") {
y <- sapply(x[["latitude"]], function(a) a[1])
oce::oce.plot.ts(time[o], y[o], ylab="Latitude", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "deepest") {
y <- sapply(x[["argos"]], function(a) max(a[["pressure"]], na.rm=TRUE))
oce::oce.plot.ts(time[o], y[o], ylab="Max Pres.", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
}
if (iitem == 1) {
if (1 == length(unique(x[["ID"]])))
axis(side=3, at=time[o], labels=x[["cycle"]][o],
cex.axis=max(0.75, 0.75*par("cex")))
}
}
}
} else if (which == "summary" && istraj) {
stop("In plot,argoFloats-method(): cannot draw summary plot for subtype = trajectories", call.=FALSE)
} else if (which == "TS") {
#cat("next is cex before call:\n");cat(vectorShow(cex))
#message("main debug=", debug)
plotArgoTS(x=x, xlim=xlim, ylim=ylim,
type=type, cex=cex, col=col, pch=pch, bg=bg,
eos=eos, TSControl=TSControl, debug=debug-1L)
argoFloatsDebug(debug, "} # plot()\n", sep="", unindent=1, style="bold")
return()
#<20230805> #return(do.call("plotArgoTS", args))
#<20230805> stop("FIXME: cannot get here")
#<20230805> argoFloatsDebug(debug, "TS plot\n", sep="")
#<20230805> if ((x[["type"]] != "argos"))
#<20230805> stop("In plot,argoFloats-method(): x must have been created by readProfiles()", call.=FALSE)
#<20230805> if (!(eos %in% c("gsw", "unesco")))
#<20230805> stop("In plot,argoFloats-method(): eos must be \"gsw\" or \"unesco\", not \"", eos, "\"", call.=FALSE)
#<20230805> # compute cycle, in case we need that
#<20230805> cycleIndex <- unlist(lapply(seq_len(x[["length"]]),
#<20230805> function(icycle)
#<20230805> rep(icycle, length(x@data$argos[[icycle]]@data$pressure))))
#<20230805> salinity <- unlist(x[["salinity", debug=debug-1L]])
#<20230805> temperature <- unlist(x[["temperature", debug=debug-1L]])
#<20230805> pressure <- unlist(x[["pressure", debug=debug-1L]])
#<20230805> ## Use byLevel to repeat the latitude and longitude values across
#<20230805> ## the depths in each profile, so that the resultant vector
#<20230805> ## will match the salinity, temperature and pressure vectors.
#<20230805> latitude <- unlist(x[["latitude", "byLevel", debug=debug-1L]])
#<20230805> longitude <- unlist(x[["longitude", "byLevel", debug=debug-1L]])
#<20230805> # Interpolate across NA longitudes (required for traj data, to get TS plot)
#<20230805> n <- length(longitude)
#<20230805> if (any(is.na(longitude)))
#<20230805> longitude <- approx(1:n, longitude, 1:n)$y
#<20230805> if (any(is.na(latitude)))
#<20230805> latitude <- approx(1:n, latitude, 1:n)$y
#<20230805> ctd <- oce::as.ctd(salinity=salinity,
#<20230805> temperature=temperature,
#<20230805> pressure=pressure,
#<20230805> latitude=latitude,
#<20230805> longitude=longitude)
#<20230805> ctd@data$cycleIndex <- cycleIndex
#<20230805> if (is.null(TSControl)) {
#<20230805> argoFloatsDebug(debug, "defaulting TSControl\n")
#<20230805> TSControl <- list(colByCycle=NULL)
#<20230805> }
#<20230805> if (is.null(TSControl$colByCycle)) {
#<20230805> argoFloatsDebug(debug, "TSControl does not contain colByCycle\n")
#<20230805> if (is.null(col)) {
#<20230805> if (which == "TS") {
#<20230805> col <- "flags"
#<20230805> } else {
#<20230805> col <- rgb(0, 0, 1, 0.5)
#<20230805> }
#<20230805> }
#<20230805> } else {
#<20230805> argoFloatsDebug(debug, "TSControl contains colByCycle\n")
#<20230805> # Ignore "col" if TSControl contains "colByCycle"
#<20230805> colByCycle <- TSControl$colByCycle
#<20230805> cycle <- unlist(x[["cycle", debug=debug-1L]])
#<20230805> lengths <- sapply(x[["argos"]], function(cc) length(cc[["pressure"]]))
#<20230805> # Increase the col length, so e.g. TSControl=list(colByCycle=1:2) will alternate colours
#<20230805> colByCycle <- rep(colByCycle, length.out=length(cycle))
#<20230805> col <- unlist(lapply(seq_along(cycle), function(i) rep(colByCycle[i], lengths[i])))
#<20230805> }
#<20230805> if (is.null(cex)) {
#<20230805> cex <- par("cex")
#<20230805> argoFloatsDebug(debug, "defaulting to cex=", cex, "\n", sep="")
#<20230805> }
#<20230805> if (is.null(pch)) {
#<20230805> pch <- 20
#<20230805> argoFloatsDebug(debug, "defaulting to pch=", pch, "\n", sep="")
#<20230805> }
#<20230805> if (col[1] == "flags") {
#<20230805> argoFloatsDebug(debug, "col is \"flags\"\n")
#<20230805> salinityFlag <- unlist(x[["salinityFlag"]])
#<20230805> temperatureFlag <- unlist(x[["temperatureFlag"]])
#<20230805> goodT <- temperatureFlag %in% c(1, 2, 5, 8)
#<20230805> goodS <- salinityFlag %in% c(1, 2, 5, 8)
#<20230805> good <- goodS & goodT
#<20230805> okT <- temperatureFlag %in% c(0)
#<20230805> okS <- salinityFlag %in% c(0)
#<20230805> ok <- okS & okT
#<20230805> col <- ifelse(good, "black", ifelse(ok, "gray", "red"))
#<20230805> if (pch == 21)
#<20230805> bg <- ifelse(good, "black", ifelse(ok, "gray", "red"))
#<20230805> }
#<20230805> argoFloatsDebug(debug, "about to plotTS()\n")
#<20230805> if (!is.null(TSControl$colByCycle)) {
#<20230805> argoFloatsDebug(debug, "colorizing by index\n")
#<20230805> cycles <- unique(cycleIndex)
#<20230805> ncycles <- length(cycles)
#<20230805> if (eos == "gsw") {
#<20230805> Slim <- range(ctd[["SA"]], na.rm=TRUE)
#<20230805> Tlim <- range(ctd[["CT"]], na.rm=TRUE)
#<20230805> } else {
#<20230805> Slim <- range(ctd[["salinity"]], na.rm=TRUE)
#<20230805> Tlim <- range(ctd[["temperature"]], na.rm=TRUE)
#<20230805> }
#<20230805> argoFloatsDebug(debug, "set Slim=", Slim, "\n")
#<20230805> argoFloatsDebug(debug, "set Tlim=", Tlim, "\n")
#<20230805> if (length(longitude) != length(salinity))
#<20230805> longitude <- rep(longitude[1], length(salinity))
#<20230805> if (length(latitude) != length(salinity))
#<20230805> latitude <- rep(latitude[1], length(salinity))
#<20230805> col <- TSControl$colByCycle
#<20230805> if (length(col) != ncycles)
#<20230805> col <- rep(col, length.out=ncycles)
#<20230805> if (length(type) != ncycles)
#<20230805> type <- rep(type, length.out=ncycles)
#<20230805> if (length(pch) != ncycles)
#<20230805> pch <- rep(pch, length.out=ncycles)
#<20230805> for (i in seq_len(ncycles)) {
#<20230805> argoFloatsDebug(debug, " handling cycle ", i, " of ", ncycles, ", which has ",
#<20230805> length(salinity[i==cycleIndex]), " points\n", sep="")
#<20230805> ctd <- oce::as.ctd(
#<20230805> salinity=salinity[i == cycleIndex],
#<20230805> temperature=temperature[i == cycleIndex],
#<20230805> pressure=pressure[i == cycleIndex],
#<20230805> latitude=latitude[i == cycleIndex],
#<20230805> longitude=longitude[i == cycleIndex])
#<20230805> if (i == 1L) {
#<20230805> plotTS(ctd, Slim=Slim, Tlim=Tlim,
#<20230805> type=type[i], col=col[i], pch=pch[i])
#<20230805> } else {
#<20230805> plotTS(ctd, add=TRUE,
#<20230805> type=type[i], col=col[i], pch=pch[i])
#<20230805> }
#<20230805> }
#<20230805> } else {
#<20230805> argoFloatsDebug(debug, "not colorizing by index\n")
#<20230805> oce::plotTS(ctd,
#<20230805> cex=cex,
#<20230805> bg=bg,
#<20230805> col=col,
#<20230805> pch=pch,
#<20230805> mar=par("mar"), mgp=par("mgp"), eos=eos,
#<20230805> type=if (is.null(type)) "p" else type, ...)
#<20230805> }
} else if (which == "QC" && !istraj) {
if (x[["type"]] != "argos")
stop("In plot,argoFloats-method(): The type of x must be \"argos\"", call.=FALSE)
IDs <- x[["ID"]]
nID <- length(unique(IDs))
if (nID != 1)
stop("In plot,argoFloats-method(): It is only possible to plot a QC of a single ID", call.=FALSE)
knownParameters <- names(x[[1]]@metadata$flags) # FIXME: is it possible that later cycles have different flags?
if (!all(names(dots) %in% c("QCControl")))
stop(names(dots)," is not a valid argument for this plot type. Try QCControl instead.")
if (is.null(QCControl)) {
QCControl <- list(parameter="temperature")
}
if (!is.list(QCControl))
stop("In plot,argoFloats-method(): QCControl must be a list")
if (!all(names(QCControl) %in% c("parameter", "dataStateIndicator")))
stop("QCControl must contain only elements named \"parameter\" and \"dataStateIndicator\", but its elements are named \"", paste(names(QCControl), collapse="\", \""), "\"")
if (!"parameter" %in% names(QCControl))
QCControl$parameter <- "temperature"
if (!"dataStateIndicator" %in% names(QCControl))
QCControl$dataStateIndicator <- FALSE
if (!is.logical(QCControl$dataStateIndicator))
stop("QCControl element 'dataStateIndicator' must be logical")
if (length(QCControl) != 2)
stop("In plot,argoFloats-method(): QCControl must contain only elements named \"parameter\" and \"dataStateIndicator\"")
if (!(QCControl$parameter %in% knownParameters))
stop("In plot,argoFloats-method(): QCControl$parameter '", QCControl$parameter, "' not found. Try one of: \"", paste(sort(knownParameters), collapse="\", \""), "\"", call.=FALSE)
qf <- function(x) {
# qf returns 100 if data are all "good" = 1 or "probably good" = 2 or "changed" = 5 or "estimated" = 8
flag <- x[[paste0(QCControl$parameter, "Flag")]]
100 * sum(1 == flag | 2 == flag | 5 == flag | 8 == flag, na.rm=TRUE) / length(flag)
}
meanf <- function(x)
mean(x[[QCControl$parameter]], na.rm=TRUE)
time <- oce::numberAsPOSIXct(unlist(lapply(x[["argos"]], function(x) x[["time"]])))
q <- unlist(lapply(x[["argos"]], qf))
m <- unlist(lapply(x[["argos"]], meanf))
omfrow <- par("mfrow")
omar <- par("mar")
par(mar=c(3, 3, 2, 1))
par(mfrow=c(if (QCControl$dataStateIndicator) 3 else 2, 1))
on.exit(par(mfrow=omfrow, mar=omar))
if (any(is.finite(q))) {
o <- order(time) # cycles are not time-ordered in index files
oce::oce.plot.ts(time[o], q[o], ylab=paste(QCControl$parameter, "% Good"), drawTimeRange=FALSE, type="l", mar=par("mar"), xaxs="i", cex.lab=0.7, cex.axis=0.7)
points(time[o], q[o], col=ifelse(q[o] < 50, "red", "black"), pch=20, cex=1)
abline(h=50, col="red", lty="dashed")
if (1 == length(unique(x[["ID"]])))
axis(side=3, at=time[o], labels=x[["cycle"]][o], cex.axis=0.7, cex.lab=0.7)
if (QCControl$dataStateIndicator) {
y <- unlist(x[["dataStateIndicator"]])
if (length(y)) {
u <- sort(unique(y))
yy <- seq_along(u)
oce::oce.plot.ts(range(time), range(yy), ylab="Data State Ind.",
drawTimeRange=FALSE, type="n", axes=FALSE,
mar=par("mar"), mgp=par("mgp"),
xaxs="i", ylim=range(0.8, (max(yy)+0.2)))
abline(h=seq_along(u), col="gray")
points(time[o], factor(y)[o])
oce::oce.axis.POSIXct(side=1, drawTimeRange=FALSE)
box()
axis(side=2, at=yy, labels=u)
} else {
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0.5, 0.5, "x lacks dataStateIndicator")
}
}
oce::oce.plot.ts(time[o], m[o], ylab=paste(QCControl$parameter, "Mean"), drawTimeRange=FALSE, type="l",
mar=par("mar"), cex.lab=0.7, cex.axis=0.7)
points(time[o], m[o], col=ifelse(q[o] < 50, "red", "black"), pch=20, cex=1)
} else {
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0, 0.5, paste(' No', QCControl$parameter, "flags available"), pos=4)
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0, 0.5, paste(' No', QCControl$parameter, "flags available"), pos=4)
}
} else if (which == "QC" && istraj) {
stop("In plot,argoFloats-method(): cannot draw QC plot for subtype = trajectories", call.=FALSE)
} else if (which == "profile") {
if (x[["type"]] != "argos")
stop("In plot,argoFloats-method(): The type of x must be \"argos\"", call.=FALSE)
if (is.null(type))
type <- "p"
if (is.null(profileControl)) {
if ("parameter" %in% names(dots)) {
warning("accepting \"parameter\" as a separate argument, until 2020-dec-01. After that, you must use e.g. profileControl=list(parameter=", dots$parameter, ")")
profileControl <- list(parameter=dots$parameter)
} else {
profileControl <- list(parameter="temperature")
}
profileControl$ytype <- "pressure"
profileControl$connect <- FALSE
}
if (!is.list(profileControl))
stop("In plot,argoFloats-method(): profileControl must be a list")
if (!"ytype" %in% names(profileControl))
profileControl$ytype <- "pressure"
if (!"parameter" %in% names(profileControl))
profileControl$parameter <- "temperature"
if (!"connect" %in% names(profileControl))
profileControl$connect <- TRUE
if (length(profileControl) != 3)
stop("In plot,argoFloats-method(): profileControl must contain only three elements, \"parameter\", \"ytype\" and \"connect\".")
if (!profileControl$ytype %in% c("pressure", "sigma0"))
stop("In plot,argoFloats-method(): profileControl$ytype must be \"pressure\" or \"sigma0\", not \"", profileControl$ytype, "\"")
N <- length(x[["argos"]])
## The known parameter names include not just the things stored in the
## data slot (of *any* of the profiles), but also some computable things. We
## will always have longitude and latitude, so there's no need to check for them.
## A core profile will always get the computable things, but bgc profiles
## may lack e.g. salinity and temperature, so all we need to check for is
## salinity, temperature, and pressure.
knownParameters <- unique(unlist(lapply(seq_len(N), function(i) names(x[[i]][["data"]]))))
if (all(c("salinity", "temperature", "pressure") %in% knownParameters)) {
knownParameters <- c(knownParameters,
"SA", "CT", "sigma0", "sigma1", "sigma2", "sigma3", "sigma4", "sigmaTheta",
"theta",
"spice",
"N2")
}
argoFloatsDebug(debug, "knownParameters: \"", paste(sort(knownParameters), collapse="\", \""), "\".\n", sep="")
if (!(profileControl$parameter %in% knownParameters))
stop("In plot,argoFloats-method(): profileControl$parameter=\"", profileControl$parameter, "\" is not in the dataset, or calculable from that dataset; try one of the following: \"", paste(sort(knownParameters), collapse='", "'), "\".", call.=FALSE)
y <- lapply(1:N, function(i) x[[i]][[profileControl$ytype]])
variable <- lapply(1:N, function(i) x[[i]][[profileControl$parameter]])
nn <- unlist(lapply(1:N, function(i) prod(dim(x[[i]][[profileControl$parameter]]))))
Y <- NULL
VARIABLE <- NULL
argoFloatsDebug(debug, "number of profiles: ", N, "\n")
for (i in seq_len(N)) {
if (nn[i] > 0) {
Y <- c(Y, NA, y[[i]])
VARIABLE <- c(VARIABLE, NA, variable[[i]])
}
}
if (type %in% c("l", "o", "b") && profileControl$connect) {
ok <- is.finite(VARIABLE) & is.finite(Y)
VARIABLE <- VARIABLE[ok]
Y <- Y[ok]
}
if (is.null(cex)) {
cex <- par("cex")
argoFloatsDebug(debug, "Profile plot defaulting to cex=", cex, "\n")
}
if (is.null(pch)) {
pch <- 20
argoFloatsDebug(debug, "Profile plot defaulting to pch=", pch, "\n")
}
plot(VARIABLE, Y,
xlim=xlim,
ylim=if (is.null(ylim)) rev(range(Y, na.rm=TRUE)) else ylim,
axes=FALSE,
ylab="", xlab="", # draw axes later, in oceanographic 'profile' locations
cex=cex,
type=if(is.null(type)) "l" else type,
col=if (is.null(col)) par("col") else col,
pch=pch,
...)
box()
axis(2)
axis(3)
ylab <- oce::resizableLabel(if (profileControl$ytype == "pressure") "p" else "sigma0", axis="y")
mtext(ylab, side=2, line=par("mgp")[1], cex=par("cex"))
# Rename for oce::resizableLabel()
if (profileControl$parameter == "salinity")
profileControl$parameter <- "S"
xlab <- oce::resizableLabel(if (profileControl$parameter == "oxygen") "oxygen umol/kg" else profileControl$parameter, axis="x")
mtext(xlab, side=3, line=par("mgp")[1], cex=par("cex"))
} else {
stop("In plot,argoFloats-method():cannot handle which=\"", which, "\"; see ?\"plot,argoFloats-method\"", call.=FALSE)
}
argoFloatsDebug(debug, "} # plot()\n", sep="", unindent=1, style="bold")
invisible(NULL)
}
)
Try the argoFloats package in your browser
Any scripts or data that you put into this service are public.
argoFloats documentation built on Oct. 18, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# vim:textwidth=128:expandtab:shiftwidth=4:softtabstop=4
colDefaults <- list(core="7", bgc="#05f076", deep="6")
#' Plot an argoFloats Object
#'
#' The action depends on the `type` of the object, and
#' this is set up by the function that created the object;
#' see \dQuote{Details}. These are basic plot styles, with
#' somewhat limited scope for customization. Since the data with
#' [argoFloats-class] objects are easy to extract, users should
#' not find it difficult to create their own plots to meet a
#' particular aesthetic. See \dQuote{Examples} and Kelley et al. (2021)
#' for more plotting examples.
#'
#' The various plot types are as follows.
#'
#' * For `which="map"`, a map of profile locations is created if subtype
#' is equal to cycles, or a rectangle highlighting the trajectory of a
#' float ID is created when subtype is equal to trajectories. This
#' only works if the `type` is `"index"` (meaning that `x` was created
#' by [getIndex()] or a subset of such an object, created with
#' [subset,argoFloats-method()]), or `argos` (meaning that
#' `x` was created with [readProfiles()]. The plot range is
#' auto-selected. If the `ocedata` package
#' is available, then its `coastlineWorldFine` dataset is used to draw
#' a coastline (which will be visible only if the plot region
#' is large enough); otherwise, if the `oce` package is available, then its
#' `coastlineWorld` dataset is used.
#' The `bathymetry` argument controls whether (and how) to draw a map underlay
#' that shows water depth. There are three possible values for `bathymetry`:
#' 1. `FALSE` (the default, via [argoDefaultBathymetry()]), meaning not to draw bathymetry;
#' 2. `TRUE`, meaning to draw bathymetry as an
#' image, using data downloaded with [oce::download.topo()].
#' Example 4 illustrates this, and also shows how to adjust
#' the margins after the plot, in case there is a need to add
#' extra graphical elements using [points()], [lines()], etc.
#' 3. A list with items controlling both the bathymetry data and its
#' representation in the plot (see Example 4). Those items are:
#'
#' 1. `source`, a mandatory value that is either
#' (a) the string `"auto"` (the default) to use
#' [oce::download.topo()] to download the data
#' or (b) a value returned by [oce::read.topo()].
#' 2. `contour`, an optional logical value (with `FALSE` as the default)
#' indicating whether to represent
#' bathymetry with contours (with depths of 100m, 200m, 500m shown,
#' along with 1km, 2km up to 10km), as opposed to an image;
#' 3. `colormap`, ignored if `contour` is `TRUE`,
#' an optional value that is either the string `"auto"` (the default)
#' for a form of GEBCO colors computed with
#' [oce::oceColorsGebco()], or a value computed with [oce::colormap()]
#' applied to the bathymetry data; and
#' 4. `palette`, ignored if `contour` is `TRUE`,
#' an optional logical value (with `TRUE` as the default)
#' indicating whether to draw a depth-color
#' palette to the right of the plot.
#' Note that the default value for `bathymetry` is set by a call
#' to [argoDefaultBathymetry()], and that this second function can only handle possibilities
#' 1 and 2 above.
#'
#' * For `which="profile"`, a profile plot is created, showing the variation of
#' some quantity with pressure or potential density anomaly, as specified by the
#' `profileControl` argument.
#'
#' * For `which="QC"`, two time-series panels are shown, with
#' time being that recorded in the individual profile in the dataset.
#' An additional argument named `parameter` must be given, to name the
#' quantity of interest. The function only works if `x` is an
#' [`argoFloats-class`] object created with [readProfiles()].
#' The top panel shows the percent of data flagged with codes
#' 1 (meaning good data), 2 (probably good), 5 (changed)
#' or 8 (estimated), as a function of time (lower axis) and
#' (if all cycles are from a single Argo float)
#' cycle number (upper axis, with smaller font).
#' Thus, low values on the top panel reveal
#' profiles that are questionable. Note that if all of data at a given time
#' have flag 0, meaning not assessed, then a quality of 0 is plotted at that
#' time. The bottom panel shows the mean value of the parameter in question
#' regardless of the flag value.
#'
#' * For `which="summary"`, one or more time-series panels are shown
#' in a vertical stack. If there is only one ID in `x`, then the cycle
#' values are indicated along the top axis of the top panel. The choice
#' of panels is set by the `summaryControl` argument.
#'
#' * For `which="TS"`, a TS plot is created, by calling [plotArgoTS()]
#' with the specified `x`, `xlim`, `ylim`, `xlab`, `ylab`,
#' `type`, `cex`, `col`, `pch`, `bg`, `eos`, and `TSControl`,
#' along with `debug-1`. See the help for [plotArgoTS()]
#' for how these parameters are interpreted.
#'
#' @param x an [`argoFloats-class`] object.
#'
#' @param which a character value indicating the type of plot. The possible
#' choices are `"map"`, `"profile"`, `"QC"`, `"summary"` and `"TS"`;
#' see \dQuote{Details}.
#'
#' @param bathymetry an argument used only if `which="map"`, to control
#' whether (and how) to indicate water depth. Note that the default
#' was `TRUE` until 2021-12-02, but was changed to `FALSE` on that
#' date, to avoid a bathymetry download, which can be a slow operations.
#' See \dQuote{Details} for details,
#' and Example 4 for a hint on compensating for the margin
#' adjustment done if an image is used to show bathymetry.
#'
#' @param geographical flag indicating the style of axes
#' for the `which="map"` case, but only if no projection is called
#' for in the `mapControl` argument. With `geographical=0` (which
#' is the default), the axis ticks are labeled with signed
#' longitudes and latitudes, measured in degrees. The signs are dropped
#' with `geographical=1`.
#' In the `geographical=4` case, signs are also dropped, but hemispheres
#' are indicated by writing `S`, `N`, `W` or `E` after axis tick labels,
#' except at the equator and prime meridian.
#' Note that this scheme mimics that used by [oce::plot,coastline-method()].
#'
#' @param xlim,ylim numerical values, each a two-element vector, that
#' set the `x` and `y` limits of plot axes, as for [plot.default()] and
#' other conventional plotting functions.
#'
#' @param xlab a character value indicating the name for the horizontal axis, or
#' `NULL`, which indicates that this function should choose an appropriate name
#' depending on the value of `which`. Note that `xlab` is not obeyed if
#' `which="TS"`, because altering that label can be confusing to the user.
#'
#' @param ylab as `xlab`, but for the vertical axis.
#'
#' @param type a character value that controls the line type, with `"p"` for
#' unconnected points, `"l"` for line segments between undrawn points, etc.;
#' see the docs for [par()], If `type` not specified, it defaults to `"p"`.
#'
#' @param cex a character expansion factor for plot symbols, or `NULL`, to get an
#' value that depends on the value of `which`.
#'
#' @param col the colour to be used for plot symbols, or `NULL`, to get an value
#' that depends on the value of `which` (see \dQuote{Details}). If `which="TS"`,
#' then the `TSControl` argument takes precedence over `col`.
#'
#' @param bg the colour to be used for plot symbol interior, for `pch`
#' values that distinguish between the interior of the symbol and the
#' border, e.g. for `pch=21`.
#'
#' @param pch an integer or character value indicating the type of plot symbol,
#' or `NULL`, to get a value that depends on the value of `which`.
#' (See [par()] for more on specifying `pch`.)
#'
#' @param eos a character value indicating the equation of state to use
#' if `which="TS"`. This must be `"gsw"` (the default) or `"unesco"`;
#' see [oce::plotTS()].
#'
#' @param mapControl a list that permits particular control of the `which="map"`
#' case. If provided, it may contain elements named `bathymetry` (which
#' has the same effect as the parameter `bathymetry`), `colLand`
#' (which indicates the colour of the land), and `projection` (which
#' may be `FALSE`, meaning to plot longitude and latitude on rectilinear axes,
#' `TRUE`, meaning to plot with [oce::mapPlot()], using Mollweide projection that
#' is suitable mainly for world-scale views, or a character value that will be
#' supplied to [oce::mapPlot()]. If a projection is used, then the positions
#' of the Argo floats are plotted with [oce::mapPoints()], rather than with
#' [points()], and if the user wishes to locate points with mouse clicks,
#' then [oce::mapLocator()] must be used instead of [locator()]. If `bathymetry`
#' is not contained in `mapControl`, it defaults to `FALSE`, and if `projection`
#' is not supplied, it defaults to `FALSE`. Note that `mapControl` takes
#' precedence over the `bathymetry` argument, if both are provided.
#' Also note that, at present, bathymetry cannot be shown with map projections.
#'
#' @param profileControl a list that permits control of the `which="profile"`
#' case. If provided, it may contain elements named `parameter` (a character value
#' naming the quantity to plot on the x axis), `ytype` (a character value equal to
#' either `"pressure"` or `"sigma0"`) and `connect` (a logical value indicating
#' whether to skip across `NA` values if the `type` argument is `"l"`, `"o"`,
#' or `"b"`).
#' If `profileControl` is not provided, it defaults to
#' `list(parameter="temperature", ytype="pressure", connect=FALSE)`. Alternatively,
#' if `profileControl` is missing any of the three elements, then they are
#' given defaults as in the previous sentence.
#'
#' @param QCControl a list that permits control of the `which="QC"`
#' case. If provided, it may contain an element named `parameter`, a character
#' value naming the quantity for which the quality-control information is
#' to be plotted, and an element named `dataStateIndicator`, a logical
#' value controlling whether to add a panel showing this quantity
#' (see Reference Table 6 of reference 1 to learn more
#' about what is encoded in `dataStateIndicator`).
#' If not provided, `QCControl` defaults to
#' `list(parameter="temperature",dataStateIndicator=FALSE)`.
#'
#' @param summaryControl a list that permits control of the `which="summary"`.
#' If provided, it should contain an element named `items`, a character vector
#' naming the items to be shown. The possible entries in this vector are
#' `"dataStateIndicator"` (see Reference Table 6 of reference 1,
#' for more information on this quantity), `"length"` (the number of levels
#' in the profile), `"deepest"` (the highest pressure recorded),
#' `"longitude"` and `"latitude"`. If `summaryControl` is not provided,
#' all of these will be shown. If all the elements of `x` have the same
#' `ID`, then the top panel will have ticks on its top axis, indicating
#' the `cycle`.
#'
#' @param TSControl a list that permits control of the `which="TS"`
#' case, and is ignored for the other cases. If `TSControl` is
#' supplied, then other arguments that normally control plots
#' are ignored. `TSControl` may have elements named
#' `drawByCycle`, `cex`, `col`, `pch`, `lty, `lwd`
#' and `type`. The first of these is a logical value indicating
#' whether to plot individual cycles with different aesthetic
#' characteristics, and the other elements describe those characteristics
#' in the usual way. By default, `type` is `"p"`, and the point
#' characteristics are given by the other parameters, which are passed to
#' [rep()] to have one value per cycle. For example,
#' `TSControl=list(drawByCycle=TRUE,col=1:2,type="l")` indicates
#' to draw lines for the individual cycles, with the first being
#' black, the second red, the third black, etc.
#' **FIXME: THIS IS NOT CODED YET.**
#'
#' @param debug an integer specifying the level of debugging.
#'
#' @param \dots extra arguments passed to the plot calls that are made
#' within this function.
#'
#' @examples
#' # Example 1: map profiles in index
#' library(argoFloats)
#' data(index)
#' plot(index)
#'
#' # Example 2: as Example 1, but narrow the margins and highlight floats
#' # within a circular region of diameter 100 km.
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0))
#' plot(index)
#' lon <- index[["longitude"]]
#' lat <- index[["latitude"]]
#' near <- oce::geodDist(lon, lat, -77.06, 26.54) < 100
#' R <- subset(index, near)
#' points(R[["longitude"]], R[["latitude"]],
#' cex=0.6, pch=20, col="red")
#' par(oldpar)
#'
#' # Example 3: TS of a built-in data file.
#' f <- system.file("extdata", "SR2902204_131.nc", package="argoFloats")
#' a <- readProfiles(f)
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(3.3, 3.3, 1, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
#' plot(a, which="TS")
#' par(oldpar)
#'\donttest{
#' # Example 4: map with bathymetry. Note that par(mar) is adjusted
#' # for the bathymetry palette, so it must be adjusted again after
#' # the plot(), in order to add new plot elements.
#' # This example specifies a coarse bathymetry dataset that is provided
#' # by the 'oce' package. In typical applications, the user will use
#' # a finer-scale dataset, either by using bathymetry=TRUE (which
#' # downloads a file appropriate for the plot view), or by using
#' # an already-downloaded file.
#' data(topoWorld, package="oce")
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 2), mgp=c(2, 0.7, 0)) # narrow margins for a map
#' plot(index, bathymetry=list(source=topoWorld))
#' # For bathymetry plots that use images, plot() temporarily
#' # adds 2.75 to par("mar")[4] so the same must be done, in order
#' # to correctly position additional points (shown as black rings).
#' par(mar=par("mar") + c(0, 0, 0, 2.75))
#' points(index[["longitude"]], index[["latitude"]],
#' cex=0.6, pch=20, col="red")
#' par(oldpar)
#'
#' # Example 5. Simple contour version, using coarse dataset (ok on basin-scale).
#' # Hint: use oce::download.topo() to download high-resolution datasets to
#' # use instead of topoWorld.
#' oldpar <- par(no.readonly=TRUE)
#' par(mar=c(2, 2, 1, 1))
#' data(topoWorld, package="oce")
#' plot(index, bathymetry=list(source=topoWorld, contour=TRUE))
#' par(oldpar)
##
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 5B. World view with Mollweide projection (Canada Day, 2020)
## jul1 <- subset(getIndex(), time=list(from="2020-09-01", to="2020-09-02"))
## par(mar=c(2, 2, 1, 1), mgp=c(2, 0.7, 0)) # narrow margins for a map
## plot(jul1, which="map", mapControl=list(projection=TRUE), bathymetry=FALSE,
## pch=20, col=4, cex=0.75)
#'
#' # Example 6. Customized map, sidestepping plot,argoFloats-method().
#' lon <- topoWorld[["longitude"]]
#' lat <- topoWorld[["latitude"]]
#' asp <- 1/cos(pi/180*mean(lat))
#' # Limit plot region to float region.
#' xlim <- range(index[["longitude"]])
#' ylim <- range(index[["latitude"]])
#' # Colourize 1km, 2km, etc, isobaths.
#' contour(x=lon, y=lat, z=topoWorld[["z"]], xlab="", ylab="",
#' xlim=xlim, ylim=ylim, asp=asp,
#' col=1:6, lwd=2, levels=-1000*1:6, drawlabels=FALSE)
#' # Show land
#' data(coastlineWorldFine, package="ocedata")
#' polygon(coastlineWorldFine[["longitude"]],
#' coastlineWorldFine[["latitude"]], col="lightgray")
#' # Indicate float positions.
#' points(index[["longitude"]], index[["latitude"]], pch=20)
#'
## # Example 7: TS plot for a particular Argo
## a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
## par(mar=c(3.3, 3.3, 1, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
## plot(a, which="TS")
##
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 7: Temperature QC plot for 25 cycles of a float in the Arabian Sea
## ais <- getIndex(filename="synthetic")
## sub <- subset(subset(ais, ID='2902123'), 50:75)
## profiles <- getProfiles(sub)
## argos <- readProfiles(profiles)
## par(mar=c(3.3, 3.3, 2, 1), mgp=c(2, 0.7, 0)) # wide margins for axis names
## plot(argos, which="QC") # defaults to temperature
## plot(argos, which="QC", QCControl=list(parameter="salinity"))
## plot(argos, which="QC", QCControl=list(parameter="salinity",dataStateIndicator=TRUE))
#'
#' # Example 7: Temperature profile of the 131st cycle of float with ID 2902204
#' a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
#' oldpar <- par(no.readonly=TRUE)
#' par(mfrow=c(1, 1))
#' par(mgp=c(2, 0.7, 0)) # mimic the oce::plotProfile() default
#' par(mar=c(1,3.5,3.5,2)) # mimic the oce::plotProfile() default
#' plot(a, which="profile")
#' par(oldpar)
#'
#' # Example 8: As Example 7, but showing temperature dependence on potential density anomaly.
#' a <- readProfiles(system.file("extdata", "SR2902204_131.nc", package="argoFloats"))
#' oldpar <- par(no.readonly=TRUE)
#' par(mgp=c(2, 0.7, 0)) # mimic the oce::plotProfile() default
#' par(mar=c(1,3.5,3.5,2)) # mimic the oce::plotProfile() default
#' plot(a, which="profile", profileControl=list(parameter="temperature", ytype="sigma0"))
#' par(oldpar)
#'
## # NOTE: removed this for CRAN submission, because it's very slow, not
## # especially informative, and would require using tempdir() for saving
## # the index, which is *not* what we want to illustrate since it goes
## # against the whole point of caching.
## # Example 10: Summary plot
## par(mar=c(2, 3.3, 2, 1))
## a <- readProfiles(getProfiles(subset(getIndex(), ID=1901584)))
## plot(a, which="summary")
#'}
#'
#' @references
#' 1. Carval, Thierry, Bob Keeley, Yasushi Takatsuki, Takashi Yoshida, Stephen Loch,
#' Claudia Schmid, and Roger Goldsmith. Argo User's Manual V3.3. Ifremer, 2019.
#' `doi:10.13155/29825`
#'
#' 2. Kelley, D. E., Harbin, J., & Richards, C. (2021). argoFloats: An R package for analyzing
#' Argo data. Frontiers in Marine Science, (8), 636922.
#' \doi{10.3389/fmars.2021.635922}
#'
#' @importFrom graphics abline axis box contour par plot.window points polygon rect text
#'
#' @importFrom grDevices extendrange gray rgb
#' @importFrom utils data
#' @importFrom oce as.ctd colormap drawPalette imagep oceColorsGebco oce.plot.ts plotTS
#'
#' @export
#'
#' @aliases plot,argoFloats-method
#'
#' @return None (invisible NULL).
#'
#' @author Dan Kelley and Jaimie Harbin
setMethod(f="plot",
signature=signature("argoFloats"),
definition=function(x, # plot,argoFloats-method
which="map",
bathymetry=argoDefaultBathymetry(),
geographical=0,
xlim=NULL, ylim=NULL,
xlab=NULL, ylab=NULL,
type="p", cex=par("cex"), col=par("fg"), pch=par("pch"), bg=par("bg"),
eos="gsw",
mapControl=NULL,
profileControl=NULL,
QCControl=NULL,
summaryControl=NULL,
TSControl=NULL,
debug=0,
...)
{
if (!requireNamespace("oce", quietly=TRUE))
stop("must install.packages(\"oce\") for plot() to work")
debug <- if (debug > 3) 3 else max(0, floor(debug + 0.5))
argoFloatsDebug(debug, "plot(x, which=\"", which, "\") {\n", sep="", unindent=1, style="bold")
dots <- list(...)
if (!inherits(x, "argoFloats"))
stop("In plot() : method is only for objects of class \"argoFloats\"", call.=FALSE)
if (length(which) != 1)
stop("\"which\" must contain only one item")
if (!geographical %in% c(0, 1, 4))
stop("In plot() : \"geographical\" must be 0, 1, or 4", call.=FALSE)
istraj <- identical(x@metadata$subtype, "trajectories")
if (which == "map") {
data("coastlineWorld", package="oce", envir=environment())
argoFloatsDebug(debug, "map plot\n", sep="")
if (!istraj) {
longitude <- x[["longitude", debug=debug]]
latitude <- x[["latitude", debug=debug]]
}
if (istraj) {
lat1 <- as.numeric(x[["latitude_min"]])
lat2 <- as.numeric(x[["latitude_max"]])
lon1 <- as.numeric(x[["longitude_min"]])
lon2 <- as.numeric(x[["longitude_max"]])
latitude <- c(lat1, lat2)
longitude <- c(lon1, lon2)
}
n <- x[["length"]]
## Find type of each cycle, for colour-coding
cycleType <- rep("core", n)
cycleType[("849" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[("862" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[("864" == x@data$index$profiler_cycleType)] <- "deep"
cycleType[!is.na(x@data$index$parameters)] <- "Bgc"
colLand <- mapControl$colLand
if (is.null(mapControl))
mapControl <- list(bathymetry=bathymetry, projection=FALSE, colLand="lightgray")
if (!"bathymetry" %in% names(mapControl))
mapControl$bathymetry <- TRUE
if (!"projection" %in% names(mapControl))
mapControl$projection <- FALSE
if (is.logical(mapControl$projection)) {
mapControl$projection <- if (mapControl$projection) "+proj=moll" else "none"
}
if (substr(mapControl$projection, 1, 6) != "+proj=" && mapControl$projection != "none")
stop("In plot,argoFloats-method(): mapControl$projection must start with \"+proj=\"", call.=FALSE)
if (mapControl$projection != "none") {
data("coastlineWorld", package="oce", envir=environment())
coastlineWorld <- get("coastlineWorld")
oce::mapPlot(coastlineWorld, col= "lightgray", projection=mapControl$projection, drawBox=FALSE)
if (!istraj) {
oce::mapPoints(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
col=if (is.null(col)) "white" else col,
pch=if (is.null(pch)) 21 else pch,
bg=if (is.null(bg)) "red" else bg,
type=if (is.null(type)) "p" else type,
...)
} else if (istraj) {
rect(lon1,lat1, lon2,lat2)
}
## warning("In plot,argoFloats-method(): projected maps do not (yet) show bathymetry", call.=FALSE)
return(invisible(NULL))
}
if (is.null(xlim))
if (!istraj) {
xlim <- extendrange(longitude)
} else if (istraj) {
xlim <- extendrange(range(lon1,lon2))
}
if (is.null(ylim))
if (!istraj) {
ylim <- extendrange(latitude)
} else if (istraj) {
ylim <- extendrange(range(lat1, lat2))
}
## Draw empty plot box, with axes, to set par("usr") for later use with bathymetry.
xlab <- if (is.null(xlab)) "" else xlab
ylab <- if (is.null(ylab)) "" else ylab
## Decode bathymetry
if (is.logical(mapControl$bathymetry)) {
drawBathymetry <- mapControl$bathymetry
bathymetry <- list(source="auto", contour=FALSE, colormap="auto", palette=TRUE)
} else if (is.list(mapControl$bathymetry)) {
drawBathymetry <- TRUE
if (!("source" %in% names(mapControl$bathymetry)))
stop("In plot() : \"bathymetry\" is a list, it must contain \"source\", at least", call.=FALSE)
if ("keep" %in% names(bathymetry))
warning("bathymetry$keep (an old argument) is no longer used\n")
if (is.null(bathymetry$contour))
bathymetry$contour <- FALSE
if (is.null(bathymetry$colormap))
bathymetry$colormap <- "auto"
if (is.null(bathymetry$palette))
bathymetry$palette <- TRUE
} else {
stop("In plot() : \"bathymetry\" must be either a logical or a list value", call.=FALSE)
}
if (!is.logical(bathymetry$palette))
stop("In plot() : \"bathymetry$palette\" must be a logical value", call.=FALSE)
argoFloatsDebug(debug, "drawBathymetry calculated to be ", drawBathymetry, "\n", sep="")
asp <- 1 / cos(pi/180*mean(range(latitude, na.rm=TRUE)))
argoFloatsDebug(debug, "asp=", asp, "\n", sep="")
if (drawBathymetry) {
argoFloatsDebug(debug, "handling bathymetry\n", sep="")
if (!requireNamespace("oce", quietly=TRUE))
stop("must install.packages(\"oce\") to plot with bathymetry")
## Handle bathymetry file downloading (or the use of a supplied value)
bathy <- NULL
if (is.character(bathymetry$source) && bathymetry$source == "auto") {
argoFloatsDebug(debug, "must either download bathymetry or use existing file\n", sep="")
if (!bathymetry$contour && bathymetry$palette) {
omar <- par("mar")
par(mar=omar + c(0, 0, 0, 2.75))
on.exit(par(mar=omar))
argoFloatsDebug(debug, "temporarily set par(mar=c(", paste(par("mar"), collapse=", "), ")) to allow for the palette\n", sep="")
}
argoFloatsDebug(debug, "using plot.window() to determine area for bathymetry download, with\n",
" extendrange(longitude)=c(", paste(extendrange(longitude), collapse=","), ")\n",
" extendrange(latitude)=c(", paste(extendrange(latitude), collapse=","), ")\n",
" xlim=c(", paste(xlim, collapse=","), ")\n",
" ylim=c(", paste(ylim, collapse=","), ")\n",
" asp=", asp, "\n", sep="")
plot.window(xlim=xlim, ylim=ylim,
xaxs="i", yaxs="i",
asp=asp,
xlab=xlab, ylab=ylab)
usr <- par("usr")
argoFloatsDebug(debug, "after using plot.window(), usr=c(", paste(round(usr, 4), collapse=", "), ")\n", sep="")
latitudeSpan <- usr[4] - usr[3]
longitudeSpan <- usr[2] - usr[1]
Dlon <- 0.05 * longitudeSpan
Dlat <- 0.05 * latitudeSpan
argoFloatsDebug(debug, "Dlon=", Dlon, " (5% of longitude span, which is ", longitudeSpan, ")\n", sep="")
argoFloatsDebug(debug, "Dlat=", Dlat, " (5% of latitude span, which is ", latitudeSpan, ")\n", sep="")
resolution <- as.integer(round(1 + 60 * latitudeSpan / 400))
if (resolution < 1)
stop("calculated resolution (=", resolution, ") cannot be under 1 minute")
argoFloatsDebug(debug, "resolution=", resolution, "\n", sep="")
## Notice the use of 0.1 and 1 digit in lon and lat, which means
## we are enlarging the field by 10km in both lon and lat.
west <- round(usr[1]-Dlon, 1)
east <- round(usr[2]+Dlon, 1)
south <- round(usr[3]-Dlat, 1)
north <- round(usr[4]+Dlat, 1)
argoFloatsDebug(debug, "will get topo in domain west=", west, ", east=", east,
", south=", south, ", north=", north, "\n")
topo <- try(oce::read.topo(oce::download.topo(west=west,
east=east,
south=south,
north=north,
resolution=resolution,
debug=debug)),
silent=FALSE)
if (inherits(topo, "try-error")) {
warning("could not download bathymetry from NOAA server\n")
drawBathymetry <- FALSE
} else {
bathy <- -topo[["z"]]
dimBathy <- dim(bathy)
argoFloatsDebug(debug, "grid size", paste(dimBathy, collapse="x"), "\n")
bathy <- as.integer(bathy)
dim(bathy) <- dimBathy
rownames(bathy) <- topo[["longitude"]]
colnames(bathy) <- topo[["latitude"]]
class(bathy) <- "bathy"
}
} else if (inherits(bathymetry$source, "bathy")) {
argoFloatsDebug(debug, "using supplied bathymetry$source\n", sep="")
bathy <- bathymetry$source
} else if (inherits(bathymetry$source, "oce") && inherits(bathymetry$source, "topo")) {
argoFloatsDebug(debug, "using oce-style topo object (converted to NOAA bathy)\n", sep="")
dim <- dim(bathymetry$source[["z"]])
## Note the negative sign, to get depth.
bathy <- matrix(-as.integer(bathymetry$source[["z"]]), nrow=dim[1], ncol=dim[2])
rownames(bathy) <- bathymetry$source[["longitude"]]
colnames(bathy) <- bathymetry$source[["latitude"]]
class(bathy) <- "bathy"
} else {
stop("cannot determine bathymetry data source")
}
argoFloatsDebug(debug, "drawBathymetry = ", drawBathymetry, " (after trying to get bathymetry)\n", sep="")
## Handle colormap
if (drawBathymetry && !inherits(bathy, "try-error") && !bathymetry$contour && is.character(bathymetry$colormap) && length(bathymetry$colormap) == 1 && bathymetry$colormap == "auto") {
argoFloatsDebug(debug, "bathy class: \"", class(bathy), "\"\n", sep="")
deepest <- max(bathy, na.rm=TRUE)
argoFloatsDebug(debug, "bathy range: ", paste(range(bathy, na.rm=TRUE), collapse=" to "), "\n")
argoFloatsDebug(debug, "setting a default colormap for ", deepest, "m to 0m depth\n", sep="")
bathymetry$colormap <- oce::colormap(zlim=c(0, deepest),
col=function(n) rev(oce::oceColorsGebco(n)))
}
## Handle optional drawing of the palette. Note space for one line of text is removed
## from the LHS of the palette and added to the RHS. This is because we know that
## there is no axis to the left of the palette, so we do not need space for one.
## We recover the stolen space by putting it back at the RHS, where it can be
## useful, especially if a map plot has another plot to its right.
if (drawBathymetry && !inherits(bathy, "try-error") && !bathymetry$contour && bathymetry$palette) {
argoFloatsDebug(debug, "drawing a bathymetry palette\n")
omai <- par("mai") # save this, since oce::drawPalette() may change it
oce::drawPalette(colormap=bathymetry$colormap, cex=if (is.null(cex)) 1 else par("cex.axis"))
on.exit(par(mai=omai))
} else {
argoFloatsDebug(debug, "not drawing a bathymetry palette, as instructed (or failed bathymetry download)\n")
}
}
argoFloatsDebug(debug, "initial xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
## Extend to world edge, if we are close and scale is large
if (diff(xlim) > 340) {
if (xlim[1] < -170) xlim[1] <- -180
if (xlim[2] > 170) xlim[2] <- 180
}
if (diff(ylim) > 160) {
if (ylim[1] < -70) ylim[1] <- -90
if (ylim[2] > 70) ylim[2] <- 90
}
argoFloatsDebug(debug, "after to-world-edge, xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
## Trim to world domain
xlim <- ifelse(xlim < -180, -180, ifelse(xlim > 180, 180, xlim))
ylim <- ifelse(ylim < -80, -90, ifelse(ylim > 80, 90, ylim))
argoFloatsDebug(debug, "after domain-trim, xlim=",
"c(", paste(xlim, collapse=","), ") and ylim=",
"c(", paste(ylim, collapse=","), ")\n", sep="")
argoFloatsDebug(debug, "about to plot, with\n",
" extendrange(longitude)=c(", paste(extendrange(longitude), collapse=","), ")\n",
" extendrange(latitude)=c(", paste(extendrange(latitude), collapse=","), ")\n",
" xlim=c(", paste(xlim, collapse=","), ")\n",
" ylim=c(", paste(ylim, collapse=","), ")\n",
" asp=", asp, "\n", sep="")
plot(xlim, ylim,
xaxs="i", yaxs="i",
asp=asp,
xlab=xlab, ylab=ylab, type="n", axes=FALSE)
## Draw box and axes. Note that we cannot just use axis() because we
## do not want to show e.g. lat beyond 90N.
usrTrimmed <- par("usr")
usrTrimmed[1] <- max(-180, usrTrimmed[1])
usrTrimmed[2] <- min(180, usrTrimmed[2])
usrTrimmed[3] <- max(-90, usrTrimmed[3])
usrTrimmed[4] <- min(90, usrTrimmed[4])
rect(usrTrimmed[1], usrTrimmed[3], usrTrimmed[2], usrTrimmed[4], lwd=1)
## Construct axes
xat <- pretty(usrTrimmed[1:2])
if (diff(range(xat)) > 300)
xat <- seq(-180, 180, 45)
yat <- pretty(usrTrimmed[3:4])
if (diff(range(yat)) > 150)
yat <- seq(-90, 90, 45)
if (geographical == 0) {
xlabels <- xat
ylabels <- yat
} else if (geographical == 1) {
xlabels <- abs(xat)
ylabels <- abs(yat)
} else if (geographical == 4) {
xlabels <- paste(abs(xat), ifelse(xat < 0, "W", ifelse(xat > 0, "E", "")), sep="")
ylabels <- paste(abs(yat), ifelse(yat < 0, "S", ifelse(yat > 0, "N", "")), sep="")
} else {
stop("In plot,argoFloats-method() : programming error: \"geographical\" must be 0, 1, or 4", call.=FALSE)
}
axis(1, at=xat, labels=xlabels, pos=usrTrimmed[3])
axis(2, at=yat, labels=ylabels, pos=usrTrimmed[1])
axis(3, at=xat, labels=FALSE, pos=usrTrimmed[4])
axis(4, at=yat, labels=FALSE, pos=usrTrimmed[2])
## argoFloatsMapAxes(geographical=geographical)
if (drawBathymetry) {
if (bathymetry$contour) {
argoFloatsDebug(debug, "indicating bathymetry with contours\n")
contour(as.numeric(rownames(bathy)),
as.numeric(colnames(bathy)),
bathy,
levels=c(100, 200, 500, seq(1e3, 10e3, 1e3)),
labels=c("100m", "200m", "500m", paste(1:10, "km", sep="")),
labcex=0.9, add=TRUE)
} else {
argoFloatsDebug(debug, "indicating bathymetry with an image\n")
oce::imagep(as.numeric(rownames(bathy)),
as.numeric(colnames(bathy)),
bathy, colormap=bathymetry$colormap, add=TRUE)
}
rect(usrTrimmed[1], usrTrimmed[3], usrTrimmed[2], usrTrimmed[4], lwd=1)
}
if (is.null(col)) {
argoFloatsDebug(debug, "defaulting col to indicate cycle type\n")
col <- rep("", n)
col[cycleType == "core"] <- colDefaults$core
col[cycleType == "Bgc"] <- colDefaults$bgc
col[cycleType == "deep"] <- colDefaults$deep
}
## Select coastline. Unlike in oce::plot,coastline-method, we base our choice
## on just the distance spanned in the north-south direction.
ocedataIsInstalled <- requireNamespace("ocedata", quietly=TRUE)
if (ocedataIsInstalled) {
usr <- par("usr")
mapSpan <- diff(range(latitude, na.rm=TRUE)) * 111 # km
##? l <- oce::geodDist(pinlon(usr[1]), pinlat(usr[3]), pinlon(usr[1]), pinlat(usr[4]))
##? r <- oce::geodDist(pinlon(usr[2]), pinlat(usr[3]), pinlon(usr[2]), pinlat(usr[4]))
##? b <- oce::geodDist(pinlon(usr[1]), pinlat(usr[1]), pinlon(usr[2]), pinlat(usr[1]))
##? t <- oce::geodDist(pinlon(usr[1]), pinlat(usr[4]), pinlon(usr[2]), pinlat(usr[4]))
##? mapSpan <- max(l, r, b, t) # largest length [km]
C <- 2 * 3.14 * 6.4e3 # circumferance of earth [km]
argoFloatsDebug(debug, "mapSpan=", mapSpan, ", C=", C, "\n")
if (mapSpan < 500) {
argoFloatsDebug(debug, "using coastlineWorldFine from ocedata package\n")
data("coastlineWorldFine", package="ocedata", envir=environment())
coastlineWorldFine <- get("coastlineWorldFine")
polygon(coastlineWorldFine[["longitude"]], coastlineWorldFine[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
} else if (mapSpan < C / 4) {
argoFloatsDebug(debug, "using coastlineWorldMedium from ocedata package\n")
data("coastlineWorldMedium", package="ocedata", envir=environment())
coastlineWorldMedium <- get("coastlineWorldMedium")
polygon(coastlineWorldMedium[["longitude"]], coastlineWorldMedium[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
} else {
argoFloatsDebug(debug, "using coastlineWorld from oce package, since the span is large\n")
data("coastlineWorld", package="oce", envir=environment())
coastlineWorld <- get("coastlineWorld")
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col= if (is.null(colLand)) "lightgray" else colLand)
}
if (is.null(pch))
pch <- 21
type <- if (is.null(type)) "p" else type
if (length(pch) == 1 && pch == 21 && !istraj) {
if (type != "l") {
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= type,
...)
} else if (type == "l") {
index1 <- subset(x, ID=unique(x[["ID"]][1]), silent=TRUE)
longitude <- index1[["longitude"]]
latitude <- index1[["latitude"]]
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= "l",
...)
for (i in unique(x[["ID"]])[2:length(unique(x[["ID"]]))]) {
index2 <- subset(x, ID=i, silent=TRUE)
points(unlist(index2[["longitude"]]), unlist(index2[["latitude"]]),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
bg=col,
type= "l",
...)
}
}
} else if (!istraj) {
points(unlist(longitude), unlist(latitude),
cex=if (is.null(cex)) 1 else cex,
pch=pch,
col=col,
type=if (is.null(type)) "p" else type,
...)
} else if (istraj) {
rect(lon1,lat1, lon2,lat2)
}
} else {
argoFloatsDebug(debug, "using coastlineWorld from oce package, since the ocedata package is not installed\n")
data("coastlineWorld", package="ocedata", envir=environment())
coastlineWorld <- get("coastlineWorld")
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col=col)
}
} else if (which == "summary" && !istraj) {
argoFloatsDebug(debug, "summary plot\n", sep="")
if (x[["length"]] < 2)
stop("In plot,argoFloats-method() : cannot draw a summary plot with only one float cycle", call.=FALSE)
if (!all(names(dots) %in% c("summaryControl")))
stop(names(dots)," is not a valid argument for this plot type. Try summaryControl instead.")
if (is.null(summaryControl))
summaryControl <- list(items=c("dataStateIndicator", "longitude", "latitude", "length", "deepest"))
if (!"items" %in% names(summaryControl))
stop("summaryControl must be a list containing a character vector named \"items\"")
if (!(summaryControl$items) %in% c("dataStateIndicator", "longitude", "latitude", "length", "deepest"))
stop("The names of the items must be dataStateIndicator, longitude, latitude, length, or deepest, not ", paste(summaryControl$items, collapse=" "))
items <- summaryControl$items
nitems <- length(items)
if (nitems) {
omfrow <- par("mfrow")
omar <- par("mar")
omgp <- par("mgp")
par(mar=c(3,3,2,1.5), mgp=c(1.75, 0.5, 0))
if (nitems < 3)
par(mfrow=c(nitems, 1))
else if (nitems == 4)
par(mfrow=c(2, 2))
else
par(mfrow=c(3, 2))
on.exit(par(mfrow=omfrow, mar=omar, mgp=omgp))
time <- as.POSIXct(unlist(x[["time"]]), origin="1970-01-01", tz="UTC")
o <- order(time)
for (iitem in seq_len(nitems)) {
argoFloatsDebug(debug, "items[", iitem, "]=\"", items[iitem], "\"\n", sep="")
if (items[iitem] == "dataStateIndicator") {
y <- unlist(x[["dataStateIndicator"]])
argoFloatsDebug(debug, oce::vectorShow(y))
if (length(y)) {
u <- sort(unique(y))
yy <- seq_along(u)
oce::oce.plot.ts(range(time), range(yy), ylab="Data State Ind.",
drawTimeRange=FALSE, type="n", axes=FALSE,
xaxs="i", ylim=range(0.8, (max(yy)+0.2)),
mar=par("mar"), mgp=par("mgp"))
abline(h=seq_along(u), col="gray")
points(time[o], factor(y)[o])
oce::oce.axis.POSIXct(side=1, drawTimeRange=FALSE)
box()
axis(side=2, at=yy, labels=u)
}
} else if (items[iitem] == "length") {
y <- sapply(x[["argos"]], function(a) length(a[["pressure"]]))
oce::oce.plot.ts(time[o], y[o], ylab="Length", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "longitude") {
y <- sapply(x[["longitude"]], function(a) a[1])
oce::oce.plot.ts(time[o], y[o], ylab="Longitude", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "latitude") {
y <- sapply(x[["latitude"]], function(a) a[1])
oce::oce.plot.ts(time[o], y[o], ylab="Latitude", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
} else if (items[iitem] == "deepest") {
y <- sapply(x[["argos"]], function(a) max(a[["pressure"]], na.rm=TRUE))
oce::oce.plot.ts(time[o], y[o], ylab="Max Pres.", drawTimeRange=FALSE,
type="p", xaxs="i",
mar=par("mar"), mgp=par("mgp"))
}
if (iitem == 1) {
if (1 == length(unique(x[["ID"]])))
axis(side=3, at=time[o], labels=x[["cycle"]][o],
cex.axis=max(0.75, 0.75*par("cex")))
}
}
}
} else if (which == "summary" && istraj) {
stop("In plot,argoFloats-method(): cannot draw summary plot for subtype = trajectories", call.=FALSE)
} else if (which == "TS") {
#cat("next is cex before call:\n");cat(vectorShow(cex))
#message("main debug=", debug)
plotArgoTS(x=x, xlim=xlim, ylim=ylim,
type=type, cex=cex, col=col, pch=pch, bg=bg,
eos=eos, TSControl=TSControl, debug=debug-1L)
argoFloatsDebug(debug, "} # plot()\n", sep="", unindent=1, style="bold")
return()
#<20230805> #return(do.call("plotArgoTS", args))
#<20230805> stop("FIXME: cannot get here")
#<20230805> argoFloatsDebug(debug, "TS plot\n", sep="")
#<20230805> if ((x[["type"]] != "argos"))
#<20230805> stop("In plot,argoFloats-method(): x must have been created by readProfiles()", call.=FALSE)
#<20230805> if (!(eos %in% c("gsw", "unesco")))
#<20230805> stop("In plot,argoFloats-method(): eos must be \"gsw\" or \"unesco\", not \"", eos, "\"", call.=FALSE)
#<20230805> # compute cycle, in case we need that
#<20230805> cycleIndex <- unlist(lapply(seq_len(x[["length"]]),
#<20230805> function(icycle)
#<20230805> rep(icycle, length(x@data$argos[[icycle]]@data$pressure))))
#<20230805> salinity <- unlist(x[["salinity", debug=debug-1L]])
#<20230805> temperature <- unlist(x[["temperature", debug=debug-1L]])
#<20230805> pressure <- unlist(x[["pressure", debug=debug-1L]])
#<20230805> ## Use byLevel to repeat the latitude and longitude values across
#<20230805> ## the depths in each profile, so that the resultant vector
#<20230805> ## will match the salinity, temperature and pressure vectors.
#<20230805> latitude <- unlist(x[["latitude", "byLevel", debug=debug-1L]])
#<20230805> longitude <- unlist(x[["longitude", "byLevel", debug=debug-1L]])
#<20230805> # Interpolate across NA longitudes (required for traj data, to get TS plot)
#<20230805> n <- length(longitude)
#<20230805> if (any(is.na(longitude)))
#<20230805> longitude <- approx(1:n, longitude, 1:n)$y
#<20230805> if (any(is.na(latitude)))
#<20230805> latitude <- approx(1:n, latitude, 1:n)$y
#<20230805> ctd <- oce::as.ctd(salinity=salinity,
#<20230805> temperature=temperature,
#<20230805> pressure=pressure,
#<20230805> latitude=latitude,
#<20230805> longitude=longitude)
#<20230805> ctd@data$cycleIndex <- cycleIndex
#<20230805> if (is.null(TSControl)) {
#<20230805> argoFloatsDebug(debug, "defaulting TSControl\n")
#<20230805> TSControl <- list(colByCycle=NULL)
#<20230805> }
#<20230805> if (is.null(TSControl$colByCycle)) {
#<20230805> argoFloatsDebug(debug, "TSControl does not contain colByCycle\n")
#<20230805> if (is.null(col)) {
#<20230805> if (which == "TS") {
#<20230805> col <- "flags"
#<20230805> } else {
#<20230805> col <- rgb(0, 0, 1, 0.5)
#<20230805> }
#<20230805> }
#<20230805> } else {
#<20230805> argoFloatsDebug(debug, "TSControl contains colByCycle\n")
#<20230805> # Ignore "col" if TSControl contains "colByCycle"
#<20230805> colByCycle <- TSControl$colByCycle
#<20230805> cycle <- unlist(x[["cycle", debug=debug-1L]])
#<20230805> lengths <- sapply(x[["argos"]], function(cc) length(cc[["pressure"]]))
#<20230805> # Increase the col length, so e.g. TSControl=list(colByCycle=1:2) will alternate colours
#<20230805> colByCycle <- rep(colByCycle, length.out=length(cycle))
#<20230805> col <- unlist(lapply(seq_along(cycle), function(i) rep(colByCycle[i], lengths[i])))
#<20230805> }
#<20230805> if (is.null(cex)) {
#<20230805> cex <- par("cex")
#<20230805> argoFloatsDebug(debug, "defaulting to cex=", cex, "\n", sep="")
#<20230805> }
#<20230805> if (is.null(pch)) {
#<20230805> pch <- 20
#<20230805> argoFloatsDebug(debug, "defaulting to pch=", pch, "\n", sep="")
#<20230805> }
#<20230805> if (col[1] == "flags") {
#<20230805> argoFloatsDebug(debug, "col is \"flags\"\n")
#<20230805> salinityFlag <- unlist(x[["salinityFlag"]])
#<20230805> temperatureFlag <- unlist(x[["temperatureFlag"]])
#<20230805> goodT <- temperatureFlag %in% c(1, 2, 5, 8)
#<20230805> goodS <- salinityFlag %in% c(1, 2, 5, 8)
#<20230805> good <- goodS & goodT
#<20230805> okT <- temperatureFlag %in% c(0)
#<20230805> okS <- salinityFlag %in% c(0)
#<20230805> ok <- okS & okT
#<20230805> col <- ifelse(good, "black", ifelse(ok, "gray", "red"))
#<20230805> if (pch == 21)
#<20230805> bg <- ifelse(good, "black", ifelse(ok, "gray", "red"))
#<20230805> }
#<20230805> argoFloatsDebug(debug, "about to plotTS()\n")
#<20230805> if (!is.null(TSControl$colByCycle)) {
#<20230805> argoFloatsDebug(debug, "colorizing by index\n")
#<20230805> cycles <- unique(cycleIndex)
#<20230805> ncycles <- length(cycles)
#<20230805> if (eos == "gsw") {
#<20230805> Slim <- range(ctd[["SA"]], na.rm=TRUE)
#<20230805> Tlim <- range(ctd[["CT"]], na.rm=TRUE)
#<20230805> } else {
#<20230805> Slim <- range(ctd[["salinity"]], na.rm=TRUE)
#<20230805> Tlim <- range(ctd[["temperature"]], na.rm=TRUE)
#<20230805> }
#<20230805> argoFloatsDebug(debug, "set Slim=", Slim, "\n")
#<20230805> argoFloatsDebug(debug, "set Tlim=", Tlim, "\n")
#<20230805> if (length(longitude) != length(salinity))
#<20230805> longitude <- rep(longitude[1], length(salinity))
#<20230805> if (length(latitude) != length(salinity))
#<20230805> latitude <- rep(latitude[1], length(salinity))
#<20230805> col <- TSControl$colByCycle
#<20230805> if (length(col) != ncycles)
#<20230805> col <- rep(col, length.out=ncycles)
#<20230805> if (length(type) != ncycles)
#<20230805> type <- rep(type, length.out=ncycles)
#<20230805> if (length(pch) != ncycles)
#<20230805> pch <- rep(pch, length.out=ncycles)
#<20230805> for (i in seq_len(ncycles)) {
#<20230805> argoFloatsDebug(debug, " handling cycle ", i, " of ", ncycles, ", which has ",
#<20230805> length(salinity[i==cycleIndex]), " points\n", sep="")
#<20230805> ctd <- oce::as.ctd(
#<20230805> salinity=salinity[i == cycleIndex],
#<20230805> temperature=temperature[i == cycleIndex],
#<20230805> pressure=pressure[i == cycleIndex],
#<20230805> latitude=latitude[i == cycleIndex],
#<20230805> longitude=longitude[i == cycleIndex])
#<20230805> if (i == 1L) {
#<20230805> plotTS(ctd, Slim=Slim, Tlim=Tlim,
#<20230805> type=type[i], col=col[i], pch=pch[i])
#<20230805> } else {
#<20230805> plotTS(ctd, add=TRUE,
#<20230805> type=type[i], col=col[i], pch=pch[i])
#<20230805> }
#<20230805> }
#<20230805> } else {
#<20230805> argoFloatsDebug(debug, "not colorizing by index\n")
#<20230805> oce::plotTS(ctd,
#<20230805> cex=cex,
#<20230805> bg=bg,
#<20230805> col=col,
#<20230805> pch=pch,
#<20230805> mar=par("mar"), mgp=par("mgp"), eos=eos,
#<20230805> type=if (is.null(type)) "p" else type, ...)
#<20230805> }
} else if (which == "QC" && !istraj) {
if (x[["type"]] != "argos")
stop("In plot,argoFloats-method(): The type of x must be \"argos\"", call.=FALSE)
IDs <- x[["ID"]]
nID <- length(unique(IDs))
if (nID != 1)
stop("In plot,argoFloats-method(): It is only possible to plot a QC of a single ID", call.=FALSE)
knownParameters <- names(x[[1]]@metadata$flags) # FIXME: is it possible that later cycles have different flags?
if (!all(names(dots) %in% c("QCControl")))
stop(names(dots)," is not a valid argument for this plot type. Try QCControl instead.")
if (is.null(QCControl)) {
QCControl <- list(parameter="temperature")
}
if (!is.list(QCControl))
stop("In plot,argoFloats-method(): QCControl must be a list")
if (!all(names(QCControl) %in% c("parameter", "dataStateIndicator")))
stop("QCControl must contain only elements named \"parameter\" and \"dataStateIndicator\", but its elements are named \"", paste(names(QCControl), collapse="\", \""), "\"")
if (!"parameter" %in% names(QCControl))
QCControl$parameter <- "temperature"
if (!"dataStateIndicator" %in% names(QCControl))
QCControl$dataStateIndicator <- FALSE
if (!is.logical(QCControl$dataStateIndicator))
stop("QCControl element 'dataStateIndicator' must be logical")
if (length(QCControl) != 2)
stop("In plot,argoFloats-method(): QCControl must contain only elements named \"parameter\" and \"dataStateIndicator\"")
if (!(QCControl$parameter %in% knownParameters))
stop("In plot,argoFloats-method(): QCControl$parameter '", QCControl$parameter, "' not found. Try one of: \"", paste(sort(knownParameters), collapse="\", \""), "\"", call.=FALSE)
qf <- function(x) {
# qf returns 100 if data are all "good" = 1 or "probably good" = 2 or "changed" = 5 or "estimated" = 8
flag <- x[[paste0(QCControl$parameter, "Flag")]]
100 * sum(1 == flag | 2 == flag | 5 == flag | 8 == flag, na.rm=TRUE) / length(flag)
}
meanf <- function(x)
mean(x[[QCControl$parameter]], na.rm=TRUE)
time <- oce::numberAsPOSIXct(unlist(lapply(x[["argos"]], function(x) x[["time"]])))
q <- unlist(lapply(x[["argos"]], qf))
m <- unlist(lapply(x[["argos"]], meanf))
omfrow <- par("mfrow")
omar <- par("mar")
par(mar=c(3, 3, 2, 1))
par(mfrow=c(if (QCControl$dataStateIndicator) 3 else 2, 1))
on.exit(par(mfrow=omfrow, mar=omar))
if (any(is.finite(q))) {
o <- order(time) # cycles are not time-ordered in index files
oce::oce.plot.ts(time[o], q[o], ylab=paste(QCControl$parameter, "% Good"), drawTimeRange=FALSE, type="l", mar=par("mar"), xaxs="i", cex.lab=0.7, cex.axis=0.7)
points(time[o], q[o], col=ifelse(q[o] < 50, "red", "black"), pch=20, cex=1)
abline(h=50, col="red", lty="dashed")
if (1 == length(unique(x[["ID"]])))
axis(side=3, at=time[o], labels=x[["cycle"]][o], cex.axis=0.7, cex.lab=0.7)
if (QCControl$dataStateIndicator) {
y <- unlist(x[["dataStateIndicator"]])
if (length(y)) {
u <- sort(unique(y))
yy <- seq_along(u)
oce::oce.plot.ts(range(time), range(yy), ylab="Data State Ind.",
drawTimeRange=FALSE, type="n", axes=FALSE,
mar=par("mar"), mgp=par("mgp"),
xaxs="i", ylim=range(0.8, (max(yy)+0.2)))
abline(h=seq_along(u), col="gray")
points(time[o], factor(y)[o])
oce::oce.axis.POSIXct(side=1, drawTimeRange=FALSE)
box()
axis(side=2, at=yy, labels=u)
} else {
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0.5, 0.5, "x lacks dataStateIndicator")
}
}
oce::oce.plot.ts(time[o], m[o], ylab=paste(QCControl$parameter, "Mean"), drawTimeRange=FALSE, type="l",
mar=par("mar"), cex.lab=0.7, cex.axis=0.7)
points(time[o], m[o], col=ifelse(q[o] < 50, "red", "black"), pch=20, cex=1)
} else {
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0, 0.5, paste(' No', QCControl$parameter, "flags available"), pos=4)
plot(0:1, 0:1, xlab="", ylab='', type="n", axes=FALSE)
box()
text(0, 0.5, paste(' No', QCControl$parameter, "flags available"), pos=4)
}
} else if (which == "QC" && istraj) {
stop("In plot,argoFloats-method(): cannot draw QC plot for subtype = trajectories", call.=FALSE)
} else if (which == "profile") {
if (x[["type"]] != "argos")
stop("In plot,argoFloats-method(): The type of x must be \"argos\"", call.=FALSE)
if (is.null(type))
type <- "p"
if (is.null(profileControl)) {
if ("parameter" %in% names(dots)) {
warning("accepting \"parameter\" as a separate argument, until 2020-dec-01. After that, you must use e.g. profileControl=list(parameter=", dots$parameter, ")")
profileControl <- list(parameter=dots$parameter)
} else {
profileControl <- list(parameter="temperature")
}
profileControl$ytype <- "pressure"
profileControl$connect <- FALSE
}
if (!is.list(profileControl))
stop("In plot,argoFloats-method(): profileControl must be a list")
if (!"ytype" %in% names(profileControl))
profileControl$ytype <- "pressure"
if (!"parameter" %in% names(profileControl))
profileControl$parameter <- "temperature"
if (!"connect" %in% names(profileControl))
profileControl$connect <- TRUE
if (length(profileControl) != 3)
stop("In plot,argoFloats-method(): profileControl must contain only three elements, \"parameter\", \"ytype\" and \"connect\".")
if (!profileControl$ytype %in% c("pressure", "sigma0"))
stop("In plot,argoFloats-method(): profileControl$ytype must be \"pressure\" or \"sigma0\", not \"", profileControl$ytype, "\"")
N <- length(x[["argos"]])
## The known parameter names include not just the things stored in the
## data slot (of *any* of the profiles), but also some computable things. We
## will always have longitude and latitude, so there's no need to check for them.
## A core profile will always get the computable things, but bgc profiles
## may lack e.g. salinity and temperature, so all we need to check for is
## salinity, temperature, and pressure.
knownParameters <- unique(unlist(lapply(seq_len(N), function(i) names(x[[i]][["data"]]))))
if (all(c("salinity", "temperature", "pressure") %in% knownParameters)) {
knownParameters <- c(knownParameters,
"SA", "CT", "sigma0", "sigma1", "sigma2", "sigma3", "sigma4", "sigmaTheta",
"theta",
"spice",
"N2")
}
argoFloatsDebug(debug, "knownParameters: \"", paste(sort(knownParameters), collapse="\", \""), "\".\n", sep="")
if (!(profileControl$parameter %in% knownParameters))
stop("In plot,argoFloats-method(): profileControl$parameter=\"", profileControl$parameter, "\" is not in the dataset, or calculable from that dataset; try one of the following: \"", paste(sort(knownParameters), collapse='", "'), "\".", call.=FALSE)
y <- lapply(1:N, function(i) x[[i]][[profileControl$ytype]])
variable <- lapply(1:N, function(i) x[[i]][[profileControl$parameter]])
nn <- unlist(lapply(1:N, function(i) prod(dim(x[[i]][[profileControl$parameter]]))))
Y <- NULL
VARIABLE <- NULL
argoFloatsDebug(debug, "number of profiles: ", N, "\n")
for (i in seq_len(N)) {
if (nn[i] > 0) {
Y <- c(Y, NA, y[[i]])
VARIABLE <- c(VARIABLE, NA, variable[[i]])
}
}
if (type %in% c("l", "o", "b") && profileControl$connect) {
ok <- is.finite(VARIABLE) & is.finite(Y)
VARIABLE <- VARIABLE[ok]
Y <- Y[ok]
}
if (is.null(cex)) {
cex <- par("cex")
argoFloatsDebug(debug, "Profile plot defaulting to cex=", cex, "\n")
}
if (is.null(pch)) {
pch <- 20
argoFloatsDebug(debug, "Profile plot defaulting to pch=", pch, "\n")
}
plot(VARIABLE, Y,
xlim=xlim,
ylim=if (is.null(ylim)) rev(range(Y, na.rm=TRUE)) else ylim,
axes=FALSE,
ylab="", xlab="", # draw axes later, in oceanographic 'profile' locations
cex=cex,
type=if(is.null(type)) "l" else type,
col=if (is.null(col)) par("col") else col,
pch=pch,
...)
box()
axis(2)
axis(3)
ylab <- oce::resizableLabel(if (profileControl$ytype == "pressure") "p" else "sigma0", axis="y")
mtext(ylab, side=2, line=par("mgp")[1], cex=par("cex"))
# Rename for oce::resizableLabel()
if (profileControl$parameter == "salinity")
profileControl$parameter <- "S"
xlab <- oce::resizableLabel(if (profileControl$parameter == "oxygen") "oxygen umol/kg" else profileControl$parameter, axis="x")
mtext(xlab, side=3, line=par("mgp")[1], cex=par("cex"))
} else {
stop("In plot,argoFloats-method():cannot handle which=\"", which, "\"; see ?\"plot,argoFloats-method\"", call.=FALSE)
}
argoFloatsDebug(debug, "} # plot()\n", sep="", unindent=1, style="bold")
invisible(NULL)
}
)
Try the argoFloats package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.