In AustralianAntarcticDivision/raadtools: Tools for Synoptic Environmental Spatial Data

## included only to make the vignette work (temporary hack)
options(default.datadir = "/rdsi/PRIVATE")

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 and https://github.com/AustralianAntarcticDivision/raadtools.

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

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

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 (fromraster::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.

library(help = raadtools)

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

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.

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.

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 theraster` 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

geb <- readtopo("gebco_08", xylim = extent(150, 220, -75, -60))

=======

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

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.

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

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)

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

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 knitrhttp://yihui.name/knitr/

##[@ncdf4; @R; @raster; @rgdal; @rgeos; @Rnews:Pebesma+Bivand:2005]

AustralianAntarcticDivision/raadtools documentation built on Feb. 14, 2024, 6:28 a.m.