Using the xtractomatic routines version 3.4

IMPORTANT - If you have used xtractomatic before

The code has undergone a major refactoring, in particular there is no longer separate code to handle bathymetry and other datasets without time. In order to have a cleaner design to go with these changes, the order of the arguments to the functions has changed, “tpos” no longer has to be passed for a dataset without time, and “xlen” and “ylen” have default values (see next section). Also in xtracto() the return now includes the parameter name, so instead of one element being just “mean” it will be “mean sst” for example. The changes in the function arguments makes them closer to that in rerddapXtracto, see next section.

IMPORTANT - Likely last version of xtractomatic

This will likely be the last version of xtractomatic, except for bug fixes or changes needed to stay compatible with new versions of R, or changes in ERDDAP and the applications it uses. Development effort will be put into improving the rerddapXtracto package, with the intent that users will eventually migrate to rerddapXtracto. As mentioned above, the changes in the function arguments will help to simplify this migration.

The Main xtractomatic functions

There are three main data extraction functions in the xtractomatic package:

• xtracto <- function(dtype, xpos, ypos, tpos = NA, xlen = 0., ylen = 0., verbose=FALSE)

• xtracto_3D <- function(dtype, xpos, ypos, tpos = NA, verbose=FALSE)

• xtractogon <- function(dtype, xpos, ypos, tpos = NA, verbose = FALSE)

There are two information functions in the xtractomatic package:

• `searchData <- function(searchList= “varname:chl”)

• getInfo <- function(dtype)

The dtype parameter in the data extraction routines specifies a combination of which dataset on the ERDDAP server to access, and as well as which parameter from that dataset. The first sections demonstrate how to use these functions in realistic settings. The Selecting a dataset section shows how to find the dtype or other parameters of interest from the available datasets.

Using the R code examples

Besides the xtractomatic packages, the examples below depend on DT, ggplot2, lubridate, mapdata, and reshape2. These can be loaded beforehand (assuming they have been installed):

library("DT")
library("ggplot2")
library("lubridate")
library("mapdata")
library("reshape2")

In order that the code snippets can 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 xtractomatic does not depend on these other packages.

There are also several R functions defined within the document that are used in other code examples. These include chlaAvg, upwell, and plotUpwell.

An xtracto example

In this section we extract data 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 show some simple plots of the extracted data.

The Marlintag38606 dataset is loaded when you load the xtractomatic library. The structure of this dataframe is shown by:

require(xtractomatic)
str(Marlintag38606)

'data.frame':   152 obs. of  7 variables:
 $ date  : Date, format: "2003-04-23" ...
 $ lon   : num  204 204 204 204 204 ...
 $ lat   : num  19.7 19.8 20.4 20.3 20.3 ...
 $ lowLon: num  204 204 204 204 204 ...
 $ higLon: num  204 204 204 204 204 ...
 $ lowLat: num  19.7 18.8 18.8 18.9 18.9 ...
 $ higLat: num  19.7 20.9 21.9 21.7 21.7 ...

The Marlintag38606 dataset consists of three variables, the date, longitude and latitude of the tagged fish, which are used as the xpos, ypos and tpos arguments in xtracto. The parameters xlen and ylen specify the spatial “radius” (actually a box, not a circle) around the point to search in for existing data. These can also be input as a time-dependent vector of values, but here a set value of 0.2 degrees is used. The example extracts SeaWiFS chlorophyll data, using a dataset that is an 8 day composite. This dataset is specified by the dtype of "swchla8day" in the xtracto function. In a later section it is explained how to access different satellite datasets other than the one shown here.

require(xtractomatic)

# 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
swchl <- xtracto("swchla8day", xpos, ypos, tpos = tpos, , xlen = .2, ylen = .2)

This command can take a while to run depending on your internet speed and how busy is the server. With a fairly decent internet connection it took ~25 seconds. The default is to run in quiet mode, if you would prefer to see output, as confirmation that the script is running, use verbose=TRUE

When the extraction has completed the data frame swchl will contain as many columns as data points in the input file (in this case 152) and will have 11 columns:

str(swchl)

'data.frame':   152 obs. of  11 variables:
 $ mean chlorophyll  : num  0.073 NaN 0.074 0.0653 0.0403 ...
 $ stdev chlorophyll : num  NA NA 0.00709 0.00768 0.02278 ...
 $ n                 : int  1 0 16 4 7 9 4 3 0 6 ...
 $ satellite date    : chr  "2003-04-19" "2003-04-27" "2003-04-27" "2003-04-27" ...
 $ requested lon min : num  204 204 204 204 204 ...
 $ requested lon max : num  204 204 204 204 204 ...
 $ requested lat min : num  19.6 19.7 20.3 20.2 20.2 ...
 $ requested lat max : num  19.8 19.9 20.5 20.4 20.4 ...
 $ requested date    : num  12165 12166 12172 12173 12174 ...
 $ median chlorophyll: num  0.073 NA 0.073 0.0645 0.031 ...
 $ mad chlorophyll   : num  0 NA 0.00297 0.00741 0.0089 ...

Plotting the results

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.

require(ggplot2)

Loading required package: ggplot2

require(mapdata)

Loading required package: mapdata
Loading required package: maps

# 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(swchl, tagData)
# 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) * 1
colnames(alldata)[1] <- 'mean'
# set limits of the map
ylim <- c(15, 30)
xlim <- c(-160, -105)
# get outline data for map
w <- map_data("worldHires", ylim = ylim, xlim = xlim)
# plot using ggplot
myColor <- colors$algae
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")

Comparing with the median

The xtracto routine extracts the environment data in a box (defined by xlen and ylen) around the given location, and calculates statistics of the data in that box. The median chla concentration may make more sense in this situation

require(ggplot2)
colnames(alldata)[10] <- 'median'
# plot using ggplot
myColor <- colors$algae
z <- ggplot(alldata,aes(x = lon,y = lat)) + 
   geom_point(aes(colour = median,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(" Median chla values at marlin tag locations")

Topography data

In addition to satellite data, xtractomatic can access topographic data. As an example we extract and plot the water depth (bathymetry) associated with the trackline.

Get the topography data along the track:

require("xtractomatic")
ylim <- c(15, 30)
xlim <- c(-160, -105)
topo <- xtracto("ETOPO180", tagData$lon, tagData$lat, xlen = .1, ylen = .1)

and plot the track:

require("ggplot2")
alldata <- cbind(topo, tagData)
alldata$lon <- alldata$lon - 360
colnames(alldata)[1] <- 'mean'
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("Depth") + 
  coord_fixed(1.3, xlim = xlim, ylim = ylim) + ggtitle("Bathymetry at marlin tag locations")

Comparing chlorophyll estimates from different satellites

This example extracts a time-series of monthly satellite chlorophyll data for the period 1998-2014 from three different satellites ( SeaWIFS, MODIS and VIIRS ), using the xtracto_3D function. The extract is for an area off the coast of California. xtracto_3D extracts data in a 3D bounding box where xpos has the minimum and maximum longitude, ypos has the minimum and maximum latitude, and tpos has the starting and ending time, given as “YYY-MM-DD”, describing the bounding box.

First we define the longitude, latitude and temporal boundaries of the box:

xpos <- c(235, 240)
ypos <- c(36, 39)
tpos <- c('1998-01-01', '2014-11-30') 

The extract will contain data at all of the longitudes, latitudes and times in the requested dataset that are within the given bounds.

If we run the xtracto_3D using the full temporal contraints defined above for SeaWiFS we get an error:

require(xtractomatic)
chlSeaWiFS <- xtracto_3D('swchlamday', xpos, ypos, tpos = tpos )

[1] "tpos  (time) has elements out of range of the dataset"
[1] "time range in tpos"
[1] "1998-01-01,2014-11-30"
[1] "time range in ERDDAP data"
[1] "1997-09-16,2010-12-16 12:00:00"

Error in xtracto_3D("swchlamday", xpos, ypos, tpos = tpos) : 
  Coordinates out of dataset bounds - see messages above

The error occurs because the time-range specified is outside of the time-bounds of the requested dataset ( SeaWIFS data stopped in 2010). The error message generated by xtracto_3D explains this, and gives the temporal boundaries of the dataset. A simple way to extract to the end of a dataset is to use “last” (see Taking advantage of “last times”):

require(xtractomatic)
tpos <- c("1998-01-16", "last")
SeaWiFS <- xtracto_3D('swchlamday', xpos, ypos, tpos = tpos)

A similar extract can be made for the MODIS dataset. We know that the time-range of the MODIS dataset is also smaller than that defined by tpos. We can determine the exact bounds of the MODIS dataset by using the getInfo function:

require(xtractomatic)
getInfo('mhchlamday')

List of 19
 $ dtypename       : chr "mhchlamday"
 $ datasetname     : chr "erdMH1chlamday"
 $ longname        : chr "Chlorophyll-a, Aqua MODIS, NPP, L3SMI, Global, Science Quality (Monthly Composite)"
 $ varname         : chr "chlorophyll"
 $ hasAlt          : logi FALSE
 $ latSouth        : logi FALSE
 $ lon360          : logi FALSE
 $ minLongitude    : num -180
 $ maxLongitude    : num 180
 $ longitudeSpacing: num 0.0417
 $ minLatitude     : num -90
 $ maxLatitude     : num 90
 $ latitudeSpacing : num -0.0417
 $ minAltitude     : num NA
 $ maxAltitude     : num NA
 $ minTime         : chr "2003-01-16"
 $ maxTime         : chr "2017-11-16"
 $ timeSpacing     : num NA
 $ infoURL         : chr "http://oceancolor.gsfc.nasa.gov/"

which gives us the min and max time of the dataset, which we can be used in the call (or use “last” as above):

require(xtractomatic)
tpos <- c('2003-01-16', "last")
MODIS <- xtracto_3D('mhchlamday', xpos, ypos, tpos = tpos)

Similarly we can extract data from the VIIRS dataset.

require(xtractomatic)
tpos <- c("2012-01-15", "last")
VIIRS <- xtracto_3D('erdVH3chlamday', xpos, ypos, tpos = tpos)

xtracto_3d returns a list of the form:

• extract$data : num [longitude, latitude, time]
• extract$varname : character string with the name of the variable
• extract$datasetname: character string of the ERDDAP dataset ID
• extract$latitude : num [latitude] latitudes of extract
• extract$longitude : num [longitude] longitude of extact
• extract$time : chr [time] times of extract
• extract$altitude : num altitude of extract

in this case the varname is chla, the ERDDAP datasetID is erdVH3chlamday, and data contains the data on the grid defined by longitude, latitude and time. The data from the different satellites can be compared by mapping the values for a given time period, which we cna do using the function plotBBox().

mapFrame <- function(longitude, latitude, chla){
  dims <- dim(chla)
  chla <- array(chla, dims[1] * dims[2])
  longitude <- longitude - 360
  chlaFrame <- expand.grid(x = longitude, y = latitude)
  chlaFrame$chla <- chla
  return(chlaFrame)
}

and also define a helper function plotFrame to plot the data:

plotFrame <- function(chlaFrame, xlim, ylim, title, logplot = TRUE){
  require(mapdata)
  require(ggplot2)
  w <- map_data("worldHires", ylim = ylim, xlim = xlim)
  myplot <- ggplot(data = chlaFrame, aes(x = x, y = y, fill = chla)) +
    geom_raster(interpolate = FALSE, na.rm = TRUE) +
    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$algae
    myplot <- myplot + scale_fill_gradientn(colours = my.col, na.value = NA, limits = c(-2, 4.5)) +
      ggtitle(title)
    }else{
      my.col <- colors$algae
     myplot <- myplot + scale_fill_gradientn(colours = my.col, na.value = NA) +
       ggtitle(title)
  }
  return(myplot)
}

We examine chlorophyll in June 2012 a period when VIIRS and MODIS both had data. Starting with VIIRS:

require(lubridate)

Loading required package: lubridate

Attaching package: ‘lubridate’

The following object is masked from ‘package:base’:

    date

xlim <- c(235, 240) - 360
ylim <- c(36, 39)
ttext <- VIIRS$time[month(VIIRS$time) == 6 & year(VIIRS$time) == 2012]
chlaFrame <- mapFrame(VIIRS$longitude,VIIRS$latitude,VIIRS$data[, , month(VIIRS$time) == 6 & year(VIIRS$time) == 2012])
chlaPlot <- plotFrame(chlaFrame, xlim, ylim, paste("VIIRS chla", ttext), logplot = FALSE)
chlaPlot

Chlorophyll values are highly skewed, with low values most places and times, and then for some locations very high values. For this reason log(chlorophyll) is usually plotted.

xlim <- c(235, 240) - 360
ylim <- c(36, 39)
chlalogFrame <- mapFrame(VIIRS$longitude,VIIRS$latitude,
                       log(VIIRS$data[, , month(VIIRS$time) == 6 & year(VIIRS$time) == 2012]))
chlalogPlot <- plotFrame(chlalogFrame, xlim, ylim, paste("VIIRS log(chla)", ttext))
chlalogPlot

And also for MODIS:

xlim <- c(235, 240) - 360
ylim <- c(36, 39)
ttext <- MODIS$time[month(MODIS$time) == 6 & year(MODIS$time) == 2012]
chlalogFrame <- mapFrame(MODIS$longitude, MODIS$latitude,
                       log(MODIS$data[, , month(MODIS$time) == 6 & year(MODIS$time) == 2012]))
chlalogPlot <- plotFrame(chlalogFrame, xlim, ylim, paste("MODIS log(chla)", ttext))
chlalogPlot

SeaWiFS stopped functioning at the end of 2010, so we look at June, 2010:

xlim <- c(235, 240) - 360
ylim <- c(36, 39)
ttext <- SeaWiFS$time[month(SeaWiFS$time) == 6 & year(SeaWiFS$time) == 2010]
chlalogFrame <- mapFrame(SeaWiFS$longitude, SeaWiFS$latitude,
                       log(SeaWiFS$data[, , month(SeaWiFS$time) == 6 & year(SeaWiFS$time) == 2010]))
chlalogPlot <- plotFrame(chlalogFrame, xlim, ylim, paste("SeaWiFS log(chla)", ttext))
chlalogPlot

We can also look at changes in log(chlorophyll) over time for a given region. We examine a higly productive region just north of San Francisco, with limits (38N, 38.5N) and (123.5W, 123W).

require(ggplot2)
require(mapdata)
xlim1 <- c(-123.5, -123) + 360
ylim1 <- c(38, 38.5)
w <- map_data("worldHires", ylim = c(37.5, 39), xlim = c(-123.5, -122.))
z <- ggplot() + geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
   theme_bw() +
   coord_fixed(1.3, xlim = c(-123.5, -122.), ylim = c(37.5, 39))
z + geom_rect(aes(xmin = -123.5, xmax = -123., ymin = 38., ymax = 38.5), colour = "black", alpha = 0.) +
  theme_bw() + ggtitle("Location of chla averaging")

We define a helper function chlaAvg which finds the latitudes and longitudes in the datasets that are closest to the given limits, extracts the subset that is within that region, and then averages over the region to create a time series.

chlaAvg <- function(longitude, latitude, chlaData, xlim, ylim, sTimes){
  xIndex <- xlim
  yIndex <- ylim
  yIndex[1] <- which.min(abs(latitude - ylim[1]))
  yIndex[2] <- which.min(abs(latitude - ylim[2]))
  xIndex[1] <- which.min(abs(longitude - xlim[1]))
  xIndex[2] <- which.min(abs(longitude - xlim[2]))
  tempData <- chlaData[xIndex[1]:xIndex[2], yIndex[1]:yIndex[2],]
  chlaAvg <- apply(tempData, 3, mean, na.rm = TRUE)
  chlaAvg <- data.frame(chla = chlaAvg, time = sTimes)
return(chlaAvg)
}

We use chlaAvg to extract the three time series over the region, combine them and plot the result.

xlim1 <- c(-123.5, -123) + 360
ylim1 <- c(38, 38.5)
# Get both the average, and the average of log transformed chl for each point in the time series 
SeaWiFSavg <- chlaAvg(SeaWiFS$longitude, SeaWiFS$latitude, SeaWiFS$data, xlim1, ylim1,SeaWiFS$time)
SeaWiFSlog <- chlaAvg(SeaWiFS$longitude, SeaWiFS$latitude, log(SeaWiFS$data), xlim1, ylim1, SeaWiFS$time)
# Run the same steps again for the MODIS and VIIRS datasets
MODISavg <- chlaAvg(MODIS$longitude, MODIS$latitude, MODIS$data, xlim1, ylim1, MODIS$time)
MODISlog <- chlaAvg(MODIS$longitude, MODIS$latitude, log(MODIS$data), xlim1, ylim1, MODIS$time)
# run the same for VIIRS
VIIRSavg <- chlaAvg(VIIRS$longitude, VIIRS$latitude, VIIRS$data, xlim1, ylim1, VIIRS$time)
VIIRSlog <- chlaAvg(VIIRS$longitude, VIIRS$latitude, log(VIIRS$data), xlim1, ylim1, VIIRS$time)

First the raw chla values:

require(ggplot2)
Chla <- rbind(VIIRSavg, MODISavg, SeaWiFSavg)
Chla$sat <- c(rep("VIIRS", nrow(VIIRSavg)), rep("MODIS", nrow(MODISavg)), rep("SeaWIF", nrow(SeaWiFSavg)))
Chla$sat <- as.factor(Chla$sat)
ggplot(data = Chla, aes(time, chla, colour = sat)) + geom_line(na.rm = TRUE) + theme_bw()

and then log(chla):

require(ggplot2)
logChla <- rbind(VIIRSlog, MODISlog, SeaWiFSlog)
logChla$sat <- c(rep("VIIRS", nrow(VIIRSlog)), rep("MODIS", nrow(MODISlog)), rep("SeaWIF",    nrow(SeaWiFSlog)))
logChla$sat <- as.factor(Chla$sat)
ggplot(data = logChla, aes(time, chla, colour = sat)) + geom_line(na.rm = TRUE)  + theme_bw()

Upwelling

ERD for many years has produced 6-hourly Bakun upwelling indices at 15 set locations. But suppose you want the values at a different location? To do this we define a function upwell that rotates Ekman transport values to obtain the offshore component (upwelling index), and also define a function plotUpwell to plot upwelling. xtracto_3D can download Ekman transport calculated from Fleet Numerical Meteorlogical and Oceanographic Center (FNMOC) pressure fields:

upwell <- function(ektrx, ektry, coast_angle){
   pi <- 3.1415927
   degtorad <- pi/180.
   alpha <- (360 - coast_angle) * degtorad
   s1 <- cos(alpha)
   t1 <- sin(alpha)
   s2 <- -1 * t1
   t2 <- s1
   perp <- s1 * ektrx + t1 * ektry
   para <- s2 * ektrx + t2 * ektry
return(perp/10)
}

plotUpwell <- function(upwelling, upDates){
require(ggplot2)
temp <- data.frame(upwelling = upwelling, time = upDates)
ggplot(temp, aes(time, upwelling)) + geom_line(na.rm = TRUE) + theme_bw()
}

We calculate 6-hourly upwelling at (37,-122) for the year 2005 and use the coast angle that ERD uses for (36N , 122W) which is 152 degrees.

require(xtractomatic)
xpos <- c(238, 238)
ypos <- c(37, 37)
tpos <- c("2005-01-01", "2005-12-31")
ektrx <- xtracto_3D("erdlasFnTran6ektrx", xpos, ypos, tpos = tpos)
ektry <- xtracto_3D("erdlasFnTran6ektry", xpos, ypos, tpos = tpos)
upwelling <- upwell(ektrx$data, ektry$data, 152)

A plot of the result:

plotUpwell(upwelling, as.Date(ektrx$time))

The one very low downwelling value distorts the plot, so values < -500 are set to NA and the data are replotted:

tempUpw <- upwelling
tempUpw[tempUpw < -500] <- NA
plotUpwell(tempUpw, as.Date(ektrx$time))

In 2005 there was a large die-off of sea animals and many felt that is was due to an anomalously delayed upwelling. Let’s compare upwelling in 2005 at that location with upwelling in 1977:

require(xtractomatic)
tpos <- c("1977-01-01", "1977-12-31")
ektrx77 <- xtracto_3D("erdlasFnTran6ektrx", xpos, ypos, tpos = tpos)
ektry77 <- xtracto_3D("erdlasFnTran6ektry", xpos, ypos, tpos = tpos)
upwelling77 <- upwell(ektrx77$data, ektry77$data, 152)
# remove the one big storm at the end
tempUpw <- upwelling77
tempUpw[tempUpw < -500] <- NA
plotUpwell(tempUpw, as.Date(ektrx77$time))

While we have not looked at all areas along the coast, at least for this location it is difficult to see upwelling in 2005 as being any more delayed than in 1977 when there was not the same affect on animals.

Comparisons of Wind Stress

A problem faced when using satellite based estimates of wind stress is that the QuikSCAT based estimates start in 1999-07-21 and end in 2009-11-19, while the ASCAT based estimates start in 2009-10-03 and are ongoing. A comparison of the data over the period of overlap would be of interest.

We examine north-south wind stress (tauy) over the period of overlap at (36N, 121W) (we go a little more offshore because the satellite data is not resolved close to the coast), and compare 3-day composites of tauy. In this region tauy is the main component of upwelling as the coastline runs almost north-south and the winds are mostly from the north or north-west.

require(xtractomatic)
xpos <- c(237, 237)
ypos <- c(36, 36)
tpos <- c("2009-10-04", "2009-11-19")
ascat <- xtracto_3D("erdQAstress3daytauy", xpos, ypos, tpos = tpos)
quikscat <- xtracto_3D("erdQSstress3daytauy", xpos, ypos, tpos = tpos)

Plotting the result

require(ggplot2)
ascatStress <- data.frame(stress = ascat$data, time = as.Date(ascat$time))
quikscatStress <- data.frame(stress = quikscat$data, time = as.Date(quikscat$time))
stressCompare <- rbind(ascatStress,quikscatStress)
stressCompare$sat <- c(rep("ascat", nrow(ascatStress)), rep("quikscat", nrow(quikscatStress)))
stressCompare$sat <- as.factor(stressCompare$sat)
ggplot(data = stressCompare, aes(time, stress, colour = sat)) + geom_line(na.rm = TRUE) + theme_bw()