Downloading STAC data using rsi when you’ve got a geographic CRS or don’t want a composite.

Methods for straying slightly off the happy path: just say NULL.
R
Tutorials
Spatial
geospatial data
R packages
Author

Mike Mahoney

Published

November 21, 2023

A quick post today, inspired by a GitHub issue.

I’ve been working recently on the new rsi package which helps you download, reproject, resample, mask, rescale, and composite data from STAC APIs.1 The standard function interface does all of these steps: it grabs all the relevant files from your STAC source, reprojects them to match your AOI and desired resolution, masks and rescales the component files, and then merges them into a composite:

library(rsi)
future::plan("multisession")

aoi <- sf::st_point(c(-74.912131, 44.080410)) |> 
  sf::st_sfc() |> 
  sf::st_set_crs(4326) |> 
  sf::st_transform(3857) |> 
  sf::st_buffer(1000)

start_date <- "2022-06-01"
end_date <- "2022-07-01"

get_landsat_imagery(
  aoi = aoi,
  start_date = start_date,
  end_date = end_date,
  output_filename = tempfile(fileext = ".tif")
) |> 
  terra::rast() |> 
  terra::plot()

What if you want to skip some of these steps? For instance, if you try to call get_landsat_imagery() with an AOI in geographic coordinates, you’ll get a warning (likely followed by an error) saying that you’re asking to resample the data to 30 degree pixels, which is probably not what you wanted:

try(
  get_landsat_imagery(
    aoi = sf::st_transform(aoi, 4326),
    start_date = start_date,
    end_date = end_date,
    output_filename = tempfile(fileext = ".tif")
  )
)
Warning: The default pixel size arguments are intended for use with projected AOIs, but `aoi` appears to be in geographic coordinates.
ℹ Pixel X size: 30. Pixel Y size: 30.
ℹ These dimensions will be interpreted in the same units as `aoi` (likely degrees), which may cause errors.
Warning in CPL_gdalwarp(source, destination, options, oo, doo, config_options,
: GDAL Error 1: Attempt to create 0x0 dataset is illegal,sizes must be larger
than zero.
Warning: Failed to download LC08_L2SP_015029_20220617_02_T1 from
2022-06-17T15:45:03.055481Z
Warning in CPL_gdalwarp(source, destination, options, oo, doo, config_options,
: GDAL Error 1: Attempt to create 0x0 dataset is illegal,sizes must be larger
than zero.
Warning: Failed to download LC09_L2SP_015029_20220609_02_T2 from
2022-06-09T15:44:23.649712Z
Warning in CPL_gdalwarp(source, destination, options, oo, doo, config_options,
: GDAL Error 1: Attempt to create 0x0 dataset is illegal,sizes must be larger
than zero.
Warning: Failed to download LC08_L2SP_015029_20220601_02_T1 from
2022-06-01T15:44:51.569374Z
Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Error
4: /tmp/Rtmp7E1urK/filedddb8dfac1e.tif: No such file or directory
Error : [rast] file does not exist: /tmp/Rtmp7E1urK/filedddb8dfac1e.tif

That’s coming from the resampling step of the function’s workflow. Can we just skip that?

Short answer: yes! If we pass NULL to the pixel_*_size arguments, we’ll skip the resampling stage and instead download our data in its native resolution:

get_landsat_imagery(
  aoi = sf::st_transform(aoi, 4326),
  start_date = start_date,
  end_date = end_date,
  pixel_x_size = NULL,
  pixel_y_size = NULL,
  output_filename = tempfile(fileext = ".tif")
) |> 
  terra::rast() |> 
  terra::plot()

This is a pattern throughout the rsi API design: if you want to skip something, pass NULL to the relevant argument. For instance (and this is where the GitHub issue comes in), if you want to not composite and instead download all the images within your spatiotemporal area of interest, we can pass NULL to the composite_function argument to skip compositing. I’ll also skip masking by passing NULL to the mask_function argument, because otherwise a handful of these images are entirely masked out:

get_landsat_imagery(
  aoi = aoi,
  start_date = start_date,
  end_date = end_date,
  output_filename = tempfile(fileext = ".tif"),
  composite_function = NULL,
  mask_function = NULL # otherwise half of these images are blank
) |> 
  lapply(terra::rast) |> 
  lapply(terra::plot)

[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

Hopefully this helps people use rsi to only perform the data wrangling steps they want!

Footnotes

  1. And calculate spectral indices from these data, and wrangle multiple rasters into a multi-band VRT – it’s a pretty neat package if I do say so myself.↩︎