R/get_sat_data.R
get_sat_data.RdReturned .tif values are 16 bit signed integers. For reflectance these have values 1 to 10000 with no data as -999 (as per the original Digital Earth Australia cube). For indices these have values of index * 10000 with no data as -32768.
get_sat_data(
x,
start_date = "2013-03-01",
end_date = paste0(as.numeric(format(Sys.Date(), "%Y")) - 1, "-12-31"),
out_dir = NULL,
source_url = "https://explorer.sandbox.dea.ga.gov.au/stac",
collections = c("ga_ls9c_ard_3", "ga_ls8c_ard_3"),
period = "P1M",
property_filter = NULL,
aggregation_func = "median",
resampling_method = "bilinear",
layers = c("green", "blue", "red"),
excludes = "nbar_",
indices = list(gdvi = c("green", "nir"), ndvi = c("nir", "red")),
mask = list(band = "oa_fmask", mask = c(2, 3)),
save_cube = FALSE,
cores = 1,
force_new = FALSE,
attempts = 1,
sleep = 5,
gdalcubes_config = list(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "5%", VSI_CACHE_SIZE =
"100000000", GDAL_HTTP_MULTIPLEX = "YES", GDAL_INGESTED_BYTES_AT_OPEN = "32000",
GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR", GDAL_HTTP_VERSION = "2",
GDAL_HTTP_MERGE_CONSECUTIVE_RANGES = "YES", GDAL_NUM_THREADS = cores),
...
)SpatRaster, or path to raster (needs to be path if running in parallel), to use as extent, resolution and crs for output.
Character or Date. e.g. "2013-03-01". These specify the temporal extent of the cube.
Character. Path into which results will be saved.
Character. URL specifying the STAC endpoint.
Character. Name of collection(s) within the STAC endpoint from which to build the cube. e.g. DEA Surface Reflectance (Sentinel-2, Collection 3) could be specified by c("ga_s2am_ard_3", "ga_s2bm_ard_3").
Character specifying the temporal cell size in ISO 8601. e.g. "P1M" specifies a time period of one month.
Passed to the property_filter argument of
gdalcubes::stac_image_collection().
Passed to the aggregation argument of
gdalcubes::cube_view().
Passed to the resampling argument of
gdalcubes::cube_view().
Character. Regular expressions, to search within the items returned by the rstac query, defining the layers (bands) from which to build the cube.
Character. Regular expressions, to search AND EXCLUDE within the items returned by the rstac query. e.g. the default 'nbar_" excludes, say, nbar_blue, but leaves nbart_blue.
Named list of pairs of layers to use to build indices within
the cube. e.g. indices = list(ndvi = c("nir", "red")) will create a layer in
the cube 'ndvi' as (nir - red) / (nir + red). Currently only relevant if
save_cube is TRUE.
Named list specifying any band, and the levels within that band, that identify pixels in an image that should be excluded from the cube. Default mask identifies cloud and cloud shadow within default source_url and collections.
Logical. If TRUE the cube will be saved as an individual
.tif file (in out_dir) per band and indice. The name of each .tif will be
layer__start_date.tif.
Numeric. Number of cores to pass to gdalcubes configuration for "GDAL_NUM_THREADS".
Logical. If a band already exists in out_dir, recreate it?
Numeric. Occasionally, rstac::get_request() or
rstac::items_fetch() returns no items (probably due to issues with
source_url). If this is the case, retry up to attempts times to return
the appropriate items.
Numeric. Time in seconds to wait between attempts.
List in the form of key = value pairs suitable for
gdalcubes::gdalcubes_set_gdal_config(). Also see
https://gdalcubes.github.io/source/concepts/config.html. These configuration
settings attempt to "improve computation times and I/O or network transfer
performance".
Passed to gdalcubes::write_tif(). Arguments, x, dir,
prefix and pack are already passed though, so attempting to use those
via dots will fail.
If save_cube = FALSE, a data cube proxy object. If save_cube = TRUE, invisible(NULL) and .tif file in out_dir for every time period, and layer or indice, requested
#library(envRaster)
out_dir <- file.path(system.file(package = "envRaster"), "examples")
cube_dir <- fs::path(out_dir, "cube")
settings_path <- fs::path(cube_dir, "settings.rds")
if(!file.exists(settings_path)) {
source(fs::path(out_dir, "name_env_tif_ex.R"))
}
settings <- rio::import(settings_path)
#> Warning: Missing `trust` will be set to FALSE by default for RDS in 2.0.0.
base <- make_base_grid(aoi # aoi is sf provided with envRaster
, out_res = 30
, out_epsg = 8059
, out_file = fs::path(cube_dir, "base.tif")
, datatype = "INT1U"
, overwrite = TRUE
)
# Can't get indices with save_cube = FALSE but can pass cube onto apply_pixel
cube <- get_sat_data(x = base
, start_date = "2023-01-01"
, end_date = "2023-12-31"
, collections = c("ga_ls9c_ard_3"
, "ga_ls8c_ard_3"
)
, property_filter = function(x) {x[["eo:cloud_cover"]] < 20}
, period = "P3M"
, layers = c("nir", "red")
, indices = NULL
, save_cube = FALSE
)
#>
|
|======================================================================| 100%
#> attempt: 1. features: 91
if(FALSE) {
# takes a minute or two
cube %>%
gdalcubes::apply_pixel("(nbart_nir - nbart_red) / (nbart_nir + nbart_red)", "ndvi") %>%
gdalcubes::animate(col = viridis::viridis
, zlim = c(0, 1)
, downsample = FALSE
)
}
# Save a cube
if(FALSE) {
# takes a minute or two
get_sat_data(x = base
, start_date = "2022-12-01"
, end_date = "2023-11-30"
, collections = settings$collection
, period = settings$period
, layers = NULL
, indices = list("ndvi" = c("nir", "red"))
, save_cube = TRUE
, out_dir = settings$out_dir
# passed to gdalcubes::write_tif
, creation_options = list("COMPRESS" = "NONE")
)
# biggest contrast in ndvi in autumn. least contrast in spring
ndvi_tifs <- fs::dir_ls(settings$out_dir, regexp = "tif$")
ndvis <- terra::rast(ndvi_tifs)
names(ndvis) <- gsub("ndvi__|\\.tif", "", basename(ndvi_tifs))
terra::global(ndvis
, fun = \(x) quantile(x, probs = 0.75) - quantile(x, probs = 0.25)
)
ndvis %>%
terra::panel()
}