R: Retrieve raster data from STAC endpoints

get_stac_data {rsi}

R Documentation

Retrieve raster data from STAC endpoints

Description

These functions retrieve raster data from STAC endpoints and optionally create composite data sets from multiple files. get_stac_data() is a generic function which should be able to download raster data from a variety of data sources, while the other helper functions have useful defaults for downloading common data sets from standard STAC sources.

Usage

get_stac_data(
  aoi,
  start_date,
  end_date,
  pixel_x_size = NULL,
  pixel_y_size = NULL,
  asset_names,
  stac_source,
  collection,
  ...,
  query_function = rsi_query_api,
  download_function = rsi_download_rasters,
  sign_function = NULL,
  rescale_bands = TRUE,
  item_filter_function = NULL,
  mask_band = NULL,
  mask_function = NULL,
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = c("merge", "median", "mean", "sum", "min", "max"),
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS",
    GDAL_HTTP_USERAGENT = "rsi (https://permian-global-research.github.io/rsi/)")
)

get_sentinel1_imagery(
  aoi,
  start_date,
  end_date,
  ...,
  pixel_x_size = 10,
  pixel_y_size = 10,
  asset_names = rsi::sentinel1_band_mapping$planetary_computer_v1,
  stac_source = attr(asset_names, "stac_source"),
  collection = attr(asset_names, "collection_name"),
  query_function = attr(asset_names, "query_function"),
  download_function = attr(asset_names, "download_function"),
  sign_function = attr(asset_names, "sign_function"),
  rescale_bands = FALSE,
  item_filter_function = NULL,
  mask_band = NULL,
  mask_function = NULL,
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "median",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

get_sentinel2_imagery(
  aoi,
  start_date,
  end_date,
  ...,
  pixel_x_size = 10,
  pixel_y_size = 10,
  asset_names = rsi::sentinel2_band_mapping$planetary_computer_v1,
  stac_source = attr(asset_names, "stac_source"),
  collection = attr(asset_names, "collection_name"),
  query_function = attr(asset_names, "query_function"),
  download_function = attr(asset_names, "download_function"),
  sign_function = attr(asset_names, "sign_function"),
  rescale_bands = FALSE,
  item_filter_function = NULL,
  mask_band = attr(asset_names, "mask_band"),
  mask_function = attr(asset_names, "mask_function"),
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "median",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

get_landsat_imagery(
  aoi,
  start_date,
  end_date,
  ...,
  platforms = c("landsat-9", "landsat-8"),
  pixel_x_size = 30,
  pixel_y_size = 30,
  asset_names = rsi::landsat_band_mapping$planetary_computer_v1,
  stac_source = attr(asset_names, "stac_source"),
  collection = attr(asset_names, "collection_name"),
  query_function = attr(asset_names, "query_function"),
  download_function = attr(asset_names, "download_function"),
  sign_function = attr(asset_names, "sign_function"),
  rescale_bands = TRUE,
  item_filter_function = landsat_platform_filter,
  mask_band = attr(asset_names, "mask_band"),
  mask_function = attr(asset_names, "mask_function"),
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "median",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

get_naip_imagery(
  aoi,
  start_date,
  end_date,
  ...,
  pixel_x_size = 1,
  pixel_y_size = 1,
  asset_names = "image",
  stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1",
  collection = "naip",
  query_function = rsi_query_api,
  download_function = rsi_download_rasters,
  sign_function = sign_planetary_computer,
  rescale_bands = FALSE,
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "merge",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

get_alos_palsar_imagery(
  aoi,
  start_date,
  end_date,
  ...,
  pixel_x_size = 25,
  pixel_y_size = 25,
  asset_names = rsi::alos_palsar_band_mapping$planetary_computer_v1,
  stac_source = attr(asset_names, "stac_source"),
  collection = attr(asset_names, "collection_name"),
  query_function = attr(asset_names, "query_function"),
  download_function = attr(asset_names, "download_function"),
  sign_function = attr(asset_names, "sign_function"),
  rescale_bands = FALSE,
  item_filter_function = NULL,
  mask_band = attr(asset_names, "mask_band"),
  mask_function = attr(asset_names, "mask_function"),
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "median",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

get_dem(
  aoi,
  ...,
  start_date = NULL,
  end_date = NULL,
  pixel_x_size = 30,
  pixel_y_size = 30,
  asset_names = rsi::dem_band_mapping$planetary_computer_v1$`cop-dem-glo-30`,
  stac_source = attr(asset_names, "stac_source"),
  collection = attr(asset_names, "collection_name"),
  query_function = attr(asset_names, "query_function"),
  download_function = attr(asset_names, "download_function"),
  sign_function = attr(asset_names, "sign_function"),
  rescale_bands = FALSE,
  item_filter_function = NULL,
  mask_band = NULL,
  mask_function = NULL,
  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
  composite_function = "max",
  limit = 999,
  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30%", VSI_CACHE_SIZE =
    "10000000", 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 = "ALL_CPUS")
)

Arguments

`aoi`	An sf(c) object outlining the area of interest to get imagery for. Will be used to get the bounding box used for calculating metrics and the output data's CRS.
`start_date`, `end_date`	Character of length 1: The first and last date, respectively, of imagery to include in metrics calculations. Should be in YYYY-MM-DD format.
`pixel_x_size`, `pixel_y_size`	Numeric of length 1: size of pixels in x-direction (longitude / easting) and y-direction (latitude / northing).
`asset_names`	The names of the assets to download. If this vector has names, then the names of the vector are assumed to be the names of assets on the STAC server, which will be renamed to the elements of the vector in the final output.
`stac_source`	Character of length 1: the STAC URL to download imagery from.
`collection`	Character of length 1: the STAC collection to download images from.
`...`	Passed to `item_filter_function`.
`query_function`	A function that takes the output from `rstac::stac_search()` and executes the request. See `rsi_query_api()` and the `query_function` slots of sentinel1_band_mapping, sentinel2_band_mapping, and landsat_band_mapping.
`download_function`	A function that takes the output from `query_function` and downloads the assets attached to those items. See `rsi_download_rasters()` for an example.
`sign_function`	A function that takes the output from `query_function` and signs the item URLs, if necessary.
`rescale_bands`	Logical of length 1: If the STAC collection implements the `raster` STAC extension, and that extension includes `scale` and `offset` values, should this function attempt to automatically rescale the downloaded data?
`item_filter_function`	A function that takes the outputs of `query_function` (usually a `STACItemCollection`) and `...` and returns a filtered `STACItemCollection`. This is used, for example, to only download images from specific Landsat platforms.
`mask_band`	Character of length 1: The name of the asset in your STAC source to use to mask the data. Set to `NULL` to not mask. See the `mask_band` slots of sentinel1_band_mapping, sentinel2_band_mapping, and landsat_band_mapping.
`mask_function`	A function that takes a raster and returns a boolean raster, where `TRUE` pixels will be preserved and `FALSE` or `NA` pixels will be masked out. See `sentinel2_mask_function()`.
`output_filename`	The filename to write the output raster to. If `composite_function` is `NULL`, item datetimes will be appended to this in order to create unique filenames. If items do not have datetimes, a sequential ID will be appended instead.
`composite_function`	Character of length 1: The name of a function used to combine downloaded images into a single composite (i.e., to aggregate pixel values from multiple images into a single value). Must be one of of "sum", "mean", "median", "min", "max". Set to `NULL` to not composite (i.e., to rescale and save each individual file independently).
`limit`	an `integer` defining the maximum number of results to return. If not informed, it defaults to the service implementation.
`gdalwarp_options`	Options passed to `gdalwarp` through the `options` argument of `sf::gdal_utils()`. The same set of options are used for all downloaded data and the final output images; this means that some common options (for instance, `PREDICTOR=3`) may cause errors if bands are of varying data types.
`gdal_config_options`	Options passed to `gdalwarp` through the `config_options` argument of `sf::gdal_utils()`.
`platforms`	The names of Landsat satellites to download imagery from. These do not correspond to the `platforms` column in `spectral_indices()`; the default argument of `c("landsat-9", "landsat-8")` corresponds to the `Landsat-OLI` value in that column.

Value

output_filename, unchanged.

Usage Tips

It's often useful to buffer your aoi object slightly, on the order of 1-2 cell widths, in order to ensure that data is downloaded for your entire AOI even after accounting for any reprojection needed to compare your AOI to the data on the STAC server.

These functions allow for parallelizing downloads via future::plan(), and for user-controlled progress updates via progressr::handlers(). If there are fewer images to download than asset_names, then this function uses lapply() to iterate through images and future.apply::future_mapply() to iterate through downloading each asset. If there are more images than assets, this function uses future.apply::future_lapply() to iterate through images.

Downloading from Planetary Computer

Certain data sets in Planetary Computer require providing a subscription key. Even for non-protected data sets, providing a subscription key grants you higher rate limits and faster downloads. As such, it's a good idea to request a Planetary Computer account, then generate a subscription key. If you set the rsi_pc_key environment variable to your key (either primary or secondary; there is no difference), rsi will automatically use this key to sign all requests against Planetary Computer.

There are currently some challenges with certain Landsat images in Planetary Computer; please see https://github.com/microsoft/PlanetaryComputer/discussions/101 for more information on these images and their current status. These files may cause data downloads to fail.

Compositing

This function can either download all data that intersects with your spatiotemporal AOI as multiple files (if composite_function = NULL), or can be used to rescale band values, apply a mask function, and create a composite from the resulting files in a single function call. Each of these steps can be skipped by passing NULL to the corresponding argument.

Masks are applied to each downloaded asset separately. Rescaling is applied to the final composite after images are combined.

A number of the steps involved in creating composites – rescaling band values, running the mask function, masking images, and compositing images – currently rely on the terra package for raster calculations. This means creating larger composites, either in geographic or temporal dimension, may cause errors. It can be a good idea to tile your aoi using sf::st_make_grid() and iterate through the tiles to avoid these errors (and to make it easier to interrupt and restart a download job).

Rescaling

If rescale_bands is TRUE, then this function is able to use the scale and offset values in the bands field of the raster STAC extension. This was implemented originally to work with the Landsat collections in the Planetary Computer STAC catalogue, but hopefully will work automatically with other data sources as well. Note that Sentinel-2 data typically doesn't use this STAC extension, and so the returned data is typically not re-scaled; divide the downloaded band values by 10000 to get reflectance values in the expected values.

Sentinel-1 Data

The get_sentinel1_imagery() function is designed to download Sentinel-1 data from the Microsoft Planetary Computer STAC API. Both the GRD and RTC Sentinel-1 collections are supported. To download RTC data, set collection to sentinel-1-rtc, and supply your subscription key as an environment variable named rsi_pc_key (through, e.g., Sys.setenv() or your .Renviron file).

AlOS PALSAR Data

The get_alos_palsar_imagery() function is designed to download ALOS PALSAR annual mosaic data from the Microsoft Planetary Computer STAC API. Data are returned as a digital number (which is appropriate for some applications and indices). To convert to backscatter (decibels) use the following formula: 10 * log10(dn) - 83.0 where dn is the radar band in digital number.

Examples


aoi <- sf::st_point(c(-74.912131, 44.080410))
aoi <- sf::st_set_crs(sf::st_sfc(aoi), 4326)
aoi <- sf::st_buffer(sf::st_transform(aoi, 5070), 100)

get_stac_data(aoi,
  start_date = "2022-06-01",
  end_date = "2022-06-30",
  pixel_x_size = 30,
  pixel_y_size = 30,
  asset_names = c(
    "red", "blue", "green"
  ),
  stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1/",
  collection = "landsat-c2-l2",
  query_function = rsi_query_api,
  sign_function = sign_planetary_computer,
  mask_band = "qa_pixel",
  mask_function = landsat_mask_function,
  item_filter_function = landsat_platform_filter,
  platforms = c("landsat-9", "landsat-8")
)

# or, mostly equivalently (will download more bands):
landsat_image <- get_landsat_imagery(
  aoi,
  start_date = "2022-06-01",
  end_date = "2022-08-30"
)

[Package rsi version 0.2.1 Index]