raadtools.Rmd

---
title: "Using the raadtools package"
author: "Michael D. Sumner"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
        fig_width: 9
        fig_height: 9
vignette: >
  %\VignetteIndexEntry{raadtools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: raadtools.bib
---


```{r,include=FALSE}
## included only to make the vignette work (temporary hack)
options(default.datadir = "/rdsi/PRIVATE")
```
```{r,setup, include=FALSE}
library(knitr)
# set global chunk options
opts_chunk$set(fig.path='figure/minimal-', fig.align='center', fig.show='hold', cache = TRUE)
options(replace.assign=TRUE,width=90)
```

# Overview 

The `raadtools` package provides a set of helper functions based on functionality in the `raster` package, to simplify fast and efficient use of time series environmental data including



* chlorophyll-a concentation
* colour images
* sea ice concentration
* sea surface current vectors
* sea surface height and height anomaly
* sea surface temperature 
* sea surface wind vectors. 

There is ongoing work to simplify extension of these tools to a wider range of data, but here we focus on the existing examples just to highlight the preferred ways of working with these data in `R`. 

# Introduction

This document introduces the `raadtools` package. Some
familiarity with `R` and the `raster` package is
assumed.

This package relies on the availability of a local file repository which must be present, see the componente packages [https://github.com/AustralianAntarcticDataCentre/raadsync](raadsync) and [https://github.com/AustralianAntarcticDivision/raadtools](raadtools). 

The very first step is to load the package, which also triggers an attempt to connect to the file repository. 

```{r,eval=TRUE}
library(raadtools)
library(raster)
```



(TO FIX: In subsequent code sections I'm also attaching raadfiles explicitly to make this document work, but it's not necessary generally. )


```{r,eval=TRUE}
library(raadfiles)
library(raadtools)

ice <- readice()
plot(ice)
```


# Simple data read

There are functions to read data based on a sequence of date-time
values.


At the simplest level, these functions may be run without setting any
arguments and will return the first available time step for the
default data source. This is handy for becoming familiar with the
kinds of data before getting really specific. Plotting the resulting
`ice} object provides a familiar polar view. These data are a
single time slice from the daily passive microwave sea ice concentration, here the very first one that is available (from `raster::getZ(ice)`).

Other read functions include `reachla`, `readcurr`,
`readsst`, `readwind` and `readtopo`. Each of
these will read a default data set in a similar way to the ice
example, see their respective help pages and explore the index page of
`raadtools` to see an overview of functionality.

Load up the package help index for an overview of all the help topics
and functions.

```{r,eval=FALSE,echo=TRUE}
library(help = raadtools)
```

Read a simple example from each of the data types mentioned above and plot.


```{r,echo=TRUE,output=FALSE}
library(raadtools)
library(raadfiles)
## chlorophyll-a is currently only patchily supported
##chla <- readchla("2009-10-01")
curr <- readcurr("2012-01-09", magonly = TRUE, lon180 = FALSE)
sst <- readsst(c("2000-01-01", "2009-06-08"))
wind <- readwind()
##topo <- readtopo("george_v_terre_adelie", xylim = extent(140, 141.5, -65, -64))

apal <- chl.pal(palette = TRUE)
op <- par(mfcol = c(2, 2))
#plot(chla, col = chl.pal(apal$breaks), breaks = apal$breaks,
#     main = "Johnson chl-a", legend = FALSE)
plot(curr, main = "ocean surface current magnitude")
plot(sst[[2]], main = "SST")
#plot(topo, main = "George V Terre Adelie bathymetry",
#     col = topo.colors(68)[-c(40:68)])
par(op)

plot(wind, nr = 2)
```

The read calls above apply slightly different options in each case,
with no arguments at all, with a vector of date-times, `xylim`
to restrict the spatial extent an `lon180` to control the
Atlantic vs. Pacific view for global data extending across the
dateline.



The `ice` object  is a `class(ice)` from the `raster`
package.


```{r,eval=TRUE}
ice
```

This object is a `RasterLayer` which consists of just a single
grid layer for the day listed at the beginning of the NSIDC data
set. The dimensions of the grid are listed, (in the order Y then X)
`dim(ice)[1:2]` and these refer to the number of cells seen
vertically on the page (`nrow(ice)`) and those seen horizontally
(`ncol(ice)`), and then the total number of cells (X * Y).


## Leveraging raster in R 

The `raster` package provides excellent facilities for working
with raster data, including images, surfaces, DEMs and time series
satellite data. From the `raster` package DESCRIPTION:

```
Reading, writing, manipulating, analyzing and modeling of
gridded spatial data. The package implements basic and high-level
functions and processing of very large files is supported.
```

This package should be consulted for documentation on the use of
Raster objects, and in particular the vignette `Raster` is a
good place to start. Use the `library` function's `help`
argument to get an overview of the available functions.

```{r,eval=FALSE}
library(raster)

## opens an index page to all of the package's functions
library(help = raster)
## opens a PDF document with extended help and examples
vignette("Raster")
```
<<<<<<< HEAD:raadtools.Rmd
=======

<<<<<<< HEAD
>>>>>>> 1fb87b2a2d061ccadc99c9905314e18e58286d86:vignettes/raadtools.Rmd
=======
>>>>>>> e63b5cf976c2c0469edc19520ea53405acfbacd0:vignettes/raadtools.Rmd
>>>>>>> master
See the Section sec:assumedsupport for more detail on the design
choices made with `raster` and the other package dependencies
used by `raadtools`. 


# Read functions in raadtools

This section discusses the details of individual read functions and
some of the rationale behind each one. Different data sets are stored
with different native coordinate systems, different pixel and temporal
resolutions, different data units, and cover different spatial and temporal ranges. As well
as all these options, the way data are stored varies tremendously
depending on the provider and the data type. These read functions all
return either a `Raster} object as defined by the
`raster` package with as complete and standard metadata as
possible.


Some data sets are inherently multi-layered for a single time step,
such as wind and current vectors and so these both have optional
arguments `dironly` and `magonly` to convert the
internal `U` and `V` vector components to directions or
magnitudes only. The read operation may be combined with a spatial
crop by using the argument `xylim` and there are other options
specific for some data sets such as `lon180` to control the
default orientation (these topics will be discussed in more detail in
later sections, otherwise see the documentation for each function
argument).

## Read topographic or bathymetric data


The `readtopo` (`readbathy` is a synonym) will return
one of several available data sets. See the help
page for the details and source of each one. As these data are not
dynamic on short time scales they do not provide a `date`
argument. The `ibcso` data set is optionally available on a
polar or longitude/latitude grid, controlled by the `polar`
argument. The Smith and Sandwell data set is returned in Pacific view
on the Mercator projection by default. Use `lon180 = FALSE` to
return Atlantic view.

<<<<<<< HEAD:raadtools.Rmd
```{r}
geb <- readtopo("gebco_08", xylim = extent(150, 220, -75, -60))
```

=======
```{r,eval=TRUE}
library(raadtools)
library(raadfiles)
geb <- readtopo("gebco_08", xylim = extent(150, 220, -75, -60))
plot(geb, col = viridis::inferno(100))
```
<<<<<<< HEAD
>>>>>>> 1fb87b2a2d061ccadc99c9905314e18e58286d86:vignettes/raadtools.Rmd
=======
>>>>>>> e63b5cf976c2c0469edc19520ea53405acfbacd0:vignettes/raadtools.Rmd
>>>>>>> master



# Plot vector data

```{r,eval=TRUE}
library(raadtools)
library(raadfiles)
## set up data in advance, to feed the plotting function
tst <- readcurr()
xylim <- extent(projectExtent(raster(extent(130, 150, -50, -30), crs = "+proj=longlat"), projection(tst)))

plotcurrents <- function(date = as.Date("1999-11-24"), ext = NULL, scale = 2, ...) {
  x <- readcurr(date)
  if (!is.null(ext)) x <- crop(x, ext)
  crds <- coordinates(x)
  plot(sqrt(x[[1L]]^2 + x[[2L]]^2), ...)
  x1 <- crds[,1]
  y1 <- crds[,2]
  x2 <- crds[,1] + values(x[[1L]]) * scale
  y2 <- crds[,2] + values(x[[2L]]) * scale
  op <- options(warn = -1)
  sink("/dev/null")
  arrows(x1, y1, x2, y2, length = 0.03, col = "#00000080")
  sink(NULL)
  options(op)
  invisible(NULL)
}

plotcurrents(date = as.Date("2000-01-05"), ext = xylim, scale = .1, asp = NA)

```

Animate a series of current plots as above.

```{r,eval=TRUE}
library(raadtools)
library(raadfiles)
dts <- seq(as.Date("2005-10-01"), length = 104, by = "1 week")
pal <- sst.pal(palette = TRUE)
cols <- pal$cols[seq(1, length(pal$cols), length = 16)]
cols <- gsub("ff$", "99", cols)
brks <- pal$breaks[seq(1, length(pal$breaks), length = 16)]

#library(animation)
#ani.start()
#for (i in seq_along(dts)) plotcurrents(dts[i], ext = xylim, scale = 2500, col = cols, breaks = brks)
#ani.stop()

```



# Extract spatio-temporal data from read functions

These have all be tested at some point, though maybe not in the current version of this document. 

More needs to be said here, and pictures . . .

```{r,eval=FALSE}

library(raadtools)


a <- structure(list(x = c(174, 168, 156, 111, 99, 64, 52, 46, -4,
                          -15, -30, -38, -47, -62, -87, -127, -145, -160, -161), y = c(-72,
                                                                                       -39, -50, -58, -35, -38, -48, -60, -48, -35, -37, -51, -68, -72,
                                                                                       -69, -54, -40, -49, -54)), .Names = c("x", "y"), row.names = c(NA,
                                                                                                                                                      -19L), class = "data.frame")

a$time <-  sort(as.Date("2005-01-01") + sample(c(0, 0, 0, 8, 20, 50), nrow(a), replace = TRUE))

extract(readssh, a)
extract(readssh, a, ssha = TRUE)
extract(readcurr, a, magonly = TRUE)
extract(readice, a)

#extract(readchla, time.resolution = "weekly")
#extract(readchla, a, time.resolution = "weekly")



##extract(readwind, a)
extract(readwind)
readwind(dironly = TRUE)



a$time <- sort(as.Date("1985-01-01") + sample(c(0, 0, 0, 10, 100, 50), nrow(a), replace = TRUE))


##x <- extract(readsst)
##extract(readsst, a)


##extract(readsst)
##extract(readsst, time.resolution = "daily")


##extract(readsst, a)


```

```{r,eval=FALSE,echo=FALSE}
data(aurora)

aurora$sst <- extract(readsst, aurora)
#aurora$chla <- extract(readchla, aurora)

## nothing for 2013 yet
##aurora$windmag <- extract(readwind, aurora[,1:3], magonly = TRUE)
##aurora$winddir <- extract(readwind, aurora[,1:3], dironly = TRUE)

##aurora$currmag <- extract(readcurr, aurora[,1:3], magonly = TRUE)
##aurora$currdir<- extract(readcurr, aurora[,1:3], dironly = TRUE)

##aurora$ssh <- extract(readssh, aurora[,1:3])
##aurora$ssha <- extract(readssh, aurora[,1:3], ssha = TRUE)

##aurora$ice <- extract(readice, aurora[,1:3])


## note this is slightly different, since there's no point in generalizing out time
##aurora$bathy <- extract(readbathy(), aurora[,1:2])
```

## Extract spatio-temporal data with spatial, temporal interpolation, and spatial aggregation

```{r,eval=FALSE}
d1 <- data.frame(lon = c(130, 125, 120, -52.27, 123.21, -25.42, 36.96, -36.28),
                 lat = c(-50, -52, -55,-66.31, -74.55, -69.66, -71.85, -68.56),
                 gmt = as.POSIXct("2006-02-08 00:00:00") + 24 *3600 * c(0, 1, 2, 3, 106, 118, 118, 190))

llproj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

library(raadtools)
library(rgdal)

## test extract vs extractxyt
sst1 <- extractxyt("oisst", d1)
sst2 <- extract(readsst, d1)

mag1 <- extractxyt("aviso", d1, magonly = TRUE)
mag2 <- extract(readcurr, d1, magonly = TRUE)

ice1 <- extractxyt("nsidc", d1)
ice2 <- extract(readice, d1)

## these diffs should be  zorro
max(abs(sst1 - sst2), na.rm = TRUE)
max(abs(mag1 - mag2), na.rm = TRUE)
max(abs(ice1 - ice2), na.rm = TRUE)


## now the goodstuff
sst3 <- extract(readsst, d1, method = "bilinear")
sst4 <- extract(readsst, d1, ctstime = TRUE)
sst5 <- extract(readsst, d1, ctstime = TRUE, method = "bilinear")

data.frame(sst = sst2, sst_xyinterp = sst3, sst_timeinterp = sst4, sst_spacetimeinterp = sst5)

## some time when we get chl-a
d2 <- data.frame(x = c(-40, -25, -6, 12, 33, 60, 89, 107),
                 y = c(-56, -55, -59, -62, -64, -62, -60, -60),
                 time = as.POSIXct("2006-12-18 00:00:00") + 24 *3600 * 7 * seq_len(8))

## resize by a factor in space so we get "more coverage"
# chla1 <- extract(readchla, d2)
# chla2 <- extract(readchla, d2, ctstime = TRUE, method = "bilinear")
# chla3 <- extract(readchla, d2, ctstime = TRUE, method = "bilinear", fact = 16)
# data.frame(chla = chla1, chla_spacetimeinterp = chla2, chla_resizespacetimeinterp = chla3)
# 
# chl <- readchla(d2$time, xylim = extent(as.matrix(d2[,1:2])))
# chl2 <- aggregate(chl, fact = 16)
# pal <- chl.pal(palette = TRUE)
# 
# ## compare the two plots, with and without spatial aggregation for "coverage"
# plot(chl, col = pal$col, breaks = pal$breaks, legend = FALSE, addfun = function() lines(d2[,1:2]))
# 
# plot(chl2, col = pal$col, breaks = pal$breaks, legend = FALSE, addfun = function() lines(d2[,1:2]))
# 

```

# Data sets and the native projection

All data sets availalable in `raadtools` are provided in their \emph{native projection`. This means that a data set always has a \emph{coordinate reference system` explicitly and that the data must be used with this coordinate system. A good example is the passive microwave sea ice concentration data, which are provided on a polar stereographic grid for both poles. Another examples is the Smith and Sandwell bathymetry which is provided in a Mercator grid.   Many data sets are provided in `longitude/latitude` but note that this is still not enough information since angular longitude/latitude can only be defined with reference to a particular \emph{datum}. 

R provides all of the tools required to work with these data. 

Coordinate reprojection: raw coordinates and point, line, or polygon objects can be transformed between different map projections using the 'rgdal' package. 

Flexible extraction methods: the `raster` package provides extraction methods to do a huge range of `query-based extractions` from gridded data. These include: 


* extract values from grids by points
* extract values from grids by objects: SpatialPoints, SpatialLines, SpatialPolygons
* extraction with function-based aggregation (for queries where more than one pixel is returned for an object, i.e. lines and polygons)
* extraction by objects with different map projections to the grid data, where the extraction will apply reprojection to the grid automatically
* extraction for multi-dimensional grids, i.e. extract values for multiple times for each point/line/polygon
* interpolation relative to a points location in a cell
* extraction return weights for each pixel, allowing weighting of an object's coverage relative to a cell
* extraction more abstractly, returning the cell numbers rather than values allowing flexible repetitive-extraction


Together all of these features allow practically any required extraction task. 


# Acknowledgements

The `raadtools` package was built with R

* R http://www.r-project.org
* R tools http://cran.r-project.org/bin/windows/Rtools/
* RStudio http://www.r-studio.com}
* Git http://git-scm.com}
* This  document was built with `knitr`http://yihui.name/knitr/

```{r,eval=FALSE,include=FALSE}
##[@ncdf4; @R; @raster; @rgdal; @rgeos; @Rnews:Pebesma+Bivand:2005]
```