stars and the R data cube for raster and vector time series

Space and time as array dimensions

The stars package is described in the README as providing classes and methods for reading, manipulating, plotting and writing spatiotemporal data cubes, to the extent that there are proper formats for doing so.

The framing is that spatiotemporal data often comes in the form of dense arrays, with space and time being array dimensions. The examples the README gives are socio-economic or demographic data, environmental variables monitored at fixed stations, raster maps, time series of satellite images with multiple spectral bands, spatial simulations, and climate or weather model output. Those are six quite different domains with one thing in common, a natural multi-dimensional structure that a data frame cannot express.

The canonical case is a cube where two dimensions are spatial raster dimensions and the third is time or band. But the README is explicit that data cubes also include higher-dimensional ones, hypercubes such as a five-dimensional cube where time, spectral band and sensor form dimensions, and lower-dimensional ones such as a plain raster image, which is a cube with the time dimension collapsed.

The last of those is a good check on whether the abstraction earns its keep. A single GeoTIFF has no time dimension, and stars still represents it as a cube with one attribute, which means the same plotting and subsetting code works whether you have one image or a thousand.

The audience is an R user with data that is genuinely a cube, and the readme's opening line, to the extent that there are proper formats for doing so, is an honest statement that this is a data model with limits rather than a universal spatial toolkit.

Raster grids that are not regular, and vector cubes

Two design decisions shape what the package can hold.

The first is that raster data do not need to be regular and aligned with North and East. The package supports regular rasters and also rotated, sheared, rectilinear and curvilinear rasters. That is a bigger claim than it sounds, because a regular raster can be indexed arithmetically while a curvilinear one cannot, so every operation that assumes a grid has to branch on the coordinate type. A package that supports all five in one object model is giving up some speed for uniformity, and the payoff is that the same code paths work regardless of how the source data happened to be gridded.

The second is the vector data cube. Where there is no pair of regularly discretised spatial dimensions, there is instead a single dimension pointing at distinct spatial feature geometries, such as polygons denoting administrative regions or points denoting sensor locations. NetCDF's CF conventions call this a discrete axis, and stars adopts the terminology.

So a cube can hold pixels or it can hold features, with the same dimensional machinery underneath. That matters for the aggregation example in the README, where a raster cube is aggregated by a vector cube of counties, taking the maximum per county. Two different kinds of spatial object meeting in one aggregation expression is a capability you do not get from separate raster and vector code paths.

stars_proxy: a cube that stays on disk

The package provides stars_proxy objects, currently only when read through GDAL, which contain only the dimensions metadata and pointers to the files on disk.

The behaviour is lazy. Reading and processing data is postponed to the moment that pixels are really needed, at plot time or when writing to disk, and it is done at the lowest spatial resolution possible that still fulfills the resolution of the graphics device. That last clause is the interesting one. A proxy does not simply defer the read, it also decides how much to read, and the deciding factor is the output device, not the data. Plotting to a small window reads less than plotting to a high resolution image, and neither reads the whole cube unless it has to.

The method list in the README shows how much of the ordinary stars vocabulary survives in the proxy world, and it is a long one. Beyond the indexing operators there is aggregate, aperm, as.data.frame, coerce, dim, droplevels, filter, hist, image, initialize, is.na, merge, mutate, plot, prcomp, predict, print, pull, rename, select, show, slice, split, st_apply, st_as_sf, st_as_stars, st_crop, st_dimensions, st_downsample, st_mosaic, st_normalize, st_redimension, st_sample, st_set_bbox, transmute and write_stars.

Two of those names are worth singling out. st_downsample and st_mosaic are the operations a proxy makes possible that a fully materialised array would make awkward, because reducing the read before it happens is far cheaper than reducing it after. And st_set_bbox is what lets a lazy cube be clipped to a region of interest without ever reading the rest of the file.

Two read paths, and GDAL doing the work

The package provides two functions to read data, read_ncdf and read_stars, where the latter reads through GDAL. The README notes that in the future both will be integrated in read_stars.

For NetCDF files, package RNetCDF is used. For reading through GDAL, package sf provides the binary linking to GDAL. So the two paths are not two ways to do the same thing; they are two different backends with different capabilities, and the choice determines which formats you can open and whether you can get a stars_proxy.

The same delegation applies to the operations. For vector and raster operations, stars uses as much as possible the routines available in GDAL and PROJ, with st_transform, rasterize, polygonize and warp named explicitly. The documentation points to the vignette on vector-raster conversions, reprojection and warping for the details.

This is the right call for an R package, because reimplementing reprojection is a bad idea and PROJ is the reference implementation. It also has a consequence for adoption, which is that the GDAL and PROJ versions your system has determine what works. A package whose correctness depends on the GDAL build underneath it is a package with an implicit dependency on the system administrator's choices.

The example that shows the object model in practice

The README walks through a hurricane precipitation example, and it is the clearest demonstration of what a stars object holds.

The first block reads a raster and slices a band:

library(dplyr)
library(stars)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
read_stars(tif) |>
  slice(index = 1, along = "band") |>
  plot()

The second imports a curvilinear grid with hourly precipitation values of a hurricane and plots the first 12 time steps, using a file from the package and reading it through GDAL subdatasets:

prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars")
(prec = read_stars(gdal_subdatasets(prec_file)[[1]]))

The printed structure of that object is the specification: three dimensions and one attribute, with x and y carrying no offset and no delta because the grid is curvilinear, a WGS 84 reference system, a time dimension running from 1 to 23 in one hour steps from 2018-09-13, and a note that the values are x and y. The x and y dimensions are 87 by 118 with longitude and latitude ranges, and the time dimension is NULL, which is the marker of a curvilinear grid.

The README then intersects the result with the counties of North Carolina, transforming the counties from NAD27 to WGS 84 to match, and taking the maximum precipitation intensity per county. That last step is the payoff of the vector cube support, because aggregating a raster by a polygon layer with a single function call is the operation that most workflows need and that separate raster and vector code makes awkward.

Editorial conclusion

stars is a good fit for an R analyst working with raster time series, satellite imagery or gridded climate output who wants a single object model for a cube rather than a matrix with a time attribute bolted on, because the package handles regular, rotated, sheared, rectilinear and curvilinear grids uniformly and keeps vector cubes in the same framework. It is a poor fit if your data is small enough for a data frame, since the overhead of a dimensional object buys you nothing there, and a poor fit if you need the full CF conventions toolbox, since the package documents only the discrete axis concept. Before adopting, check that your data can be read through one of the two documented paths, read_ncdf or read_stars via GDAL, and remember that stars_proxy objects are currently only created when reading through GDAL, so a NetCDF workflow does not get the lazy path.

Frequently asked questions

Official sources

1. Issues
2. License: Apache-2.0
3. Project website
4. README
5. r-spatial/stars on GitHub