Repository files navigation

Spatiotemporal data often comes in the form of dense arrays, with spaceand time being array dimensions. Examples include

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

This R package provides classes and methods for reading, manipulating,plotting and writing such data cubes, to the extent that there areproper formats for doing so.

The canonical data cube most of us have in mind is that where twodimensions represent spatial raster dimensions, and the third time (orband), as e.g. shown here:

By data cubes however we also consider higher-dimensional cubes(hypercubes) such as a five-dimensional cube where in addition to time,spectral band and sensor form dimensions:

or lower-dimensional cubes such as a raster image:

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

Raster data do not need to be regular and aligned with North/East, andpackagestars supports besidesregular alsorotated,sheared,rectilinear andcurvilinear rasters:

Vector data cubes arise when we do not have two regularly discretizedspatial dimensions, but a single dimension that points to distinctspatial feature geometries, such as polygons (e.g. denotingadministrative regions):

or points (e.g. denoting sensor locations):

NetCDF’s CF-convention calls this adiscreteaxis.

stars provides two functions to read data:read_ncdf andread_stars, where the latter reads through GDAL. (In the future, bothwill be integrated inread_stars.) For reading NetCDF files, packageRNetCDF is used, for reading through GDAL, packagesf provides thebinary linking to GDAL.

For vector and raster operations,stars uses as much as possible theroutines available in GDAL and PROJ (e.g. st_transform,rasterize,polygonize,warp). Read more about this in the vignette onvector-raster conversions, reprojection,warping.

Packagestars providesstars_proxy objects (currently only when readthrough GDAL), which contain only the dimensions metadata and pointersto the files on disk. These objects work lazily: reading and processingdata is postponed to the moment that pixels are really needed (at plottime, or when writing to disk), and is done at the lowest spatialresolution possible that still fulfills the resolution of the graphicsdevice. More details are found in thestars proxyvignette.

The following methods are currently available forstars_proxy objects:

methods(class="stars_proxy")#  [1] [               [[<-            [<-             adrop#  [5] aggregate       aperm           as.data.frame   c#  [9] coerce          dim             droplevels      filter# [13] hist            image           initialize      is.na# [17] Math            merge           mutate          Ops# [21] plot            prcomp          predict         print# [25] pull            rename          select          show# [29] slice           slotsFromS3     split           st_apply# [33] st_as_sf        st_as_stars     st_crop         st_dimensions<-# [37] st_downsample   st_mosaic       st_normalize    st_redimension# [41] st_sample       st_set_bbox     transmute       write_stars# see '?methods' for accessing help and source code

In the following, a curvilinear grid with hourly precipitation values ofa hurricane is imported and the first 12 time steps are plotted:

prec_file= system.file("nc/test_stageiv_xyt.nc",package="stars")(prec= read_stars(gdal_subdatasets(prec_file)[[1]]))# stars object with 3 dimensions and 1 attribute# attribute(s):#                                         Min. 1st Qu. Median     Mean 3rd Qu.# Total_precipitation_surface... [kg/m^2]    0       0   0.75 4.143009    4.63#                                           Max.# Total_precipitation_surface... [kg/m^2] 163.75# dimension(s):#      from  to                  offset   delta  refsys# x       1  87                      NA      NA  WGS 84# y       1 118                      NA      NA  WGS 84# time    1  23 2018-09-13 19:00:00 UTC 1 hours POSIXct#                                  values x/y# x    [87x118] -80.61 [°],...,-74.88 [°] [x]# y      [87x118] 32.44 [°],...,37.62 [°] [y]# time                               NULL# curvilinear grid# or: (prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE))sf::read_sf(system.file("gpkg/nc.gpkg",package="sf"),"nc.gpkg")|>   st_transform(st_crs(prec))->nc# transform from NAD27 to WGS84nc_outline= st_union(st_geometry(nc))plot_hook=function() plot(nc_outline,border='red',add=TRUE)prec|>  slice(index=1:12,along="time")|>  plot(downsample= c(3,3,1),hook=plot_hook)

and next, intersected with with the counties of North Carolina, wherethe maximum precipitation intensity was obtained per county, andplotted:

a= aggregate(prec,by=nc,FUN=max)plot(a,max.plot=23,border='grey',lwd=.5)

We can integrate over (reduce) time, for instance to find outwhen themaximum precipitation occurred. The following code finds the time index,and then the corresponding time value:

index_max=function(x) ifelse(all(is.na(x)),NA, which.max(x))b= st_apply(a,"geom",index_max)b|> mutate(when= st_get_dimension_values(a,"time")[b$index_max])|>  select(when)|>  plot(key.pos=1,main="time of maximum precipitation")

With packagecubble, we can make a glyph map to see the magnitude andtimings of county maximum precipitation:

library(cubble)library(ggplot2)a|> setNames("precip")|>  st_set_dimensions(2,name="tm")|>units::drop_units()|>  as_cubble(key=id,index=tm)->a.cba.cb|>  face_temporal()|>  unfold(long,lat)|>  mutate(tm= as.numeric(tm))|>  ggplot(aes(x_major=long,x_minor=tm,y_major=lat,y_minor=precip))+  geom_sf(data=nc,inherit.aes=FALSE)+  geom_glyph_box(width=0.3,height=0.1)+  geom_glyph(width=0.3,height=0.1)

Packagegdalcubes can be used to create data cubes (or functions fromthem) from image collections, sets of multi-band images with varying

• spatial resolution
• spatial extent
• coordinate reference systems (e.g., spread over multiple UTM zones)
• observation times

and does this by resampling and/or aggregating over space and/or time.It reuses GDAL VRT’s and gdalwarp for spatial resampling and/or warping,and handles temporal resampling or aggregation itself.

ncdfgeom reads and writes vector data cubes from and to netcdf filesin a standards-compliant way.

Packagesraster and its successor,terra are powerful packages forhandling raster maps and stacks of raster maps both in memory and ondisk, but do not address

• non-raster time series,
• multi-attribute rasters time series
• rasters with mixed type attributes (e.g., numeric, logical, factor,POSIXct)
• rectilinear or curvilinear rasters

A list ofstars commands matching existingraster commands is foundin thiswiki.A list of translations in the opposite direction (fromstars toraster orterra) still needs to be made.

A comment on the differences betweenstars andterra is foundhere.

• blog posts:first,second,third, andnewerblog posts
• vignettes
• the originalR Consortiumproposal.

This project has been realized with financialsupportfrom the