diff --git a/DESCRIPTION b/DESCRIPTION index b4bc837..0e3d12c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: rerddapXtracto Type: Package Title: Extracts Environmental Data from ERD's ERDDAP Web Service -Version: 0.2.0 +Version: 0.3.0 Authors@R: person("Roy", "Mendelssohn", email = "roy.mendelssohn@noaa.gov", role = c("aut","cre")) Description: The xtractomatic package contains three functions that access environmental data from ERD's ERDDAP service. The rxtracto function extracts @@ -15,8 +15,6 @@ Description: The xtractomatic package contains three functions that access to be extracted. URL: https://github.com/rmendels/rerddapXtracto BugReports: https://github.com/rmendels/rerddapXtracto/issues -Depends: - R (>= 3.3.0) License: CC0 LazyData: TRUE Imports: @@ -24,16 +22,12 @@ Imports: methods, ncdf4, parsedate, + plotdap, rerddap, sp, stats Suggests: - akima, - dplyr, - ggfortify, ggplot2, knitr, - lubridate, mapdata, - xts -RoxygenNote: 5.0.1 +RoxygenNote: 6.0.1 diff --git a/R/checkBounds.R b/R/checkBounds.R index f8a0726..8dc7916 100644 --- a/R/checkBounds.R +++ b/R/checkBounds.R @@ -20,14 +20,33 @@ checkBounds <- function(dataCoordList, dimargs) { dimLen <- length(names(dataCoordList)) for (i in (1:dimLen)) { cIndex <- which(names(dataCoordList)[i] == names(dimargs)) - if ((min(dimargs[[cIndex]]) < min(dataCoordList[[i]])) | (max(dimargs[[cIndex]]) > max(dataCoordList[[i]]))) { - print('dimension coordinate out of bounds') - print(paste0('dimension name: ', names(dimargs)[cIndex])) - print(paste('given coordinate bounds', dimargs[cIndex])) - returnCode <- 1 - print(paste('ERDDAP datasets bounds', min(dataCoordList[i]), max(dataCoordList[i]))) + if (names(dimargs)[cIndex] == 'time') { + min_dimargs <- min(as.numeric(dimargs[[cIndex]])) + max_dimargs <- max(as.numeric(dimargs[[cIndex]])) + temp_time <- parsedate::parse_iso_8601(dataCoordList[[i]]) + min_coord <- min( as.numeric(temp_time)) + max_coord <- max( as.numeric(temp_time)) + } else { + min_dimargs <- min(dimargs[[cIndex]]) + max_dimargs <- max(dimargs[[cIndex]]) + min_coord <- min(dimargs[[cIndex]]) + max_coord <- max(dimargs[[cIndex]]) + } + if ((min_dimargs < min_coord) | (max_dimargs > max_coord)) { + if (names(dimargs)[cIndex] == 'time') { + print('dimension coordinate out of bounds') + print(paste0('dimension name: ', names(dimargs)[cIndex])) + print(paste('given coordinate bounds', min_dimargs, max_dimargs)) + returnCode <- 1 + print(paste('ERDDAP datasets bounds', as.Date(min_coord), as.Date( max_coord))) + } else { + print('dimension coordinate out of bounds') + print(paste0('dimension name: ', names(dimargs)[cIndex])) + print(paste('given coordinate bounds', min_dimargs, max_dimargs)) + returnCode <- 1 + print(paste('ERDDAP datasets bounds', min_coord, max_coord)) + } } - } diff --git a/R/checkInput.R b/R/checkInput.R index d006bb3..5788b0d 100644 --- a/R/checkInput.R +++ b/R/checkInput.R @@ -49,8 +49,8 @@ checkInput <- function(dataInfo, parameter, urlbase, callDims) { } # check that the base url ends in / lenURL <- nchar(urlbase) - if (substr(urlbase, lenURL, lenURL) == '/') { - urlbase <- substr(urlbase, 1, (lenURL - 1)) + if (substr(urlbase, lenURL, lenURL) != '/') { + urlbase <- paste0(urlbase, '/') } # check that urlbase connects to an ERDDAP diff --git a/R/getFIleCoords.R b/R/getFIleCoords.R index 0022140..1cd4f47 100644 --- a/R/getFIleCoords.R +++ b/R/getFIleCoords.R @@ -16,7 +16,7 @@ getfileCoords <- function(datasetID, dataCoords, urlbase) { # to start do brute force way with for loop coordList <- list() for (i in 1:length(dataCoords)) { - myURL <- paste0(urlbase, '/griddap/', datasetID, '.csv?', dataCoords[i], '[0:1:last]') + myURL <- paste0(urlbase, 'griddap/', datasetID, '.csv?', dataCoords[i], '[0:1:last]') coordVals <- utils::read.csv(myURL, skip = 2, header = FALSE, stringsAsFactors = FALSE) coordVals <- coordVals[, 1] coordList[[i]] <- coordVals diff --git a/R/plotBox.R b/R/plotBox.R new file mode 100644 index 0000000..1e1f090 --- /dev/null +++ b/R/plotBox.R @@ -0,0 +1,84 @@ +#' plot result of xtracto_3D or rxtracto_3D +#' +#' \code{plotBox} is a function to plot the results from +#' rxtracto() and xtracto() +#' +#' @export +#' @param resp data frame returned from rxtracto() or xtracto() +#' @param plotColor the color to use in plot from rerddap +#' @param time a function to map multi-time to one, or else identity for animation +#' @param animate if multiple times, if TRUE will animate the maps +#' @param myFunc function of one argument to transform the data +#' @param maxpixels maximum numbe rof pixels to use in making the map - controls resolution +#' @return a plotdap plot +#' +#' @examples +#' tagData <- Marlintag38606 +#' xpos <- tagData$lon +#' ypos <- tagData$lat +#' tpos <- tagData$date +#' zpos <- rep(0., length(xpos)) +#' urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +#' swchlInfo <- rerddap::info('erdSWchla8day') +#' swchl <- rxtracto(swchlInfo, parameter = 'chlorophyll', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, xlen = .2, ylen = .2) +#' plotBox(xpos, ypos, swchl, plotColor = 'chlorophyll') + +plotBox <- function(resp, plotColor = 'viridis', time = NA, animate = FALSE, name = NA, myFunc = NA, maxpixels = 10000){ + require(rerddap) + require(plotdap) + if (!is.function(time)) { + time <- function(x) mean(x, na.rm = TRUE) + } + if (is.function(myFunc)) { + resp[[1]] <- myFunc(resp[[1]]) + } + if (!is.na(name)) { + names(resp)[1] <- name + } + paramName = names(resp)[1] + myStruct <- meltnc(resp) + myStruct <- structure( + myStruct, + class = c("griddap_nc", "nc", "data.frame") + ) + p <- plotdap::plotdap() + parameter1 <- as.formula(paste('~', paramName)) + myList <- list(p, myStruct, parameter1, plotColor, time, animate, maxpixels ) + names(myList) <- c('plot', 'grid', 'var', 'fill', 'time', 'animate', 'maxpixels') + myplot <- do.call(plotdap::add_griddap, myList) + myplot +} + +meltnc <- function(resp ){ + ## modified from rerddap::ncdf4_get + rows = length(resp[[1]]) + if (is.null(resp$time)) { + exout <- do.call("expand.grid", list(longitude = resp$longtiude, latitude = resp$latitude)) + meta <- dplyr::arrange_(exout, names(exout)[1]) + } else { + time <- as.character(resp$time) + time <- suppressWarnings(rep(time, each = rows/length(resp$time))) + lat <- rep(rep(resp$latitude, each = length(resp$longitude)), + length(resp$time)) + lon <- rep(rep(resp$longitude, times = length(resp$latitude)), + times = length(resp$time)) + meta <- data.frame(time, lat, lon, stringsAsFactors = FALSE) + } + + # make data.frame + df <- as.vector(resp[[1]]) + df <- data.frame(df) + names(df) <- names(resp)[1] + alldf <- if (NROW(meta) > 0) cbind(meta, df) else df + + # Fool plotdap that there is a summary + summary_time = list(vals = as.numeric(resp$time)) + summary_lons <- list(vals = resp$longitude) + summary_lats <- list(vals = resp$latitude) + dims <- list(time = summary_time, longitude = summary_lons, latitude = summary_lats) + summary <- list(dims) + names(summary) <- 'dim' + # output + list(summary = summary, data = alldf) +} + diff --git a/R/plotTrack.R b/R/plotTrack.R new file mode 100644 index 0000000..8cdaa8c --- /dev/null +++ b/R/plotTrack.R @@ -0,0 +1,53 @@ +#' plot result of xtracto or rxtracto +#' +#' \code{plotTrack} is a function to plot the results from +#' rxtracto() and xtracto() +#' @export +#' @param xcoord passed to rxtracto() or xtracto() +#' @param ycoord passed to rxtracto() or xtracto() +#' @param resp data frame returned from rxtracto() or xtracto() +#' @param plotColor the color to use in plot from rerddap +#' @param myFunc function of one argument to transform the data +#' @return a plotdap plot +#' +#' @examples +#' tagData <- Marlintag38606 +#' xpos <- tagData$lon +#' ypos <- tagData$lat +#' tpos <- tagData$date +#' zpos <- rep(0., length(xpos)) +#' urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +#' swchlInfo <- rerddap::info('erdSWchla8day') +#' swchl <- rxtracto(swchlInfo, parameter = 'chlorophyll', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, xlen = .2, ylen = .2) +#' plotTrack(xpos, ypos, swchl, plotColor = 'chlorophyll') + +plotTrack <- function(xcoord, ycoord, resp, plotColor = 'viridis', name = NA, myFunc = NA, shape = 20, size = .5){ + require(rerddap) + require(plotdap) + ind <- which(xcoord > 180) + xcoord[ind] <- xcoord[ind] - 360 + if (is.function(myFunc)) { + resp[[1]] <- myFunc(resp[[1]]) + } + myDataFrame = data.frame(xcoord, ycoord, resp[[1]]) + nameLen <- nchar(names(resp)) + if (is.na(name)) { + paramName <- substr(names(resp)[1], 6, nameLen) + }else{ + paramName = name + } + names(myDataFrame) <- c('longitude', 'latitude', paramName) +myStruct <- structure( + myDataFrame, + class = c("tabledap", "data.frame") + ) +p <- plotdap::plotdap() +paramName1 <- as.formula(paste('~', paramName)) +myList <- list(p, myStruct, paramName1, plotColor, shape, size) +names(myList) <- c('plot', 'table', 'var', 'color', 'shape', 'size') +#plotCmd <- paste0('add_tabledap(plotdap(), myStruct, ~', paramName, +# ', color = ', deparse(plotColor), ', shape =20, size= .5)') +#myPlot <- eval(parse(text = plotCmd)) +myPlot <- do.call(plotdap::add_tabledap, myList) +myPlot +} diff --git a/R/rxtracto.R b/R/rxtracto.R index 5be8286..3b1a835 100644 --- a/R/rxtracto.R +++ b/R/rxtracto.R @@ -35,7 +35,7 @@ #' \item column 11 = median absolute deviation of data within search radius #' } #' @examples -#' urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +#' urlbase <- 'https://upwell.pfeg.noaa.gov/erddap' #' dataInfo <- rerddap::info('erdMBsstd8day') #' parameter <- 'sst' #' xcoord <- c(230, 231) @@ -72,7 +72,7 @@ -rxtracto <- function(dataInfo, parameter = NULL, xcoord=NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xlen = 0., ylen = 0., zlen = 0., xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE) { +rxtracto <- function(dataInfo, parameter = NULL, xcoord=NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xlen = 0., ylen = 0., zlen = 0., xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'https://upwell.pfeg.noaa.gov/erddap', verbose = FALSE) { # Check Passed Info ------------------------------------------------------- rerddap::cache_setup(temp_dir = TRUE) @@ -127,8 +127,8 @@ rxtracto <- function(dataInfo, parameter = NULL, xcoord=NULL, ycoord = NULL, zco tcoordLim <- NULL if (!is.null(working_coords$tcoord1)) { isoTime <- dataCoordList$time - udtTime <- parsedate::parse_date(isoTime) - tcoord1 <- parsedate::parse_date(working_coords$tcoord1) + udtTime <- parsedate::parse_iso_8601(isoTime) + tcoord1 <- parsedate::parse_iso_8601(working_coords$tcoord1) tcoordLim <- c(min(tcoord1), max(tcoord1)) } @@ -199,12 +199,12 @@ latSouth <- working_coords$latSouth # the call will be the same as last time, so no need to repeat out_dataframe[i,] <- oldDataFrame } else { - griddapCmd <- makeCmd(urlbase, xName, yName, zName, tName, parameter, + griddapCmd <- makeCmd(dataInfo, urlbase, xName, yName, zName, tName, parameter, erddapCoords$erddapXcoord, erddapCoords$erddapYcoord, erddapCoords$erddapTcoord, erddapCoords$erddapZcoord, verbose ) - extract <- eval(parse(text = griddapCmd)) + extract <- do.call(rerddap::griddap, griddapCmd ) if (length(extract) == 0) { print(griddapCmd) diff --git a/R/rxtracto_3D.R b/R/rxtracto_3D.R index d5335bb..f194ad7 100644 --- a/R/rxtracto_3D.R +++ b/R/rxtracto_3D.R @@ -29,7 +29,7 @@ #' \item extract$time - the times of the extracts #' } #' @examples -#' urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +#' urlbase <- 'https://upwell.pfeg.noaa.gov/erddap' #' dataInfo <- rerddap::info('erdMBsstd8day') #' parameter <- 'sst' #' xcoord <- c(230, 235) @@ -74,7 +74,7 @@ #' yName = yName, zName = zName) #' -rxtracto_3D <- function(dataInfo, parameter = NULL, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'https://upwell.pfeg.noaa.gov/erddap', verbose=FALSE) { +rxtracto_3D <- function(dataInfo, parameter = NULL, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'https://upwell.pfeg.noaa.gov/erddap/', verbose=FALSE) { # Check Passed Info ------------------------------------------------------- @@ -120,8 +120,8 @@ if (!is.null(working_coords$tcoord1)) { isoTime <- dataCoordList$time udtTime <- parsedate::parse_date(isoTime) tcoord1 <- removeLast(isoTime, working_coords$tcoord1) - tcoord1 <- parsedate::parse_date(tcoord1) - tcoordLim <- c(min(tcoord1), max(tcoord1)) + tcoord1 <- parsedate::parse_iso_8601(tcoord1) + tcoordLim <- tcoord1 } dimargs <- list(xcoordLim, ycoordLim, zcoordLim, tcoordLim) @@ -145,7 +145,7 @@ erddapCoords <- erddapList$erddapCoords -griddapCmd <- makeCmd(urlbase, xName, yName, zName, tName, parameter, +griddapCmd <- makeCmd(dataInfo, urlbase, xName, yName, zName, tName, parameter, erddapCoords$erddapXcoord, erddapCoords$erddapYcoord, erddapCoords$erddapTcoord, erddapCoords$erddapZcoord, verbose ) @@ -154,7 +154,7 @@ griddapCmd <- makeCmd(urlbase, xName, yName, zName, tName, parameter, # Get the data ------------------------------------------------------------ -griddapExtract <- eval(parse(text = griddapCmd)) +griddapExtract <- do.call(rerddap::griddap, griddapCmd ) diff --git a/R/rxtractogon.R b/R/rxtractogon.R index baf2480..039fb57 100644 --- a/R/rxtractogon.R +++ b/R/rxtractogon.R @@ -53,7 +53,7 @@ -rxtractogon <- function(dataInfo, parameter, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE) { +rxtractogon <- function(dataInfo, parameter, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'https://upwell.pfeg.noaa.gov/erddap', verbose = FALSE) { rerddap::cache_setup(temp_dir = TRUE) diff --git a/R/utils.R b/R/utils.R index 16a37bd..f892ce7 100644 --- a/R/utils.R +++ b/R/utils.R @@ -28,6 +28,8 @@ strextract1 <- function(str, pattern) regmatches(str, regexpr(pattern, str)) strtrim1 <- function(str) gsub("^\\s+|\\s+$", "", str) +### End utilites from rerddap + findERDDAPcoord <- function(dataCoordList, isotime, udtime, xcoordLim, ycoordLim, tcoordLim, zcoordLim, xName, yName, tName, zName) { newxIndex <- rep(NA_integer_, 2) @@ -73,32 +75,44 @@ findERDDAPcoord <- function(dataCoordList, isotime, udtime, xcoordLim, ycoordLim } -makeCmd <- function(urlbase, xName, yName, zName, tName, parameter, +makeCmd <- function(dataInfo, urlbase, xName, yName, zName, tName, parameter, erddapXcoord, erddapYcoord, erddapTcoord, erddapZcoord, verbose ) { - myCallOpts <- "" - if (!(urlbase == "http://upwell.pfeg.noaa.gov/erddap")) { - myCallOpts <- paste0(", url='", urlbase,"/'") + + myCallOpts <- list(dataInfo) + myCallOptsNames <- list('x') + if (!(urlbase == "https://upwell.pfeg.noaa.gov/erddap/")) { + myCallOpts$url <- urlbase + myCallOptsNames <- c(myCallOptsNames, 'url') } if (verbose) { - myCallOpts <- paste0(myCallOpts,", callopts = httr::verbose()") + myCallOpts$callopts <- httr::verbose() + myCallOptsNames <- c(myCallOptsNames, 'callopts') } - griddapCmd <- 'rerddap::griddap(dataInfo,' if (!is.na(erddapXcoord[1])) { - griddapCmd <- paste0(griddapCmd, xName,'=c(',erddapXcoord[1],',', erddapXcoord[2],'),') - } - if (!is.na(erddapYcoord[1])) { - griddapCmd <- paste0(griddapCmd, yName,'=c(',erddapYcoord[1],',', erddapYcoord[2],'),') + myCallOpts$xName <- erddapXcoord + myCallOptsNames <- c(myCallOptsNames, xName) } + if (!is.na(erddapYcoord[1])) { + myCallOpts$yName <- erddapYcoord + myCallOptsNames <- c(myCallOptsNames, yName) + } if (!is.na(erddapZcoord[1])) { - griddapCmd <- paste0(griddapCmd, zName,'=c(',erddapZcoord[1],',', erddapZcoord[2],'),') + myCallOpts$zName <- erddapZcoord + myCallOptsNames <- c(myCallOptsNames, zName) } if (!is.na(erddapTcoord[1])) { - griddapCmd <- paste0(griddapCmd, tName,'=c("',erddapTcoord[1],'","', erddapTcoord[2],'"),') + myCallOpts$tName <- erddapTcoord + myCallOptsNames <- c(myCallOptsNames, tName) } - griddapCmd <- paste0(griddapCmd,'fields="', parameter,'",read = FALSE', myCallOpts,')') + myCallOpts$fields = parameter + myCallOptsNames <- c(myCallOptsNames, 'fields') + myCallOpts$read <- FALSE + myCallOptsNames <- c(myCallOptsNames, 'read') + + names(myCallOpts) <- myCallOptsNames - return(griddapCmd) + return(myCallOpts) } removeLast <- function(isotime, tcoord1) { diff --git a/README.md b/README.md index 6c5960c..18535a8 100644 --- a/README.md +++ b/README.md @@ -1,38 +1,53 @@ -# rerddapXtracto +# rerddapXtracto (Version 0.3.0) rerddapXtracto - R package for accessing environmental data using rerddap (**For Testing Purposes Only**) -** Version. Still for testing purposes only ** -** For real work - use at own risk ** -** Suggest doing small tests on any extract ** -** For rxtracto(), do say the first 3 points with "verbose = TRUE" ** -** Make certain the URL calls amke sense ** -** In final output, make certain data limits ** -** and limits requested from dataset make sense ** +- ** Version. Still for testing purposes only ** +- ** For real work - use at own risk ** +- ** Suggest doing small tests on any extract ** +- ** For rxtracto(), do say the first 3 points with "verbose = TRUE" ** +- ** Make certain the URL calls amke sense ** +- ** In final output, make certain data limits ** +- ** and limits requested from dataset make sense ** -`rerddapXtracto` is an R package developed to subset and extract satellite and other oceanographic related data from a remote ERDDAP server. The program can extract data for a moving point in time along a user-supplied set of longitude, latitude and time points; in a 3D bounding box; or within a polygon (through time). +`rerddapXtracto` is an R package developed to subset and extract satellite and other oceanographic related data from a remote ERDDAP server. The program can extract data for a moving point in time along a user-supplied set of longitude, latitude, time and depth (new in thos version) points; in a 3D bounding box; or within a polygon (through time). -These functions differ from those in the xtractomatic package in that they use the `rerddap` package to access gridded data on any ERDDAP server, but they require the user to provide initial information about the data to be extracted. +New in this version is that a track can now move in (x, y, z, t) space if appropriae for the dataset being accessed. And two plotting functions have been added, `plotTrack()` and `plotBox()` that make use of the `plotdap` package. See the new [rerdapXtracto vignette](https://rmendels.github.io/UsingrerddapXtracto.nb.html). A lot of the code has been reworked, in particular the handling of time, and in the formation of the requests to `rerddap`. -There are three main data extraction functions in the `xtractomatic` package: -- `rxtracto <- function(dataInfo, parameter = NULL, xcoord=NULL, ycoord=NULL, zcoord = NULL, tcoord=NULL, xlen=0., ylen=0., xName='longitude', yName='latitude', zName='altitude', tName='time', urlbase='http://upwell.pfeg.noaa.gov/erddap')` +There are three main data extraction functions in the `rerddapXtracto` package: -- `rxtracto_3D <- function(dataInfo, parameter = NULL, xcoord=NULL, ycoord=NULL, zcoord = NULL, tcoord=NULL, xName='longitude', yName='latitude', zName='altitude', tName='time', urlbase='http://upwell.pfeg.noaa.gov/erddap')` +- `rxtracto <- function(dataInfo, parameter = NULL, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xlen = 0., ylen = 0., zlen = 0., xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE)` -- `rxtractogon <- function(dataInfo, parameter, xcoord=NULL, ycoord=NULL, zcoord = NULL, tcoord=NULL, xName='longitude', yName='latitude', zName='altitude', tName='time', urlbase='http://upwell.pfeg.noaa.gov/erddap')` +- `rxtracto_3D <- function(dataInfo, parameter = NULL, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE)` +- `rxtractogon <- function(dataInfo, parameter, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE)` +and two functions for producing maps: -`rerddapXtracto` uses the `rerddap`, `ncdf4` and `sp` packages , and these packages (and the packages imported by these packages) must be installed first or `rerddapXtracto` will fail to install. +- `plotTrack <- function(xcoord, ycoord, resp, plotColor = 'viridis', name = NA, myFunc = NA, shape = 20, size = .5)` + +- `plotBox <- function(resp, plotColor = 'viridis', time = NA, animate = FALSE, name = NA, myFunc = NA, maxpixels = 10000)` + + + +`rerddapXtracto` uses the `rerddap`, `ncdf4` , `parsedata`, `plotdap` and `sp` packages , and these packages (and the packages imported by these packages) must be installed first or `rerddapXtracto` will fail to install. ```{r install,eval=FALSE} install.packages("rerddap", dependencies = TRUE) install.packages("ncdf4") +install.packages("parsedata") install.packages("sp") ``` -The `xtractomatic` package at the moment can be installed from Github using the devtools package: +To install the plotdap package from Github: + +```{r plotdap, eval = FALSE} +install.packages("devtools") +devtools::install_github('ropensci/plotdap') +``` + +The `rerddapXtracto` package at the moment can be installed from Github using the devtools package: ```{r install,eval=FALSE} install.packages("devtools") diff --git a/man/plotBox.Rd b/man/plotBox.Rd new file mode 100644 index 0000000..9fa9308 --- /dev/null +++ b/man/plotBox.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotBox.R +\name{plotBox} +\alias{plotBox} +\title{plot result of xtracto_3D or rxtracto_3D} +\usage{ +plotBox(resp, plotColor = "viridis", time = NA, animate = FALSE, + name = NA, myFunc = NA, maxpixels = 10000) +} +\arguments{ +\item{resp}{data frame returned from rxtracto() or xtracto()} + +\item{plotColor}{the color to use in plot from rerddap} + +\item{time}{a function to map multi-time to one, or else identity for animation} + +\item{animate}{if multiple times, if TRUE will animate the maps} + +\item{myFunc}{function of one argument to transform the data} + +\item{maxpixels}{maximum numbe rof pixels to use in making the map - controls resolution} +} +\value{ +a plotdap plot +} +\description{ +\code{plotBox} is a function to plot the results from +rxtracto() and xtracto() +} +\examples{ +tagData <- Marlintag38606 +xpos <- tagData$lon +ypos <- tagData$lat +tpos <- tagData$date +zpos <- rep(0., length(xpos)) +urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +swchlInfo <- rerddap::info('erdSWchla8day') +swchl <- rxtracto(swchlInfo, parameter = 'chlorophyll', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, xlen = .2, ylen = .2) +plotBox(xpos, ypos, swchl, plotColor = 'chlorophyll') +} diff --git a/man/plotTrack.Rd b/man/plotTrack.Rd new file mode 100644 index 0000000..5fc1099 --- /dev/null +++ b/man/plotTrack.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotTrack.R +\name{plotTrack} +\alias{plotTrack} +\title{plot result of xtracto or rxtracto} +\usage{ +plotTrack(xcoord, ycoord, resp, plotColor = "viridis", name = NA, + myFunc = NA, shape = 20, size = 0.5) +} +\arguments{ +\item{xcoord}{passed to rxtracto() or xtracto()} + +\item{ycoord}{passed to rxtracto() or xtracto()} + +\item{resp}{data frame returned from rxtracto() or xtracto()} + +\item{plotColor}{the color to use in plot from rerddap} + +\item{myFunc}{function of one argument to transform the data} +} +\value{ +a plotdap plot +} +\description{ +\code{plotTrack} is a function to plot the results from +rxtracto() and xtracto() +} +\examples{ +tagData <- Marlintag38606 +xpos <- tagData$lon +ypos <- tagData$lat +tpos <- tagData$date +zpos <- rep(0., length(xpos)) +urlbase <- 'http://upwell.pfeg.noaa.gov/erddap' +swchlInfo <- rerddap::info('erdSWchla8day') +swchl <- rxtracto(swchlInfo, parameter = 'chlorophyll', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, xlen = .2, ylen = .2) +plotTrack(xpos, ypos, swchl, plotColor = 'chlorophyll') +} diff --git a/vignettes/UsingrerddapXtracto.Rmd b/vignettes/UsingrerddapXtracto.Rmd index a638404..75e3ed4 100644 --- a/vignettes/UsingrerddapXtracto.Rmd +++ b/vignettes/UsingrerddapXtracto.Rmd @@ -19,7 +19,9 @@ library(rerddapXtracto) `rerddapXtracto` is an R package developed to subset and extract satellite and other oceanographic related data from any ERDDAP server using the R package [rerddap](https://cran.r-project.org/web/packages/rerddap/index.html) developed by Scott Chamberlain and the wonderful people at [rOpenSci](https://ropensci.org). ERDDAP is a simple to use yet powerful web data service developed by Bob Simons. `rerddapXtracto` extends the `rerddap` package by being able to extract data for a moving point in time along a user-supplied set of longitude, latitude and time points; and also extracting data within a polygon (through time). `rerddapXtracto` extends the functionality of the R package [xtractomatic](https://cran.r-project.org/web/packages/xtractomatic/index.html) by being able to work with (hopefully) most gridded datasets available from any ERDDAP server. The disadavantage compared to `xtractomatic` is that the user has to do more work to obtain information about the dataset to be accessed, while in `xtractomatic` that information is built in. -This version has one main change - that is in both `rxtracto()` and `rxtracto_3D()` the zcoord is not limited to be at a set location. That means for `rxtracto_3D()` that if the zCoord needs to be given for any reason, then it must be of length two, and for `rxtracto()` if the zCoord needs to be given for any reason, it must be of the same length as the other coordinates, and can also have a "zlen"", like "xlen" and "ylen", that defines a bounding box within which to make the extract. The advantage of this is it allows `rxtracto()` to mkae extracts moving in (x, y, z, t) space. +This version has several major changes. In both `rxtracto()` and `rxtracto_3D()` the zcoord is not limited to be at a set location. That means for `rxtracto_3D()` that if the zCoord needs to be given for any reason, then it must be of length two, and for `rxtracto()` if the zCoord needs to be given for any reason, it must be of the same length as the other coordinates, and can also have a "zlen"", like "xlen" and "ylen", that defines a bounding box within which to make the extract. The advantage of this is it allows `rxtracto()` to make extracts moving in (x, y, z, t) space. + +Second there are now two functions, `plotTrack()` for tracks and `plotBox()` for grids, that produce quick maps of the output using the R package `plotdap`. All of the examples have been changed to use these functions. Bounding box extracts, as in `rxtracto_3D()`, can be done just using `rerddap`, but the `rerddap` function `griddap()` returns a "melted" version of the data, while `rxtracto_3D()` retains the grid structure. Also, `rxtractogon()` uses `rxtracto_3D()` and therefore `plotBox()` can be used to plot the the results from `rxtractogon()`. These functions should also work with the output from the R package `xtractomatic`. ### The Main xtractomatic functions @@ -31,7 +33,13 @@ There are three main data extraction functions in the `rerddapXtracto` package: - `rxtractogon <- function(dataInfo, parameter, xcoord = NULL, ycoord = NULL, zcoord = NULL, tcoord = NULL, xName = 'longitude', yName = 'latitude', zName = 'altitude', tName = 'time', urlbase = 'http://upwell.pfeg.noaa.gov/erddap', verbose = FALSE)` -The functions are similar to but not identical to the functions in `xtractomatic`. The main differences are having to obtain information about the dataset first using the function `rerddap::info()`, and possibly having to give the names of the coordinate variables, as these can't be assumed (for example the zcoord could be in sigma coordinates). More specifically: +and two functions for producing maps: + +- `plotTrack <- function(xcoord, ycoord, resp, plotColor = 'viridis', name = NA, myFunc = NA, shape = 20, size = .5)` + +- `plotBox <- function(resp, plotColor = 'viridis', time = NA, animate = FALSE, name = NA, myFunc = NA, maxpixels = 10000)` + +The data extraction functions are similar to but not identical to the functions in `xtractomatic`. The main differences are having to obtain information about the dataset first using the function `rerddap::info()`, and possibly having to give the names of the coordinate variables, as these can't be assumed (for example the zcoord could be in sigma coordinates). More specifically: - dataInfo: the return from an `rerddap::info()` call to a dataset on an ERDDAP server - parameter: character string containing the name of the parameter to extract @@ -52,10 +60,10 @@ With all due respect to the [Chambers Brothers](https://www.youtube.com/watch?v= `rerddapXtracto` uses the R packages `rerddap`, `ncdf4`, `sp` and `parsedate`, and these packages (and the packages imported by these packages) must be installed first or `rerddapXtracto` will fail to install. ```{r install, eval = FALSE} -install.packages("rerddap", dependencies = TRUE) install.packages("ncdf4", dependencies = TRUE) -install.packages("sp", dependencies = TRUE) install.packages("parsedate", dependencies = TRUE) +install.packages("rerddap", dependencies = TRUE) +install.packages("sp", dependencies = TRUE) ``` The `rerddapXtracto` package is not available through CRAN at the moment, but it is available from [Github](https://github.com/rmendels/rerddapXtracto) @@ -73,31 +81,53 @@ Once installed, to use `rerddapXtracto`: library("rerddapXtracto") ``` +To install the plotdap package from Github: + +```{r plotdap, eval = FALSE} +install.packages("devtools") +devtools::install_github('ropensci/plotdap') +``` + +Note that plotdap depends on a number of packages that must be installed. These include the packages ggplot2, raster and sf. To use the animation features, `gganimate` must be installed, see [gganimate](https://github.com/dgrtwo/gganimate). `gganimate` requires a version of ImageMagick to be installed on your computer. + If the other R libraries (`rerddap`, `ncdf4`,`sp`, and `parsedate`) have been installed they will be found and do not need to be explicitly loaded. ### Using the R code examples -Besides the `rerddapXtracto` package, the examples below depend on the R packages `akima`, `dplyr`, `ggplot2`, `lubridate`, `mapdata`, and `xts`. These can be loaded beforehand (assuming they have been installed): +Once installed, to use rerddapXtracto: -```{r loadpacks, eval = FALSE} -library("akima") -library("dplyr") +```{r, eval = FALSE} +library("rerddapXtracto") +``` + +and to use the plotting functions: + +```{r, eval = FALSE} library("ggplot2") -library("lubridate") -library("mapdata") -library("xts") +library("plotdap") +library("sf") ``` -In order that the code snippets be more stand-alone, the needed libraries are always `required()` in the examples. -It should be emphasized that these other packages are used to manipulate and plot the data in R, other packages could be used as well. The use of `rerddapXtracto` does not depend on these other packages. +## Getting Started -There are also several R functions defined within the document that are used in other code examples. These include `mapFrame()`, and `plotFrame()`. +The plotting functions are new, and there are some fine points that neeed to be understood if they are to be used properly, in particular `plotBox()`. Both plotTrack() and `plotBox()` rearrange the output so that the plotdap functions add_tabledap() and `add_griddap()` think that the output is from `rerddap`, and then make the appropriate `plotdap` call. When the data that is passed to `add_griddap()` has mutiple time periods, there are two options. The first option is to set the parameter “time” to a function that reduces the data to one dimension in the time coordinate (such as the mean), or else to set “time” equal to “identity” and set “animate” to be “TRUE” which will produce a time animation of the results. The function `plotBox()` works the same way, except that the default function is `mean(na.rm = TRUE)`. The following link to examples that show how to use different features of the plotting functions: + +- [Setting the color palette](#colorPalette) shows how to use the “plotColor” option. The “plotColor” parameter can be the name of any of the colors included in the rerddap color pallete. These colors are based on the cmocean colormaps designed by Kristen Thyng (see http://matplotlib.org/cmocean/ and https://github.com/matplotlib/cmocean), which were initally developed for Python, but a version of the colormaps is used in the oce package by Dan Kelley and Clark Richards and that is also what is used in rerddap. + +- [Plot one time period](#plot1) example shows how to manipulate an existing output from `rxtracto_3D()` or `rextractogon()` to plot just one time period. + +- [Transform the data](#transform) example shows how to use the “myFunc” option to transform the data before plotting. The function has to be a function of a single argument. This example also shows how to use the “name” option to chnage the name displayed on the color bar. In this example, we want depth to go downwards in the colorbar, and the name given changed from “altitude”, which is the name on ERDDAP, to the name “Depth”. + +- [Name](#name) example shows how to change the name on the colorbar. + +- [Modify the graph](#modify) shows how to use the `plotdap` functon `add_ggplot()` to modify a graph once it has been generated. + +- [Animate](#animate) shows how to animate a grid with multiple time periods. -## Getting Started The first several examples reproduce some of the examples in the `xtractomatic` vignette, hopefully to make clear how the functions in the two packages relate. One change from the `xtractomatic` vignette is that the plots use the `cmocean` colormaps designed by Kristen Thyng (see http://matplotlib.org/cmocean/ and https://github.com/matplotlib/cmocean). These colormaps were initally developed for Python, but a version of the colormaps is used in the `oce` package by Dan Kelley and Clark Richards and that is what is used here. @@ -123,7 +153,7 @@ In this section we extract data along a trackline found in the chlorophyll around that point. Positions where there was a tag location but no chlorophyll values are also shown. +We plot the track line with the locations colored according to the mean of the satellite `chlorophyll` around that point. Positions where there was a tag location but no `chlorophyll` values are also shown. This example shows the use of the “plotColor” parameter to use the “chlorophyll” color palette. -```{r meantrackPlot, fig.align = 'center', fig.width = 6, fig.height = 3} + +```{r meantrackPlot, fig.align = 'center', fig.width = 6, fig.height = 5, message = FALSE} require("ggplot2") -require("mapdata") -# First combine the two dataframes (the input and the output) into one, -# so it will be easy to take into account the locations that didn’t -# retrieve a value. - -alldata <- cbind(tagData, swchl) - -# adjust the longitudes to be (-180, 180) -alldata$lon <- alldata$lon - 360 -# Create a variable that shows if chla is missing -alldata$missing <- is.na(alldata$'mean chlorophyll') * 1 -alldata$mean <- as.array(alldata$'mean chlorophyll') -# set limits of the map -ylim <- c(15, 30) -xlim <- c(-160, -105) -myColor <- colors$chlorophyll -# get outline data for map -w <- map_data("worldHires", ylim = ylim, xlim = xlim) -# plot using ggplot -z <- ggplot(alldata,aes(x = lon, y = lat)) + - geom_point(aes(colour = mean, shape = factor(missing)), size = 2.) + - scale_shape_manual(values = c(19, 1)) -z + geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + - scale_colour_gradientn(colours = myColor,limits = c(0., 0.32), "Chla") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + ggtitle("Mean chla values at marlin tag locations") +require("plotdap") +require("sf") + +myPlot <- plotTrack(xpos, ypos, swchl, plotColor = 'chlorophyll') +myPlot ``` ### Topography data -The second example from the `xtractomatic` vignette is accessing topographic data along the Marlin track. +The second example from the `xtractomatic` vignette is accessing topographic data along the Marlin track. This example alos shows how to pass a function to `plotTrack` to transform the data before plotting, and to change the name shown on the colorbar. -```{r topotag, fig.align = 'center', fig.width = 6, fig.height = 4, warning = FALSE} + +```{r topotag, fig.align = 'center', fig.width = 6, fig.height = 5, warning = FALSE} require("ggplot2") -require("mapdata") +require("plotdap") require("rerddap") require("rerddapXtracto") +require("sf") ylim <- c(15, 30) xlim <- c(-160, -105) topoInfo <- rerddap::info('etopo360') topo <- rxtracto(topoInfo, parameter = 'altitude', xcoord = xpos, ycoord = ypos, xlen = .1, ylen = .1) -alldata <- cbind(tagData, topo) -alldata$lon <- alldata$lon - 360 -alldata$mean <- as.array(alldata$'mean altitude') -# get outline data for map -w <- map_data("worldHires", ylim = ylim, xlim = xlim) -z <- ggplot(alldata,aes(x = lon, y = lat)) + - geom_point(aes(colour = mean), size = 2.) + - scale_shape_manual(values = c(19, 1)) -z + geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + - scale_colour_gradient( name = "Depth") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + ggtitle("Bathymetry at marlin tag locations") - +myFunc = function(x) -x +topoPlot <- plotTrack(xpos, ypos, topo, plotColor = 'density', name = 'Depth', myFunc = myFunc) +topoPlot ``` Again, note the differences from the `xtractomatic` example: @@ -217,6 +219,23 @@ Again, note the differences from the `xtractomatic` example: - the user has to know the parameter name desired for that datasetID, in this case 'altitude' - only the xcoord and ycoord values need to be given, as 'etopo360' is a 2-D dataset, and this is handled more gracefully then in `xtractomatic`. +The following is an artificial example showing a track moving in (x, y, z, t) space. Since the times of the model output change, the actual times are retrieved, and the last three times used in the example. + +```{r extract3D} +require("rerddap") +urlBase <- "http://erddap.marine.ie/erddap/" +parameter <- "Sea_water_temperature" +dataInfo <- rerddap::info("IMI_CONN_3D", url = urlBase) +#get the actual last 3 times, and extract from data frame +dataInfo1 <- read.csv("https://erddap.marine.ie/erddap/griddap/IMI_CONN_3D.csv0?time[last-2:1:last]",stringsAsFactors = FALSE, header = FALSE, row.names = NULL) +sstTimes <- dataInfo1[[1]] +sstLats <- c(53.505758092414446, 53.509303546859805, 53.51284900130517) +sstLons <- c(-10.25975390624996, -10.247847656249961, -10.23594140624996) +sstDepths <- c(2, 6, 10) +sstTrack <- rxtracto(dataInfo, parameter = parameter, xcoord = sstLons, ycoord = sstLats, tcoord = sstTimes, zcoord = sstDepths, xlen = .05, ylen = .05, zlen = 0., zName = 'altitude', urlbase = urlBase) +str(sstTrack) +``` + ## Using `rxtracto_3D` @@ -224,9 +243,9 @@ The function `rxtracto_3D()` adds no new capabilites to `rerddap`, but it does r ### Obtaining VIIRS chlorophyll data +We examine VIIRS chlorophyll for the “latest” data as of when the vignette was generated: -```{r VIIRSchla} -require("lubridate") +```{r VIIRSchla, warning = FALSE, message = FALSE} require("rerddap") require("rerddapXtracto") @@ -236,7 +255,6 @@ tpos <- c("last", "last") tpos <- c("2017-04-15", "2017-04-15") VIIRSInfo <- rerddap::info('erdVH3chlamday') VIIRS <- rxtracto_3D(VIIRSInfo, parameter = 'chla', xcoord = xpos, ycoord = ypos, tcoord = tpos) -VIIRS$time <- lubridate::as_date(VIIRS$time) ``` `rxtracto_3d()` returns a list of the form: @@ -250,52 +268,18 @@ VIIRS$time <- lubridate::as_date(VIIRS$time) The coordinate names of the structure are based on the names given in the `rxtracto_3d()` call, so may differ between datasets. A similar call to `rerddap::griddap()` will either return the data "pre-melted" (that is long-form) or only get the netcdf file and have the user read in the data. There are tradeoffs to having the data "pre-melted", for consistency we maintain a structure similar to that in `xtracto_3D()`, plus this works even if the xcoord and ycoord are not longitude and latitude, where the "pre-melting" in `rerddap::griddap()` fails. -We can map the results using a "helper" function `mapFrame()` that "melts" the data into "longform": -```{r mapFrame} -mapFrame <- function(longitude, latitude, chla) { - dims <- dim(chla) - chla <- array(chla, dims[1] * dims[2]) - chlaFrame <- expand.grid(x = longitude, y = latitude) - chlaFrame$chla <- chla - return(chlaFrame) -} -``` - -and also define a helper function `plotFrame()` to plot the data: - -```{r plotFrame} -plotFrame <- function(chlaFrame, xlim, ylim, title, logplot = TRUE) { - require("ggplot2") - require("mapdata") - w <- map_data("worldHires", ylim = ylim, xlim = xlim) - myplot <- ggplot(data = chlaFrame, aes(x = x, y = y, fill = chla)) + - geom_raster(interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) - if (logplot) { - my.col <- colors$chlorophyll - myplot <- myplot + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(-2, 4)) + - ggtitle(title) - }else{ - myplot < -myplot + scale_fill_gradientn(colours = my.col, na.value = NA) + - ggtitle(title) - } - return(myplot) -} -``` - -We examine chlorophyll the "latest" data as of when the vignette was generated: +We can map the data using `plotBox()`: + ```{r VIIRSLogPlot, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE} -xlim <- c(-125, -120) -ylim <- c(36, 39) -chlalogFrame <- mapFrame(VIIRS$longitude, VIIRS$latitude, - log(VIIRS$chla[, , 1])) -chlalogPlot <- plotFrame(chlalogFrame, xlim, ylim, "log(VIIRS chla)") +require("ggplot2") +require("plotdap") +require("sf") +myFunc <- function(x) log(x) +chlalogPlot <- plotBox(VIIRS, plotColor = 'chlorophyll', myFunc = myFunc) chlalogPlot ``` @@ -320,31 +304,40 @@ sanctchl <- rxtractogon(dataInfo, parameter = parameter, xcoord = xpos, ycoord = str(sanctchl) ``` -The extract (see `str(sanctchl)`) contains two time periods of chlorophyll masked for data only in the sanctuary boundaries. A plot of the the second time period: +The extract (see `str(sanctchl)`) contains two time periods of chlorophyll masked for data only in the sanctuary boundaries. This example shows how to pull out only a single time period to be used in `plotBox()`. -```{r mbnmsChlaPlot, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE} + +```{r mbnmsChlaPlot, fig.width = 6, fig.height = 3, fig.align = 'center', warning = FALSE} require("ggplot2") -require("mapdata") -xlim <- c(-123.5, -121.) -ylim <- c(35, 38) -mbnmsFrame <- mapFrame(sanctchl$longitude,sanctchl$latitude, log(sanctchl$chla[, , 2])) -my.col <- colors$chlorophyll -w <- map_data("worldHires", ylim = ylim, xlim = xlim) -myplot <- ggplot() + geom_path(data = mbnms,aes(x = Longitude, y = Latitude), colour = "black") -myplot <- myplot + - geom_raster(data = mbnmsFrame, aes(x = x, y = y, fill = chla), interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(-1, 3)) + - ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + - ggtitle(paste("log(Chla) in MBNMS", sanctchl$time[2])) -myplot +require("plotdap") +require("sf") +myFunc <- function(x) log(x) +sanctchl1 <- sanctchl +sanctchl1$chla <- sanctchl1$chla[, , 2] +sanctchl1$time <- sanctchl1$time[2] +sanctchlPlot <- plotBox(sanctchl1, plotColor = 'chlorophyll', myFunc = myFunc) +sanctchlPlot +``` + +This extract can be used to show the ability to animate the output through time: + + +```{r animate, eval = FALSE} +require("gganimate") +#> Loading required package: gganimate +require("ggplot2") +require("plotdap") +require("sf") +myFunc <- function(x) log(x) +sanctchlPlot <- plotBox(sanctchl, plotColor = 'chlorophyll', myFunc = myFunc, time = identity, animate = TRUE) ``` +![Sanctuary Animation](ani.gif) + The MBNMS is famous for containing the Monterey Canyon, which reaches depths of up to 3,600 m (11,800 ft) below surface level at its deepest. `rxtractogon()` can extract the bathymetry data for the MBNMS from the ETOPO dataset: -```{r mbnmsBathy} +```{r mbnmsBathy, warning = FALSE} require("rerddap") dataInfo <- rerddap::info('etopo180') xpos <- mbnms$Longitude @@ -355,32 +348,24 @@ str(bathy) Mapping the data to show the canyon: -```{r mbnmsBathyPlot, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE} +```{r mbnmsBathyPlot, fig.width = 5, fig.height = 5, fig.align = 'center', warning = FALSE, message = FALSE} require("ggplot2") require("mapdata") -xlim <- c(-123.5, -121.) -ylim <- c(35, 38) -mbnmsFrame <- mapFrame(bathy$longitude, bathy$latitude, bathy$depth[, , 1]) -w <- map_data("worldHires", ylim = ylim, xlim = xlim) -myplot <- ggplot() + geom_path(data = mbnms, aes(x = Longitude, y = Latitude), colour = "black") -myplot <- myplot + geom_raster(data = mbnmsFrame, aes(x = x, y = y, fill = chla),interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + scale_fill_gradient(na.value = NA, name = "Depth") + - ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + ggtitle("MBNMS Bathymetry") -myplot +myFunc = function(x) -x +bathyPlot <- plotBox(bathy, plotColor = 'density', myFunc = myFunc, name = 'Depth') +bathyPlot ``` ## Temperature at 70m in the north Pacific from the SODA model output -This is an example of an extract from a 4-D dataset (results from the "Simple Ocean Data Assimilation (SODA)" model), and illustrate the case where the z-coordinate does not have the default name "altitude". Water temperature at 70m depth is extracted for the North Pacific Ocean. +This is an example of an extract from a 4-D dataset (results from the "Simple Ocean Data Assimilation (SODA)" model), and illustrate the case where the z-coordinate does not have the default name "altitude". Water temperature at 70m depth is extracted for the North Pacific Ocean east of the dateline. ```{r soda70} require("rerddap") dataInfo <- rerddap::info('erdSoda331oceanmday') -xpos <- c(135.25, 240.25) +xpos <- c(185.25, 240.25) ypos <- c(20.25, 60.25) zpos <- c(76.80285, 76.80285) tpos <- c('2010-12-15', '2010-12-15') @@ -388,28 +373,15 @@ soda70 <- rxtracto_3D(dataInfo, parameter = 'temp', xcoord = xpos, ycoord = ypos str(soda70) ``` -Since the data cross the dateline, it is necessary to use the new "world2Hires" continental outlines in the package "mapdata" which is Pacific Ocean centered. Unfortunatley there is a small problem where the outlines from certain countries wrap and mistakenly appear in plots, and those countries must be removed, see code below. ```{r soda70Plot, fig.width = 6, fig.height = 3, fig.align = 'center', warning = FALSE} require("ggplot2") -require("mapdata") -xlim <- c(135, 240) -ylim <- c(20, 60) -soda70Frame <- mapFrame(soda70$longitude,soda70$latitude, soda70$temp[, , 1, 1]) -my.col <- colors$temperature -## Must do a kludge to remove countries that wrap and mess up the plot -w1 <- map("world2Hires", xlim = c(135, 240), ylim = c(20, 60), fill = TRUE, plot = FALSE) -remove <- c("UK:Great Britain", "France", "Spain", "Algeria", "Mali", "Burkina Faso", "Ghana", "Togo") -w <- map_data("world2Hires", regions = w1$names[!(w1$names %in% remove)], ylim = ylim, xlim = xlim) -myplot <- ggplot() + - geom_raster(data = soda70Frame, aes(x = x, y = y, fill = chla), interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(-3,30), name = "temperature") + - ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + - ggtitle(paste("temperature at 70 meters depth from SODA for", soda70$time[1])) -myplot +require("plotdata") +require("sf") +sodaPlot <- plotBox(soda70, plotColor = 'temperature', name = 'temp_at_70m', maxpixels = 30000) +sodaPlot + ``` @@ -433,22 +405,21 @@ str(NAtlSSS) ```{r NAtlSSSplot, fig.width = 6, fig.height = 3, fig.align = 'center', warning = FALSE} require("ggplot2") -require("mapdata") -xlim <- c(-17.99375, -1.00625) -ylim <- c(48.00625, 57.50625) -NAtlSSSFrame <- mapFrame(NAtlSSS$longitude, NAtlSSS$latitude, NAtlSSS$sea_surface_salinity[, , 1]) -my.col <- colors$salinity -w <- map_data("worldHires", ylim = ylim, xlim = xlim) -myplot <- ggplot() + - geom_raster(data = NAtlSSSFrame, aes(x = x, y = y, fill = chla), interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(34, 36), name = "salinity") + - ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + - ggtitle(paste("salinity", NAtlSSS$time[1])) -myplot +require("plotdap") +NAtlSSSPlot <- plotBox(NAtlSSS, plotColor = 'salinity', name = "salinity", maxpixels = 30000) +NAtlSSSPlot ``` +A lot of the details in the ocean are hidden in the plot above, because there are some low salinity values right close to shore. The plot can be modified using the plotdap function add_ggplot() so that only values between (32, 36) are plotted, and to change the colorbar to reflect this. + + +```{r NAtlSSSplot1, fig.width = 6, fig.height = 3, fig.align = 'center', warning = FALSE} +require("ggplot2") +require("plotdap") +add_ggplot(NAtlSSSPlot, scale_colour_gradientn(colours = colors$salinity, na.value = NA, limits = c(32, 36)), scale_fill_gradientn(colours = colors$salinity, na.value = NA, limits = c(32, 36))) +``` + + ### IFREMER The French agency IFREMER also has an ERDDAP server. We obtain salinity data at 75 meters from the Global Ocean, Coriolis Observation Re-Analysis CORA4.1 model off the west coast of the United States. @@ -457,7 +428,7 @@ The French agency IFREMER also has an ERDDAP serv require("rerddap") urlBase <- "http://www.ifremer.fr/erddap/" parameter <- "PSAL" -ifrTimes <- c("2013-05-15", "2013-05-15") +ifrTimes <- c("2013-09-15", "2013-09-15") ifrLats <- c(30., 50.) ifrLons <- c(-140., -110.) ifrDepth <- c(75., 75.) @@ -466,35 +437,13 @@ ifrPSAL <- rxtracto_3D(dataInfo, parameter = parameter, xcoord = ifrLons, ycoord str(ifrPSAL) ``` -The `ggplot2` function `geom_raster()` is not designed for unevenly spaced coordinates, as are the latitudes from this model. The function `interp()` from the package `akima` is used to interpolate the data which are then plotted. - +Plotting th results using `plotBox()`: ```{r ifrPSALplot, fig.width = 6, fig.height = 3, fig.align='center', warning = FALSE} -## ggplot2 has trouble with unequal y's - require("akima") - require("dplyr") - require("ggplot2") - require("mapdata") - xlim <- c(-140, -110) - ylim <- c(30, 51) -## ggplot2 has trouble with unequal y's - my.col <- colors$salinity - tempData1 <- ifrPSAL$PSAL[, , 1, 1] - tempData <- array(tempData1 , 61 * 54) - tempFrame <- expand.grid(x = ifrPSAL$longitude, y = ifrPSAL$latitude) - tempFrame$temp <- tempData - tempFrame1 <- dplyr::filter(tempFrame, !is.nan(temp)) - myinterp <- akima::interp(tempFrame1$x, tempFrame1$y, tempFrame1$temp, xo = seq(min(tempFrame1$x), max(tempFrame1$x), length = 61), yo = seq(min(tempFrame1$y), max(tempFrame1$y), length = 54)) - myinterp1 <- expand.grid(x = myinterp$x, y = myinterp$y) - myinterp1$temp <- array(myinterp$z, 61 * 54) - w <- map_data("worldHires", ylim = ylim, xlim = xlim) - myplot <- ggplot() + - geom_raster(data = myinterp1, aes(x = x, y = y, fill = temp), interpolate = FALSE) + - geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + - theme_bw() + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(32, 35), name = "salinity") + - ylab("latitude") + xlab("longitude") + - coord_fixed(1.3, xlim = xlim, ylim = ylim) + ggtitle(paste("salinity at 75 meters",ifrPSAL$time[1] )) - myplot +require("ggplot2") +require("plotdap") +ifrPSALPlot <- plotBox(ifrPSAL, plotColor = 'salinity', name = "salinity", maxpixels = 30000) +ifrPSALPlot ``` diff --git a/vignettes/ani.gif b/vignettes/ani.gif new file mode 100644 index 0000000..c9514cc Binary files /dev/null and b/vignettes/ani.gif differ