---
title: "Using the xtractomatic routines"
author: "Cara Wilson and Roy Mendelssohn"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Usingxtractomatic}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
```{r initialize, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(xtractomatic)
```

## Introduction 

`xtractomatic` is an <span style="color:blue">R</span> package developed to subset and extract satellite and other oceanographic related data from a remote 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).  The `xtractomatic` functions were originally developed for the marine biology tagging community, to match up environmental data available from satellites (sea-surface temperature, sea-surface chlorophyll, sea-surface height, sea-surface salinity, vector winds) to track data from various tagged animals or shiptracks (`xtracto`). The package has since been extended to include the routines that extract data a 3D bounding box (`xtracto_3D`) or within a polygon (`xtractogon`).  The `xtractomatic`  package accesses  data that are served through the <span style="color:blue">ERDDAP</span> (Environmental Research Division Data Access Program) server at the NOAA/SWFSC Environmental Research Division in Santa Cruz, California. The <span style="color:blue">ERDDAP</span> server can also be directly accessed at <http://coastwatch.pfeg.noaa.gov/erddap>. <span style="color:blue">ERDDAP</span> is a simple to use yet powerful web data service developed by Bob Simons.  

This is the third version of `xtractomatic`.  The original version was created by Dave Foley for <span style="color:blue">Matlab</span> and was modified for <span style="color:blue">R</span> by Cindy Bessey and Dave Foley.  That version used a precurser to <span style="color:blue">ERDDAP</span> called <span style="color:blue">BobDAP</span>.  <span style="color:blue">BobDAP</span> was less flexible and provided access to fewer datasets than does <span style="color:blue">ERDDAP</span>.  The second version of xtractomatic was written by Roy Mendelssohn to use <span style="color:blue">ERDDAP</span> rather than <span style="color:blue">BobDAP</span>.  This version, also written by Roy Mendelssohn, has many improvements including access to many more datasets,  improved speed, better error checking, and crude search capabilites, and is the first version as a true <span style="color:blue">R</span> package.

### If you have used xtractomatic before

If you have existing <span style="color:blue">R</span> code that uses older versions of the `xtractomatic` functions, the code will likely require some simple and straightforward modifications to be used with this version.  See the section [What has Changed and Steps Going Forward](#new)

### The Main xtractomatic functions

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

- `xtracto <- function(xpos, ypos, tpos, dtype, xlen, ylen, verbose=FALSE)`

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

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


There also are two information functions in the `xtractomatic` package: 

- `searchData <- function(searchList=list(list("varname","chl"))) ` 

- `getInfo <- function(dtype)`

The `dtype` parameter in the data extraction routines specifies a combination of which dataset on the <span style="color:blue">ERDDAP</span> 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](#dataset) section shows how to find the `dtype` or other parameters of interest from the available datasets.

## Setting up 

`xtractomatic` uses the `httr`, `ncdf4` and `sp` packages, and these packages must be installed first or `xtractomatic` will fail to install (the Windows version of `ncdf4` is not available through CRAN because of some quirks in CRAN policy, but it can be obtained at <http://cirrus.ucsd.edu/~pierce/ncdf/>).   

```{r install,eval=FALSE}
install.packages("httr", dependencies = TRUE)
install.packages("ncdf4",dependencies = TRUE) 
install.packages("sp", dependencies = TRUE)
```

The `xtractomatic` package is available from <span style="color:blue">CRAN</span> and can be installed by:

```{r installCRAN,eval=FALSE}
install.packages("xtractomatic")
```


or the development version is available from [Github](https://github.com/rmendels/xtractomatic and can be installed from <span style="color:blue">Github</span>,

```{r installGit,eval=FALSE}
install.packages("devtools")
devtools::install_github("rmendels/xtractomatic")
```

Once installed, to use `xtractomatic` do

```{r load,eval=FALSE}
library("xtractomatic")
```

If the other libraries  (`httr`, `ncdf4`, and `sp`) have been installed they will be found and do not need to be explicitly loaded.

### Using the R code examples

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

```{r loadpacks,eval=FALSE}
library("ggplot2")
library("ggfortify") 
library("lubridate")
library("mapdata")
library("RColorBrewer")
library("reshape2")
library("xts")
```

In order that the code snippets can be more stand-alone, the needed libraries are always `required` in the examples.  The package `ggfortify` at the time of this writing is not on `CRAN`.  It can be installed using the `devtools` package:

```{r ggfortify,eval=FALSE}
devtools::install_github("sinhrks/ggfortify")
```

It should be emphasized that these other packages are used to manipulate and plot the data in <span style="color:blue">R</span>, other packages could be used as well.  The use of `xtractomatic` does not depend on these other packages.

There are also several <span style="color:blue">R</span> functions defined within the document that are used in other code examples.  These include `mapFrame`, `plotFrame`, `chlaAvg`, `upwell`, and `plotUpwell`.

The plots use `ggplot2`.  A basic tutorial on `ggplot2` is given by Dr. Eric Anderson of the National Marine Fisheries Service at [ggplot2_lecture1](http://eriqande.github.io/rep-res-web/lectures/ggplot_1_pole_vault_example.html#ggplot2-intro-lecture) and [ggplot2_lecture2](http://eriqande.github.io/rep-res-web/lectures/ggplot_2_reshape_facets_stats.html#ggplot2-second-lecture). There also is the tutorial at [ggplot2_Handout](http://www.ceb-institute.org/bbs/wp-content/uploads/2011/09/handout_ggplot2.pdf) though it is getting a little dated for some features of `ggplot2`. Anderson also has a simple [guide to maps in R](http://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html#map-making-in-R), other guides on making maps using `ggplot2` include  [spatial.ly](http://spatial.ly/2013/12/introduction-spatial-data-ggplot2/) and [zevross](http://zevross.com/blog/2014/07/16/mapping-in-r-using-the-ggplot2-package/).  Eric Anderson also has a [nice page](http://eriqande.github.io/2014/12/24/plotting-sst-with-ggplot.html) on combining raster data with maps using `ggplot2`.


## Getting Started

### An xtracto example 

In this section we extract data along a trackline found in the <span style="color:red">Marlintag38606</span> 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 <span style="color:red">Marlintag38606</span> dataset is loaded when you load the `xtractomatic` library. The structure of this dataframe is shown by: 

```{r MarlinTag}
require(xtractomatic)
str(Marlintag38606)
```

The  <span style="color:red">Marlintag38606</span> dataset consists of three variables, the <span style="color:blue">date</span>, <span style="color:blue">longitude</span> and <span style="color:blue">latitude</span> of the tagged fish, which are used as the <span style="color:blue">xpos</span>, <span style="color:blue">ypos</span> and <span style="color:blue">tpos</span> arguments in `xtracto`. The parameters <span style="color:blue">xlen</span> and <span style="color:blue">ylen</span> 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 <span style="color:blue">SeaWiFS chlorophyll data</span>, using a dataset that is an 8 day composite. This dataset is specified by the <span style="color:blue">dtype</span> 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.  


```{r getMarlinChl}
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(xpos, ypos, tpos, "swchla8day", 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 extaction 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: 

```{r swchlStr}
str(swchl)
```

### Plotting the results

We plot the track line with the locations colored according to the mean of the satellite <span style="color:blue">chlorophyll</span> around that point. Positions where there was a tag location but no <span style="color:blue">chlorophyll</span> values are also shown.

```{r meantrackPlot, fig.align='center',fig.width=6,fig.height = 4}
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)*1
# 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
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_gradient(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 <span style="color:blue">median chla</span> concentration may make more sense in this situation

```{r mediantrackPlot, fig.align='center',fig.width=6,fig.height = 4}
require(ggplot2)
# plot using ggplot
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_gradient(limits=c(0.,0.32),"Chla") + 
  coord_fixed(1.3,xlim = xlim, ylim = ylim) + ggtitle(" Median chla values at marlin tag locations")
```

The <span style="color:blue">median chla</span> is smaller than the <span style="color:blue">mean chla</span> (a maximum of  0.273 compared to 0.323).

### 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.

```{r topotag, fig.align='center',fig.width=6,fig.height = 4, warning = FALSE}
require("xtractomatic")
require("ggplot2")
ylim<-c(15,30)
xlim<-c(-160,-105)
topo <- xtracto(tagData$lon,tagData$lat,tagData$date,"ETOPO360",.1,.1)
alldata <- cbind(tagData, topo)
alldata$lon<-alldata$lon-360
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")

```


## Reading in tagged data 

The above example uses data already converted to <span style="color:blue">R</span> format, but clearly the utility of the `xtracto` function lies in using it with one's own dataset.  This dataset was read in from a text file called <span style="color:blue">Marlin-tag38606.txt</span>, which is also made available when the `xtractomatic` package is loaded.  To be able to refer to this file: 

```{r marlinTxt,results = "hide"}

system.file("extdata", "Marlin-tag38606.txt", package = "xtractomatic")

```

The first 5 lines of the file can be viewed: 
```{r}
datafile <- system.file("extdata", "Marlin-tag38606.txt", package = "xtractomatic")
system(paste("head -n5 ",datafile))

```

and read in as follows: 

```{r}

Marlingtag38606 <- read.csv(datafile, head=TRUE, stringsAsFactors = FALSE,sep="\t")
str(Marlingtag38606)
```

The file Marlingtag38606.txt is a "tab" separated file, not a comma seperated file, so the `sep` option in the `read.csv` function above indicates that to the function, and `stringsAsFactors = FALSE` tells the function not to treat the dates as factors. 

To be useful the date field, now formatted for example as `4/23/2003`, needs to be converted to an  <span style="color:blue">R</span> date field:

```{r,eval=FALSE}
Marlintag38606$date <- as.Date(Marlintag38606$date, format='%m/%d/%Y')
```

The format of the date has to be given because the dates in the file are not in one of the formats recognized by `as.Date`. A small y (%y) indicates the year has 2 digits (97) while a capital Y (%Y) indicates the year has 4 digits (1997). If your file is in the EU format with commas as decimal points and semicolons as field separators then use read.csv2 rather than read.csv.   

This dataset includes the range of the 95% confidence interval for both the latitude and the longitude, which are the last four columns of data.  In the example above the search "radius" was read into xtracto was a constant value of 0.2 degrees, but a vector can also be read in.  To use the limits given in the file, create the vectors as follows:

```{r,eval=FALSE}
xrad <- abs(Marlintag38606$higLon - Marlintag38606$lowLon)
yrad <- abs(Marlintag38606$higLat - Marlintag38606$lowLat)
```

and then use these for the `xlen` and `ylen` inputs in `xtracto`: 

```{r,eval=FALSE}
xpos<-Marlintag38606$lon 
ypos<-Marlintag38606$lat 
tpos<-Marlintag38606$date
swchl <- xtracto(xpos, ypos, tpos, "swchla8day", xrad, yrad)
```

## Using `xtracto_3D` 

### Comparing chlorophyll estimates from different satellites

This example extracts a time-series of monthly satellite <span style="color:blue">chlorophyll</span> data for the period 1998-2014 from three different satellites ( <span style="color:red">SeaWIFS</span>,  <span style="color:red">MODIS</span> and <span style="color:red">VIIRS</span> ), 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 <span style="color:blue">longitude</span>,  `ypos` has the minimum and maximum <span style="color:blue">latitude</span>, and `tpos` has the starting and ending <span style="color:blue">time</span>, given as "YYY-MM-DD", describing the bounding box.

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

```{r setLims3D} 
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 <span style="color:red">SeaWiFS</span> we get an error: 

```{r SeaWiFSBad, error=TRUE, purl = FALSE}
require(xtractomatic)
chlSeaWiFS<-xtracto_3D(xpos, ypos, tpos, 'swchlamday')
```

The error occurs because the time-range specified is outside of the time-bounds of the requested dataset ( <span style="color:red">SeaWIFS</span> 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"](#lastTimes)): 


```{r SeaWiFS}
require(xtractomatic)
require(lubridate)
tpos<-c("1998-01-16","last")
SeaWiFS<-xtracto_3D(xpos, ypos, tpos, 'swchlamday')
SeaWiFS$time<-as.Date(SeaWiFS$time)
```

A similar extract can be made for the <span style="color:red">MODIS</span> dataset. We know that the time-range of the <span style="color:red">MODIS</span> dataset is also smaller than that defined by `tpos`. We can determine the exact bounds of the <span style="color:red">MODIS</span> dataset by using the `getInfo` function: 

```{r MODISgetInfo}
require(xtractomatic)
getInfo('mhchlamday')
```

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

```{r MODISchla,error=TRUE, purl = FALSE}
require(xtractomatic)
require(lubridate)
tpos<-c('2003-01-16',"last")
MODIS<-xtracto_3D(xpos, ypos, tpos,'mhchlamday')
MODIS$time<-as.Date(MODIS$time)
```

Similarly we can extract data from the <span style="color:red">VIIRS</span> dataset. 

```{r VIIRSchla}
require(xtractomatic)
require(lubridate)
tpos<-c("2012-01-15","last")
VIIRS<-xtracto_3D(xpos, ypos, tpos, 'erdVH2chlamday')
VIIRS$time <- as.Date(VIIRS$time)
```

`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 <span style="color:blue">varname</span> is <span style="color:blue">chla</span>, the <span style="color:blue">ERDDAP</span> `datasetID` is <span style="color:blue">erdVH2chlamday</span>, and <span style="color:blue">data</span> contains the data on the grid defined by <span style="color:blue">longitude</span>, <span style="color:blue">latitude</span> and <span style="color:blue">time</span>.  The data from the different satellites can be compared by mapping the values for a given time period. To do so we define a helper function `mapFrame` to reshape the data to be used in `ggplot2`.

```{r mapFrame}
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:

```{r plotFrame}
plotFrame<-function(chlaFrame, xlim, ylim,title,logplot=TRUE){
  require(mapdata)
  require(ggplot2)
  require(RColorBrewer)
  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 <- colorRampPalette(rev(brewer.pal(11, "RdBu")))(5.5)  
    myplot<-myplot + scale_fill_gradientn(colours = my.col, na.value = NA,limits=c(-2,4.5)) +
      ggtitle(title)
    }else{
     my.col <- colorRampPalette(rev(brewer.pal(11, "RdBu")))((diff(range(chlaFrame$chla,na.rm=TRUE)))) 
     myplot<-myplot + scale_fill_gradientn(colours = my.col, na.value = NA) +
       ggtitle(title)
  }
  return(myplot)
}
```

We examine <span style="color:blue">chlorophyll</span> in June 2012 a period when VIIRS and MODIS both had data.  Starting with <span style="color:red">VIIRS</span>:

```{r VIIRSPlot, fig.width=5,fig.height = 5, fig.align='center'}
require(lubridate)
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
```

<span style="color:blue">Chlorophyll</span> values are highly skewed, with low values most places and times, and then for some locations very high values.  For this reason <span style="color:blue">log(chlorophyll)</span> is usually plotted.

```{r VIIRSLogPlot, fig.width=5,fig.height = 5, fig.align='center'}
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 <span style="color:red">MODIS</span>:

```{r MODISLogPlot, fig.width=5,fig.height = 5, fig.align='center'}
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
```

<span style="color:red">SeaWiFS</span> stopped functioning at the end of 2010, so we look at June, 2010:

```{r SeaWiFSLogPlot, fig.width=5,fig.height = 5, fig.align='center'}
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 <span style="color:blue">log(chlorophyll)</span> 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).

```{r chlaAvgset, fig.width=5,fig.height = 3, fig.align='center'}
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 <span style="color:blue">latitudes</span> and <span style="color:blue">longitudes</span> 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.

```{r chlaAvg}
chlaAvg<-function(longitude,latitude,chlaData,xlim,ylim,dStart){
  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<-ts(chlaAvg, frequency = 12, start = dStart)
return(chlaAvg)
}

```



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

```{r AverageChla}
require(ggplot2)
require(ggfortify)
# 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,c(1998,1))
SeaWiFSlog <- chlaAvg(SeaWiFS$longitude,SeaWiFS$latitude,log(SeaWiFS$data),xlim1,ylim1,c(1998,1))
# Run the same steps again for the MODIS and VIIRS datasets
MODISavg <- chlaAvg(MODIS$longitude,MODIS$latitude,MODIS$data,xlim1,ylim1,c(2003,1))
MODISlog <- chlaAvg(MODIS$longitude,MODIS$latitude,log(MODIS$data),xlim1,ylim1,c(2003,1))
# run the same for VIIRS
VIIRSavg <- chlaAvg(VIIRS$longitude,VIIRS$latitude,VIIRS$data,xlim1,ylim1,c(2012,1))
VIIRSlog <- chlaAvg(VIIRS$longitude,VIIRS$latitude,log(VIIRS$data),xlim1,ylim1,c(2012,1))
```

First the raw <span style="color:blue">chla</span> values:

```{r plotChlaTS, fig.width=8,fig.height = 4, fig.align='center',warning=FALSE}
require(ggplot2)
require(ggfortify)
Chla<-cbind(VIIRSavg,MODISavg,SeaWiFSavg)
autoplot(Chla, facets = FALSE) + theme_bw() + ggtitle("Chla time series off Point Reyes")
```

and then <span style="color:blue">log(chla)</span>:

```{r plotlogChlaTS, fig.width=8,fig.height = 4, fig.align='center',warning=FALSE}
require(ggplot2)
require(ggfortify)
logChla<-cbind(VIIRSlog,MODISlog,SeaWiFSlog)
autoplot(logChla, facets = FALSE) + theme_bw() + ggtitle("log(Chla) time series off Point Reyes")
```


### Upwelling

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


```{r upCalc}
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)
}
```

```{r plotUpwellingFun}
plotUpwell<-function(upwelling,upDates){
require(ggplot2)
require(ggfortify)
require(xts)
temp<-xts(upwelling,order.by=upDates)
autoplot(temp) + theme_bw() + ylab("Upwelling")
}
```

We calculate 6-hourly upwelling at (37,-122) for the year 2005 and use the coast angle that <span style="color:red">ERD</span> uses for (36N , 122W) which is 152 degrees.

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

A plot of the result:

```{r plotUpw1, fig.width=7,fig.height = 4, fig.align='center',warning=FALSE}
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:
 
```{r plotUpw3, fig.width=7,fig.height = 4, fig.align='center',warning=FALSE}
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 <span style="color:blue">upwelling</span>.  Let's compare <span style="color:blue">upwelling</span> in 2005 at that location with <span style="color:blue">upwelling</span> in 1977:

```{r plotUpw77, fig.width=7,fig.height = 4, fig.align='center',warning=FALSE}
require(xtractomatic)
tpos<-c("1977-01-01","1977-12-31")
ektrx77<-xtracto_3D(xpos,ypos,tpos,"erdlasFnTran6ektrx")
ektry77<-xtracto_3D(xpos,ypos,tpos,"erdlasFnTran6ektry")
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 <span style="color:blue">wind stress</span> is that the <span style="color:red">QuikSCAT</span> based estimates start in 1999-07-21 and end in 2009-11-19, while the <span style="color:red">ASCAT</span> 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 <span style="color:blue">north-south wind stress (tauy)</span> 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 <span style="color:blue">tauy</span>. In this region <span style="color:blue">tauy</span> is the main component of <span style="color:blue">upwelling</span> as the coastline runs almost north-south and the winds are mostly from the north or north-west.

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

Plotting the result

```{r plotstress, fig.width=6,fig.height = 4, fig.align='center',warning=FALSE}
require(ggplot2)
require(ggfortify)
require(xts)
ascatStress<-xts(ascat$data, order.by = as.Date(ascat$time))
quikscatStress<-xts(quikscat$data, order.by = as.Date(quikscat$time))
stressCompare<-cbind(ascatStress,quikscatStress)
colnames(stressCompare)<-c("ascat","quikscat")
autoplot(stressCompare, facets = FALSE) + theme_bw()
```

While the two <span style="color:blue">wind stress</span> series are close, there is a significant difference in the negative (southward) value on October 28, that is there will be a signifcant difference in the estimated (positive) upwelling value.  Combining the series into a longer series would be desirable, but this suggests that there could be problems with inter-year comparisons due to the different sources of the data.

### Taking advantage of "last times" {#lastTimes}

The <span style="color:blue">ERDDAP</span> server understands a time of "last" as being the last time period available for that dataset, and permits arithmetic using "last", that is "last-5" is a valid time and will be 5 time periods before the last time (note that since different datasets have different time periods this difference is dataset dependent). For example to extract the last six months of <span style="color:blue">SST</span> from <span style="color:red">POES AVHRR</span> we  look at the monthly <span style="color:blue">agsstmday</span> dataset:

```{r getPoes}
require(xtractomatic)
xpos<-c(235,240)
ypos<-c(36,39)
tpos<-c("last-5","last")
poes<-xtracto_3D(xpos,ypos,tpos,"agsstamday")

```

Plotting the result using faceting in ggplot2:

```{r plotPoes,fig.width=7.5, fig.height = 4, fig.align='center'}
require(ggplot2)
require(reshape2)
require(mapdata)

melt_sst <- function(longitude,latitude,mTime,sst) {
  dimnames(sst) <- list(long = longitude, lat = latitude)
  ret <- melt(sst, value.name = "sst")
  cbind(date = mTime, ret)
}
xpos<-xpos-360
longitude<-poes$longitude-360
latitude<-poes$latitude
imon<-1:6
tmp <- lapply(imon, function(x) melt_sst(longitude,latitude,poes$time[x],poes$data[,,x]))
allmons <- do.call(rbind, tmp)
w <- map_data("worldHires", ylim = ypos, xlim = xpos)
ggplot(data = allmons, aes(x = long, y = lat, fill = sst)) + 
    geom_polygon(data = w, aes(x=long, y = lat, group = group), fill = "grey80") +
    geom_raster(interpolate = TRUE) +
    scale_fill_gradientn(colours = rev(rainbow(10)), na.value = NA) +
    theme_bw() +
    coord_fixed(1.3,xlim = xpos, ylim = ypos) + 
    facet_wrap(~ date, nrow = 2)

```

### Extracting Multiple, non-consecutive time periods

What if we want to extract data from non-consecutive time periods, for example the last six Decembers for an area? <span style="color:blue">SST</span> was unusually warm off of Catalina Island in Dec 2014, and we show this by comparing the same day in December for six years, 2009-2014:

```{r getCatalina}

require(reshape2)

allyears=NULL
xpos<-c(238,243)
ypos<-c(30,35)
for (year in 2009:2014){
  tpos<-rep(paste(year,'-12-20',sep=""),2)
  tmp <- xtracto_3D(xpos,ypos,tpos,'jplMURSST')
#ggplot is sensitve that the grid is absolutely regular
  tmp$latitude <- seq(range(tmp$latitude)[1], range(tmp$latitude)[2], length.out=length(tmp$latitude))
  tmp$longitude <- seq(range(tmp$longitude)[1], range(tmp$longitude)[2], length.out=length(tmp$longitude))  
  dimnames(tmp$data) <- list(long = tmp$longitude, lat = tmp$latitude) 
  ret <- melt(tmp$data, value.name="sst")
  ret <- cbind(date=tmp$time, ret)
  allyears<-rbind(allyears,ret)
}
allyears$long <- allyears$long-360 
```

Plotting the result using faceting in ggplot2:


```{r plotCatalina,fig.width=7.5, fig.height = 4, fig.align='center'}
require(ggplot2)
require(mapdata)
xpos<-xpos-360
#long<-allyears$longitude-360
#lat<-allyears$latitude
coast <- map_data("worldHires", ylim = ypos, xlim = c(-123.5,-120.))

ggplot(data = allyears, aes(x = long, y = lat, fill = sst)) + 
  geom_polygon(data = coast, aes(x=long, y = lat, group = group), fill = "grey80") +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colours = rev(rainbow(10)), na.value = NA) +
  theme_bw() +
  coord_fixed(1.3,xlim = xpos, ylim = ypos) + 
  facet_wrap(~ date, nrow = 2)

```

This example uses the <span style="color:red">MUR (Multi-scale Ultra-high Resolution)</span> <span style="color:blue">SST</span>, a high quality ultra-high resolution <span style="color:blue">SST</span> produced by the <span style="color:red">MISST</span> project.

## Using `xtractogon` 

The `xtractogon` function extracts a time-series of satellite data that are within a user supplied polygon. As an example, the boundary points of the  Monterey Bay National Marine Sanctuary are available in the `mbnms` dataset which are loaded with the `xtractomatic` package.  The `mbnms` dataset containes two variables, <span style="color:blue">Latitude</span> and <span style="color:blue">Longitude</span>, which define the boundaries of the sanctuary:  

```{r mbnms}
require(xtractomatic)
str(mbnms)
```

The example below extracts a time-series  of monthly <span style="color:red">VIIRS</span> chlorophyll data within the MBNMS. 
 
```{r mbnmsChla}
require(xtractomatic)
tpos<-c("2014-09-01","2014-10-01")
#tpos <-as.Date(tpos)
xpos<-mbnms$Longitude
ypos<-mbnms$Latitude
sanctchl <- xtractogon(xpos,ypos,tpos,'erdVH2chlamday')
str(sanctchl)
```

The extract (see`str(sanctchl)`) contains two time periods of <span style="color:blue">chlorophyll</span> masked for data only in the sanctuayr boundaries.  A plot of the the second time period: 

```{r mbnmsChlaPlot, fig.width=5,fig.height = 5, fig.align='center'}
require(RColorBrewer)
require(ggplot2)
require(mapdata)
xlim = c(-123.5,-121.)
ylim = c(35,38)
mbnmsFrame<-mapFrame(sanctchl$longitude+360,sanctchl$latitude,log(sanctchl$data[,,2]))
my.col <- colorRampPalette(rev(brewer.pal(11, "RdBu")))(5.5)  
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(-2,4.5)) +
    ylab("latitude") + xlab("longitude") +
    coord_fixed(1.3,xlim = xlim, ylim = ylim) + 
    ggtitle(paste("log(Chla) in MBNMS",sanctchl$time[2]))
myplot
```


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. `xtractogon` can extract the bathymetry data for the MBNMS from the <span style="color:red">ETOPO</span> dataset:

```{r mbnmsBathy}
require(xtractomatic)
tpos<-c("2014-09-01","2014-10-01")
#tpos <-as.Date(tpos)
xpos<-mbnms$Longitude
ypos<-mbnms$Latitude
bathy <- xtractogon(xpos,ypos,tpos,'ETOPO180')
str(bathy)
```

Mapping the data to show the canyon:

```{r mbnmsBathyPlot, fig.width=5,fig.height = 5, fig.align='center'}
require(ggplot2)
require(RColorBrewer)
require(mapdata)
xlim = c(-123.5,-121.)
ylim = c(35,38)
mbnmsFrame<-mapFrame(bathy$longitude+360,bathy$latitude,bathy$data[,,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
```


## What happens when you request an extract

When you make an `xtractomatic` request, particularly for track data using the `xtracto` function, it is important to understand what is extracted, because the remote dataset requested likely will have a different temporal and spatial resolution then the local dataset.

Specifically, let `longitude`, `latitude` and `time` be the coordinate system of the remote <span style="color:blue">ERDDAP</span> dataset, and let `xpos`, `ypos` and `tpos` be the bounds of a request.  Then the <span style="color:blue">ERDDAP</span> request is based on the nearest grid point of the <span style="color:blue">ERDDAP</span> dataset:

```{r nearGrid,eval=FALSE}
latitude[which.min(abs(latitude - ypos[1]))]  # minimum latitude
latitude[which.min(abs(latitude - ypos[2]))]  # maximum latitude
longitude[which.min(abs(longitude- xpos[1]))] # minimum longitude
longitude[which.min(abs(longitude - xpos[2]))] # maximum longitude
isotime[which.min(abs(time- tpos[1]))] # minimum time
isotime[which.min(abs(time - tpos[2]))] # maximum time
```

where  `tpos` and `time` have been converted to an <span style="color:blue">R</span> date format so that it is a number rather than a string. For example, the FNMOC 6-hourly Ekman transports are on a 1-degree grid. A request for the data at a longitude of 220.2 and a latitude of 38.7 will return the result at a longtiude of 220 and a latitude of 39.

## Selecting a dataset {#dataset}

As seen in the examples above, the functions  `xtractomatic`, `xtracto`, `xtracto_3D`, and `xtractogon` all use the parameter <span style="color:blue">dtype</span> to specify the combination of dataset and parameter to use. The datasets accessed by `xtractomatic` are a subset of the almost 900 datasets available on ERD’s ERDDAP server at http://coastwatch.pfeg.noaa.gov/erddap/index.html. The datasets chosen to be accessible by `xtractomatic` are those that are the most useful to the majority of our users. This version of `xtractomatic` accesses 153 datasets. Information about a particular dataset is availble via the `getInfo` function. `getInfo` takes either the numbered value of a dataset, or the datatype name (the parameter <span style="color:blue">dtypename</span>), which is the first information returned by `getInfo`. `getInfo` does not return the number of the dataset, the use of the <span style="color:blue">dtypename</span>  is encouraged rather than the number associated with the dataset. It is strongly recommended that <span style="color:blue">dtype</span> not be passed as a number, as the numbers can change and will not be used in the next version.

There are a large number of temporally composited datasets available. The temporal composite time interval is usually indicated in the <span style="color:blue">dtypename</span> (for example "phssta8day" is an 8-day composite), <span style="color:blue">datasetname</span> (for example "erdPHssta8day") and <span style="color:blue">longname</span> (for example "SST, Pathfinder Ver 5.0, Day and Night, Global, Science Quality (8 Day Composite)") variables associated with a dataset. The time interval between data points in the time-series is indicated by the <span style="color:blue">timespacing</span> parameter. For example note the difference in timespacing for <span style="color:blue">dtypename</span>="mhchla8day" and <span style="color:blue">dtypename</span>="mbchla8day", which are both 8-day composites.

```{r timespacing}
getInfo('mhchla8day') 
getInfo('mbchla8day') 
```

The difference is that the latter dataset (`mbchla8day`) is a running 8-day composite, so there is a value for every day, while the former dataset (`mhchla8day`) calculates composites for non-overlapping 8-day periods.

The `searchData` function can be used to search for a subset of that satisfy certain criteria. `searchData` requires a list of lists, where each element of the first item of the list has to be one of `dtypename, datasetname, longname, varname`, and the other is the string to search for in that field. For example to find all of the datasets that have a `varname` that contains `chl` and a `datasetname` that contains `mday` (most monthly datasets have `mday` in their name):

```{r searchData }
list1<-list('varname','chl')
list2<-list('datasetname','mday')
mylist<-list(list1,list2)
searchData(mylist)
```

which returns among others the three <span style="color:blue">chlorophyll</span> datasets extracted in the `xtracto_3D` section. Note that since the function requires a "list of lists", if only doing the first search the code would be:

```{r searchData2, eval=FALSE }
mylist<-list(list('varname','chl'))
searchData(myList)
```

`searchData` iteratively refines the matches in the order given by the list of list, thus it can be viewed as doing a logical "AND" on the matches.

The <span style="color:blue">longname</span> of a dataset often contains a lot of this information. Alternatively you can generate a list of the <span style="color:blue">datasetname</span> and <span style="color:blue">longname</span> of all of the datasets accessed by `xtractomatic` (we only show the first five lines of output):

```{r erddapStruct}
# to actually use do:
# cat(paste(xtractomatic:::erddapStruct[,1],':' ,xtractomatic:::erddapStruct[,3],"\n"))
cat(paste0(xtractomatic:::erddapStruct[1:5,1],': ' ,xtractomatic:::erddapStruct[1:5,3],"\n"))
```

and then use `getInfo` with the <span style="color:blue">datasetname</span> to see all of the information on that dataset

### Topographic datasets

`xtractomatic` can extract data from two topographic datasets (<span style="color:red">ETOPO</span>)  , with <span style="color:blue">dtypename</span> of "ETOPO180" and "ETOPO360", which have longitudes ranging from -180 to 180 and 0 to 360 respectively. Currently the `getInfo` and `searchData` commands do not work with these two datasets, but `xtracto`, `xtracto_3D`, and `xtractogon` do. All of these routines require input of a time vector (tpos). The bathymetry datasets are timeless so it does not matter what you input for this parameter, however, it must be given.



## What has Changed and Steps Going Forward {#new}

In this and all past version of `xtractomatic`, information about the datasets was built into the functions.  In the present case the dataset information is stored in a dataframe called `erddapStruct` which is auotmatically loaded when you loaded the package.

This has the disadvantage that it is tied to specific datasets and a specifc <span style="color:blue">ERDDAP</span> server, rather than working with any dataset served by any <span style="color:blue">ERDDAP</span>.  

The dataset was specified by the `dtype`, which could be either a number or a character string, and references the information in that dataframe.  The `dtype` names in this version for the most part are same, **BUT THE `dtype` NUMBERS HAVE CHANGED**. If you have used `xtractomatic` before you need to check that the `dtype` information that you are passing is what you think it is.  You can do this using the supplied `getInfo` or `searchData` functions or peruse the `erddapStruct` dataframe.
There are many more datasets available in this version, adn now "last" can be used as time to get the last time period available for that dataset.

### `dtype` as number is deprecated

It is strongly recommended that you move away from passing `dtype` as a number, as the numbers can change and will not be in the next version.  Any <span style="color:blue">ERDDAP</span> dataset has an <span style="color:blue">ERDDAP</span> "dataset ID", as well one or more variables associated with that "dataset ID". Some datasets, such as the <span style="color:blue">ASCAT winds</span> have multiple variables associated with the dataset:

<http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdQAstress1day.html>.

This information is contained in `erddapStruct$datasetname` and `erddapStruct$varname`. Since a single dataset may have more than one variable,  `dtype` as character string is a way to have a unique dataset id and variable name in a single variable.  In the next version, no information about the datasets will be built into the package, and `xtractomatic` will rely on the search and query capabilities of <span style="color:blue">ERDDAP</span>. Then the request will be based on a list of the form  `c(dataset id, varname1,  varname2)`.  This is the approach being used in the <span style="color:blue">ROpenSci</span> projects `rnoaa` and `rerddap`, and we will be making use of parts their code.
