use gdalcubes to subset and rescale single rasters on the fly (gdalwarp)

gdalcubes_raster.R

url <- "/vsicurl/https://minio.carlboettiger.info/public-biodiversity/carbon/Irrecoverable_Carbon_2010/Irrecoverable_C_Total_2010.tif"

library(gdalcubes)
library(spData)
library(sf)
library(stars)
india <- spData::world |> dplyr::filter(name_long=="India")


# subset by bounding box, reproject coordinates, rescale dx/dy
gdal_read <- function(x, box = NULL, crs = NULL, dx= NULL, dy = NULL) {
  
  # We can use the native box, crs, and/or scale if not provided. 
  # Though if all are null, we might be better with read_stars...
  proxy <- stars::read_stars(x)
  
  
  # dx,dy can't be non-NULL if crs is not given
  if( (!is.null(dx) || !is.null(dy)) && is.null(crs)) {
    stop("if dx/dy values are given, a crs must be specified!")
  }

  # guess crs and dx/dy
  if(is.null(crs)) {
    crs <- st_crs(proxy)
    id <- gsub('.*ID\\["EPSG",(\\d{4})\\]\\]$', "\\1", crs$wkt)
    crs <- paste0("EPSG:", id)
    ## dx 
    dim <- st_dimensions(proxy)
    dx <- dim$x$delta
    dy <- abs(dim$y$delta)
  }
  
  # if crs is given but dx/dy aren't, the native CRS could be different.
  # in principle we could transform ourselves, but lazy ask user to do it.
 if( (is.null(dx) || is.null(dy)) && !is.null(crs)) {
    stop(paste("please specify the dx/dy values in CRS", crs))
  }

  ## bounding box can also be a polygon.  Handles re-projecting it too.
  if(is.null(box)) {
    box <- st_bbox(proxy)
  } else {
    box <- sf::st_bbox(box) |> 
      st_as_sfc() |>
      sf::st_transform(st_crs(crs)) |>
      st_bbox()
  }
  # NOTE: temporal domain is dummy var here
  view <- cube_view(srs = crs,
                    extent = list(t0 = "2010-01-01", t1 = "2010-01-01",
                                  left = box[1], right = box[3],
                                  top = box[4], bottom = box[2]),
                    dx = dx, dy=dy, dt = "P1Y")
  
  img <- gdalcubes::create_image_collection(x, date_time = "2010-01-01")
  raster_cube(img, view)
}


# Read in a raster url with polygon subset & CRS transform into meters
z <- gdal_read(url, india, "EPSG:24378", 10000, 10000)
plot(z)

# Subset but don't downscale:
z <- gdal_read(url, india)
plot(z)

here's using the new by_util in terra

url <- "/vsicurl/https://minio.carlboettiger.info/public-biodiversity/carbon/Irrecoverable_Carbon_2010/Irrecoverable_C_Total_2010.tif"

india <- spData::world |> dplyr::filter(name_long=="India")

library(terra)  ## we need >= 	1.7-46 for by_util
## has yucky extent 
target <- project(rast(india), "EPSG:24378")
res(target) <- c(10000, 10000)
## a bit convoluted, but if we want nice clean alignment to the resolution (I have the logic in vaster but I still have to think hard)
##target <- rast(align(ext(project(rast(india), "EPSG:24378")), 10000), crs = "EPSG:24378", res = 10000)

## we've set to the polygon extent/crs and set a resolution
project(rast(url), target, by_util = TRUE)


## crop(rast(url), india) on its own is probably as fast as this anyway
## subset but don't downscale (rast twice so we only have the grid specification)
project(rast(url), crop(rast(rast(url)), india), by_util = TRUE)