The {stars} package is awesome for working with data cubes and higher-dimensional raster GIS datasets. Here's a pretend example I've created called landsat7_sweep with the following dimensions:
- x, y are the spatial dimensions
- band is the wavelength bands collected by landsat7
- time is the time dimension
library("stars")
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
library("tidyverse")
landsat7_data <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
landsat7_data_2015_01_01 <- landsat7_data
landsat7_data_2015_01_02 <- landsat7_data
landsat7_sweep <- c(landsat7_data_2015_01_01, landsat7_data_2015_01_02) %>%
st_redimension() %>%
st_set_dimensions(4, values = c(as.Date("2015-01-01"), as.Date("2015-01-02")), name = "time")
landsat7_sweep
#> stars object with 4 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#> L7_ETMs.tif.L7_ETMs.tif.1
#> Min. : 47.00
#> 1st Qu.: 65.00
#> Median : 76.00
#> Mean : 77.34
#> 3rd Qu.: 87.00
#> Max. :255.00
#> dimension(s):
#> from to offset delta refsys point values x/y
#> x 1 349 288776 28.5 UTM Zone 25, Southern Hem... FALSE NULL [x]
#> y 1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE NULL [y]
#> band 1 6 NA NA NA NA NULL
#> time 1 2 2015-01-01 1 days Date NA NULL
Real-world examples of these datasets require much more data preparation:
- downloading many files
- iteratively reading in .tif files
- converting attributes to dimensions
I would love a friendly way to export landsat7_sweep for others to use. Unfortunately, write_stars() isn't suitable for two reasons:
-
As noted in the stars website, write_stars() can only export single attribute stars objects
-
write_stars() fails to export objects with more than 3 dimensions
landsat7_sweep %>%
write_stars("landsat7_sweep")
#> Error in sf::gdal_write(obj, ..., file = dsn, driver = driver, options = options, : dimensions don't match
What advice would you give for providing landsat7_sweep as a dataset for others?
Is the best bet to use .RData files created as follows:
save(landsat7_sweep, file = "landsat7_sweep.RData")