Getting Started
There are some fine points that need to be understood to properly use
the plotting functions, in particular plotBBox(). Both
plotTrack() and plotBBox() rearrange the
output so that the functions plotdap::add_tabledap() and
plotdap::add_griddap() think that the output is from
rerddap, and then make the appropriate plotdap
call. When the data that are passed to add_griddap() has
multiple 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. If an animation is requested
and the option “cumulative” is set to be “TRUE”, then the animation will
be cumulative. This is a nice feature for displaying tracks through
time. The function plotBBox() 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 shows how to use the “plotColor” option. The “plotColor” parameter can be the name of any of the colors included in the
cmoceanpackage. Thecmoceancolormaps were designed by Kristen Thyng (see https://matplotlib.org/cmocean/ and https://github.com/matplotlib/cmocean), and were initially developed for Python)Plot one time period example shows how to manipulate an existing output from
rxtracto_3D()orrextractogon()to plot just one time period.Transform the data 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 change 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 example shows how to change the name on the colorbar.
Modify the graph shows how to use the
plotdapfunctionadd_ggplot()to modify a graph once it has been generated.Animate a track shows how to animate a track over time.
Animate a grid shows how to animate a grid with multiple time periods.
The first step is to obtain information about the dataset of interest from the ERDDAP server being used. The needed information include:
- baseURL of the ERDDAP server (the default is https://upwell.pfeg.noaa.gov/erddap)
- the datasetID of the data to be accessed (for example the 1-day composite of chlorophyll from MODIS Aqua has a datasetID of ‘erdMBchla1day’ - see https://upwell.pfeg.noaa.gov/erddap/griddap/erdMBchla1day.html)
- the parameter name (in this case “chlorophyll”)
In order for rerddapXtracto to have this information, as
well as the coordinate variables and their limits, and the parameter
names, a call must be made to the function rerddap::info()
for the appropriate datasetID and baseURL:
require("rerddap")
## base URL does not need to given because it is the default one
dataInfo <- info('erdMBchla1day')
dataInfothis returns the following information:
<ERDDAP info> erdMBchla1day
Base URL: https://upwell.pfeg.noaa.gov/erddap/
Dimensions (range):
time: (2006-01-01T12:00:00Z, 2020-10-27T12:00:00Z)
altitude: (0.0, 0.0)
latitude: (-45.0, 65.0)
longitude: (120.0, 320.0)
Variables:
chlorophyll:
Units: mg m-3An rxtracto example
In this section data is extracted along a trackline found in the Marlintag38606 dataset, which is the track of a tagged marlin in the Pacific Ocean (courtesy of Dr. Mike Musyl of the Pelagic Research Group LLC), and then simple plots of the extracted data are made. Since this can be a long extract, a progress bar is displayed:
require("rerddap")
require("rerddapXtracto")
# First we will copy the Marlintag38606 data into a variable
# called tagData so that subsequent code will be more generic.
tagData <- Marlintag38606
xpos <- tagData$lon
ypos <- tagData$lat
tpos <- tagData$date
zpos <- rep(0., length(xpos))
swchlInfo <- rerddap::info('erdSWchla8day')
swchl1 <- rxtracto(swchlInfo, parameter = 'chlorophyll',
xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos,
xlen = .2, ylen = .2, progress_bar = TRUE)Plotting the results
The track line with the locations colored according to the mean of
the satellite chlorophyll around that point can be plotted
using plotTrack(). 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
cmocean “algae” color palette.
require("ggplot2")
require("plotdap")
myPlot <- plotTrack(swchl1, xpos, ypos, tpos, plotColor = 'algae')
myPlot
Animating the track
To make a cumulative animation of the track:
myPlot <- plotTrack(swchl1, xpos, ypos, tpos, plotColor = 'algae',
animate = TRUE, cumulative = TRUE)
Output from rxtracto() can readily be converted to a
form for use in the “tidyverse” by converting the output to a dataframe,
to wit:
Topography data
This example demonstrates how to pass a function to
plotTrack() to transform the data before plotting, how to
change the name shown on the colorbar, and how to call
plotTrack() if the dataset does not have a time
coordinate.
require("ggplot2")
require("plotdap")
require("rerddap")
require("rerddapXtracto")
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)
myFunc = function(x) -x
topoPlot <- plotTrack(topo, xpos, ypos, NA, plotColor = 'dense',
name = 'Depth', myFunc = myFunc)
topoPlot
Track moving in (x, y, z, t) space.
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.
require("rerddap")
urlBase <- "https://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')
str(sstTrack)This returns the following structure information:
List of 13
$ mean Sea_water_temperature : num [1:3] 11.7 11.7 11.6
$ stdev Sea_water_temperature : num [1:3] 0.0407 0.0555 0.0859
$ n : int [1:3] 493 491 484
$ satellite date : chr [1:3] "2020-10-31T22:00:00Z" "2020-10-31T23:00:00Z" "2020-11-01T00:00:00Z"
$ requested lon min : num [1:3] -10.3 -10.3 -10.3
$ requested lon max : num [1:3] -10.2 -10.2 -10.2
$ requested lat min : num [1:3] 53.5 53.5 53.5
$ requested lat max : num [1:3] 53.5 53.5 53.5
$ requested z min : num [1:3] 2 6 10
$ requested z max : num [1:3] 2 6 10
$ requested date : chr [1:3] "2020-10-31T22:00:00Z" "2020-10-31T23:00:00Z" "2020-11-01T00:00:00Z"
$ median Sea_water_temperature: num [1:3] 11.7 11.7 11.6
$ mad Sea_water_temperature : num [1:3] 0.0208 0.0361 0.0986
- attr(*, "row.names")= chr [1:3] "1" "2" "3"
- attr(*, "class")= chr [1:2] "list" "rxtractoTrack"Crossing the dateline
The following is an artificial example of a track that crosses the date-line, using the MUR(Multi-scale Ultra-high Resolution) SST analysis:
dataInfo <- rerddap::info('jplMURSST41mday')
parameter <- 'sst'
xcoord <- c(179.7, 179.8, 179.9, 180., 180.1, 180.2, 180.3, 180.4)
ycoord <- c(40, 40, 40, 40, 40, 40, 40, 40)
tcoord <- c('2018-03-16', '2018-03-16', '2018-03-16','2018-03-16','2018-03-16',
'2018-03-16','2018-03-16','2018-03-16')
xlen <- .05
ylen <- .05
extract <- rxtracto(dataInfo, parameter = parameter, xcoord = xcoord,
ycoord = ycoord, tcoord = tcoord,
xlen = xlen, ylen = ylen)
str(extract)The returned structure is:
List of 13
$ mean sst : num [1:8] 11.1 11.1 11.1 11.2 11.1 ...
$ stdev sst : num [1:8] 0.01192 0.00602 0.01025 0.00876 0.01446 ...
$ n : int [1:8] 30 30 35 25 30 30 30 35
$ satellite date : chr [1:8] "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" ...
$ requested lon min: num [1:8] 180 180 180 180 180 ...
$ requested lon max: num [1:8] 180 180 180 180 180 ...
$ requested lat min: num [1:8] 40 40 40 40 40 ...
$ requested lat max: num [1:8] 40 40 40 40 40 ...
$ requested z min : logi [1:8] NA NA NA NA NA NA ...
$ requested z max : logi [1:8] NA NA NA NA NA NA ...
$ requested date : chr [1:8] "2018-03-16" "2018-03-16" "2018-03-16" "2018-03-16" ...
$ median sst : num [1:8] 11.1 11.1 11.1 11.2 11.1 ...
$ mad sst : num [1:8] 0.01416 0.0052 0.01149 0.00887 0.01744 ...
- attr(*, "row.names")= chr [1:8] "1" "2" "3" "4" ...
- attr(*, "class")= chr [1:2] "list" "rxtractoTrack"ERDDAP Interpolate service
ERDDAP has a new “Interpolate” service
which takes a table of coordinate values from a gridded dataset, up to
100 entries at a time, and will interpolate the value at each point
based on quite a large list of different interpolation methods, as well
as the number of neighbors to use in the interpolation (see https://coastwatch.pfeg.noaa.gov/erddap/convert/interpolate.html.
As the “Interpolate” service includes the mean, median and standard
deviation of the neighboring points, the service can return similar
results as does the usual rxtracto() but much more quickly,
as 100 observations are processed at a time. This greatly reduces the
time it takes to process very large tracks. A downside is that neighbors
are defined in terms of pixels, not in terms of coordinate units as is
used in the usual rxtracto() call. For very high-resolution
datasets even the largest number of neighbors provided by the service
may not be large enough to capture the uncertainty in the actual
location of an observation in the track.
The “Interpolate” service requires that the ERDDAP server being accessed is at least
version 2.10. As of this writing, the default ERDDAP in rerddap is at https://upwell.pfeg.noaa.gov/erddap, which is not.
Therefore, as in the example below, you need to make certain of the
version of ERDDAP used in the call to
rerddap:info().
The “interpolate” option is implemented in a new, optional argument
to rxtracto(), “interp = NULL”. Since the argument is
optional, and has a default value of “NULL”, existing calls to
rxtracto() will not be affected by the change. To use the
interpolation option, “interp” must be a set equal to a character array
of size 2, that is
For example, to use the ‘Mean’ function with 16 neighbors this would be:
Note carefully that “Mean” is capitalized, and that the number of neighbors is passed as a character string, not as a number. Also the values of “xlen”, “ylen”, and “zlen”, if given, are ignored. So to use the “interp” option in the first Marlin example, the call would be:
my_interp <- c('Mean', '16')
#
# use the coastwatch ERDDAP as it is a more recent version
#
swchlInfo <- rerddap::info('erdSWchla8day',
url = 'https://coastwatch.pfeg.noaa.gov/erddap/')
swchl1 <- rxtracto(swchlInfo, parameter = 'chlorophyll',
xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos,
interp = my_interp, progress_bar = TRUE)Tracks along projected data and Polar Bears
rxtracto() can now extract tracks from projected ERDDAP datasets if the proper data
transformations are done first. Here we cover one possible way of
transforming the data in order to work with rxtracto() by
extracting sea ice concentration along a polar bear track. The polar
bear track is in the included ‘PB_Argos’ dataset, which is described in
the paper “Auger-Méthé, M., Newman, K., Cole, D., Empacher, F., Gryba,
R., King, A. A., Leos-Barajas, V., Mills Flemming, J., Nielsen, A.,
Petris, G., and Thomas, L.. 2021. A guide to state–space modeling of
ecological time series. Ecological Monographs 91(4):e01470.
10.1002/ecm.1470” and is used by permission of the authors. The dataset
was obtained from https://datadryad.org/stash/dataset/doi:10.5061/dryad.4qrfj6q96.
One caveat is that while the vector parameters “xlen”, “ylen” and
“zlen” can still be used with projected data, to be sensible these
values must also be in projected units, and for many projections in
order to have a consistent box the values will have to vary with each
latitude-longitude pair. With the present version of
rxtracto() it is up to the user to make this calculation,
though the steps below should provide some idea of how to do so. In the
example, “xlen”, “ylen” are zero.
Given a “track” with latitude and longitude values and a projected dataset, the following steps are necessary:
Determine the projection used in the ERDDAP dataset.
Transform the latitudes and longitudes into this projection.
Make the extract
If desired, transform the projected coordinates back into latitude and longitude.
The dataset to be used is sea ice concentration from the NOAA/National Snow and Ice Data Center availabe from the Polarwatch ERDDAP at https://polarwatch.noaa.gov/erddap/griddap/nsidcG02202v4nh1day.html .
The first step is to extract the first 100 datapoints from the PB_Argos dataset.
PB_Argos_subset <- PB_Argos[1:100, ]
# datetime is not in format for ERDDAP
PB_Argos_subset$DateTime <- lubridate::as_datetime(PB_Argos_subset$DateTime)
head(PB_Argos_subset, 4)
# A tibble: 4 × 4
DateTime QualClass Lat Lon
<dttm> <chr> <dbl> <dbl>
1 2009-04-20 17:01:39 B 70.4 -132.
2 2009-04-20 17:23:00 A 71.0 -131.
3 2009-04-20 18:12:15 A 71.0 -131.
4 2009-04-20 20:43:17 A 70.9 -131.The dataset info can be used to find the dataset projection (note this uses EPSG codes for the CRS, which can be found at https://spatialreference.org/ref/epsg/):
url= 'https://polarwatch.noaa.gov/erddap/'
datasetid = 'nsidcG02202v4nh1day'
dataInfo <- rerddap::info(datasetid, url)
proj_crs_code_index <- which(dataInfo$alldata$NC_GLOBAL$attribute_name == "proj_crs_code" )
proj_crs_code <- dataInfo$alldata$NC_GLOBAL$value[proj_crs_code_index]
proj_crs_code
[1] "EPSG:3411"Next we will create a temporary dataframe with the latitude and longitude values, convert that dataframe to an ‘sf’ object, and then transform the coordinates:
temp_df <- data.frame(Lat = PB_Argos_subset$Lat, Lon = PB_Argos_subset$Lon)
# transform PB_Argos_subset to sf object
# EPSG:4326 is basic lat-lon coordinates
temp_df <- sf::st_as_sf(temp_df, coords = c("Lon", "Lat"), crs = 4326)
# project data
temp_df <- sf::st_transform(temp_df, crs = 3411)
# get projection coordinates
coordinates <- sf::st_coordinates(temp_df)Now things are set to perform the extract:
xcoord <- coordinates[, 1]
ycoord <- coordinates[, 2]
# R datetimes are passed as number, require ISO character string
tcoord <- as.character(PB_Argos_subset$DateTime)
parameter <- 'cdr_seaice_conc'
extract <- rxtracto(dataInfo,
xName="xgrid",
yName="ygrid",
tName="time",
parameter=parameter,
xcoord = xcoord,
ycoord = ycoord,
tcoord = tcoord
)
str(extract)
List of 13
$ mean cdr_seaice_conc : num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
$ stdev cdr_seaice_conc : num [1:100] NA NA NA NA NA NA NA NA NA NA ...
$ n : int [1:100] 1 1 1 1 1 1 1 1 1 1 ...
$ satellite date : chr [1:100] "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" ...
$ requested x min : num [1:100] -2143958 -2077605 -2077124 -2082576 -2036149 ...
$ requested x max : num [1:100] -2143958 -2077605 -2077124 -2082576 -2036149 ...
$ requested y min : num [1:100] -118477 -131404 -131955 -130368 -120720 ...
$ requested y max : num [1:100] -118477 -131404 -131955 -130368 -120720 ...
$ requested z min : logi [1:100] NA NA NA NA NA NA ...
$ requested z max : logi [1:100] NA NA NA NA NA NA ...
$ requested date : chr [1:100] "2009-04-20 17:01:39" "2009-04-20 17:23:00" "2009-04-20 18:12:15" "2009-04-20 20:43:17" ...
$ median cdr_seaice_conc: num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
$ mad cdr_seaice_conc : num [1:100] 0 0 0 0 0 0 0 0 0 0 ...
- attr(*, "row.names")= chr [1:100] "1" "2" "3" "4" ...
- attr(*, "class")= chr [1:2] "list" "rxtractoTrack"Since a given location may not fall on a gridpoint of the dataset and extracts are performed at the nearest grid point, the locations of the extract in terms of latitude and longitude can be obtained by reversing the steps above:
