fix reading raster tiles with temporal dimensions

NEWS.md

@@ -8,6 +8,7 @@
 ## Bug Fixes
 
 - `biome` and `ecoregions` now properly handle 0-length tibbles (#196)
+- logic for reading raster tiles with temporal dimensions improved
 
 ## New features
 
@@ -17,6 +18,8 @@
 ## Internal
 
 - `calc_indicators()` now includes a check for 0-length tibbles (#199)
+- .read_raster_source now uses a single logic to cover all cases (e.g. single tiles,
+  tiled rasters with and without temporal dimension, single temporal rasters)
 
 
 # mapme.biodiversity 0.4.0

R/calc_indicators.R

@@ -123,58 +123,32 @@
 
 
 .read_raster_source <- function(x, tindex) {
-  all_bboxes <- lapply(1:nrow(tindex), function(i) paste(as.numeric(st_bbox(tindex[i, ])), collapse = " "))
-  is_stacked <- length(unique(unlist(all_bboxes))) == 1
 
-  if (is_stacked) { # current resource/extent all have the same bounding box
+  geoms <- tindex[["geom"]]
+  unique_geoms <- unique(geoms)
+  grouped_geoms <- match(geoms, unique_geoms)
+  names(grouped_geoms) <- tindex[["location"]]
+  grouped_geoms <- sort(grouped_geoms)
 
-    filenames <- basename(tindex$location)
-    out <- terra::rast(tindex$location)
-    names(out) <- filenames
-  } else {
-    is_unique <- length(unique(unlist(all_bboxes))) == nrow(tindex)
-
-    if (is_unique) { # all tiles have a different bounding box
-      target_files <- tindex$location[unlist(st_intersects(x, tindex))]
+  n_tiles <- length(unique(grouped_geoms))
+  n_timesteps <- unique(table(grouped_geoms))
 
-      if (length(target_files) == 0) {
-        warning("No intersection with resource.")
-        return(NULL)
-      } else if (length(target_files) == 1) {
-        out <- terra::rast(target_files)
-      } else {
-        # create a vrt for multiple targets
-        vrt_name <- tempfile("vrt", fileext = ".vrt")
-        out <- terra::vrt(target_files, filename = vrt_name)
-      }
-    } else { # some tiles share the same bboxes, and others do not, needs proper merging
-      # We assume here that the tiles present in tileindex have a temporal dimension.
-      # Thus each timestep should end up in its own layer. Different tiles from
-      # the same timestep should be spatially merged. We want to avoid merging
-      # different tile from different timesteps. We thus assume some regularity
-      # in how the name of a raster file expresses its temporal dimension.
-      # With this assumption, we can expect the files in tindex to be ordered.
-      # Thus we retrive the index of all files sharing the same bbox and assume
-      # that they belong to different timesteps. The files in between these
-      # indices thus belong to the previous timestep and we can merge these
-      # as a vrt and later join the bands. We always assign the name of the
-      # first file as the layername.
-      unique_bboxes <- unique(unlist(all_bboxes))
-      layer_index <- which(all_bboxes == unique_bboxes[[1]])
-      temporal_gap <- layer_index[2] - layer_index[1] - 1
-      out <- lapply(layer_index, function(j) {
-        target_files <- tindex$location[j:(j + temporal_gap)]
-        org_filename <- basename(target_files[1])
-        filename <- tools::file_path_sans_ext(org_filename)
-        vrt_name <- tempfile(pattern = sprintf("vrt_%s.vrt", filename))
-        tmp <- terra::vrt(target_files, filename = vrt_name)
-        names(tmp) <- org_filename
-        tmp
-      })
-      out <- do.call(c, out)
-    }
+  if (length(n_timesteps) > 1) {
+    stop("Did not find equal number of tiles per timestep.")
   }
 
+  out <- lapply(1:n_timesteps, function(i){
+    index <- rep(FALSE, n_timesteps)
+    index[i] <- TRUE
+    filenames <- names(grouped_geoms[index])
+    layer_name <- tools::file_path_sans_ext(basename(filenames[1]))
+    vrt_name <- tempfile(pattern = sprintf("vrt_%s", layer_name), fileext = ".vrt")
+    tmp <- terra::vrt(filenames, filename = vrt_name)
+    names(tmp) <- layer_name
+    tmp
+  })
+  out <- do.call(c, out)
+
   # crop the source to the extent of the current polygon
   cropped <- try(terra::crop(out, terra::vect(x)))
   if (inherits(cropped, "try-error")) {
@@ -211,7 +185,7 @@
   n_rows <- sapply(results[index_tbl], nrow)
   if (any(n_rows == 0)) {
     stop(paste("0-length tibbles returned for some assets.\n",
-    "Make sure the indicator function returns NA if it cannot be calculated for an asset."))
+               "Make sure the indicator function returns NA if it cannot be calculated for an asset."))
   }
 
   # case all assets returned tibbles

R/get_resources.R

@@ -110,16 +110,7 @@ get_resources <- function(x, resources, ...) {
   if (selected_resource[[resource]]$type == "raster") {
     tindex_file <- file.path(rundir, paste0("tileindex_", resource, ".gpkg"))
     if (file.exists(tindex_file)) file.remove(tindex_file)
-
-    footprints <- lapply(unique(downloaded_files), function(file) {
-      tmp <- rast(file)
-      footprint <- st_as_sf(st_as_sfc(st_bbox(tmp)))
-      st_geometry(footprint) <- "geom"
-      footprint$location <- sources(tmp)
-      footprint
-    })
-
-    footprints <- do.call(rbind, footprints)
+    footprints <- .make_footprints(downloaded_files)
     write_sf(footprints, dsn = tindex_file)
     downloaded_files <- tindex_file
   }
@@ -131,3 +122,14 @@ get_resources <- function(x, resources, ...) {
   attributes(x) <- atts
   x
 }
+
+.make_footprints <- function(raster_files) {
+  footprints <- lapply(unique(raster_files), function(file) {
+    tmp <- rast(file)
+    footprint <- st_as_sf(st_as_sfc(st_bbox(tmp)))
+    st_geometry(footprint) <- "geom"
+    footprint$location <- sources(tmp)
+    footprint
+  })
+  do.call(rbind, footprints)
+}

tests/testthat/test-calc_indicator.R

@@ -1,7 +1,7 @@
 test_that("calc_indicator works", {
   aoi <- read_sf(
     system.file("extdata", "gfw_sample.gpkg",
-      package = "mapme.biodiversity"
+                package = "mapme.biodiversity"
     )
   )
 
@@ -92,7 +92,7 @@ test_that("calc_indicator works", {
 test_that("Parallelization works", {
   aoi <- read_sf(
     system.file("extdata", "gfw_sample.gpkg",
-      package = "mapme.biodiversity"
+                package = "mapme.biodiversity"
     )
   )
 
@@ -216,7 +216,7 @@ test_that(".bind_assets works correctly", {
 test_that(".prep works correctly", {
   x <- read_sf(
     system.file("extdata", "gfw_sample.gpkg",
-      package = "mapme.biodiversity"
+                package = "mapme.biodiversity"
     )
   )
 
@@ -295,3 +295,46 @@ test_that(".prep works correctly", {
     "Resource type 'sth' currently not supported"
   )
 })
+
+
+test_that(".read_raster_source works correctly", {
+
+  dummy <- terra::rast()
+  dummy_splitted <- aggregate(dummy, fact = c(ceiling(nrow(dummy) / 4), ceiling(ncol(dummy) / 4)))
+  dummy_splitted[] <- 1:16
+  polys <- terra::as.polygons(dummy_splitted) %>% st_as_sf()
+  dummies <- lapply(1:nrow(polys), function(i) crop(dummy_splitted, polys[i, ]))
+  temp_loc <- tempfile()
+  dir.create(temp_loc, showWarnings = FALSE)
+  purrr::walk(1:length(dummies), function(i) {
+    writeRaster(dummies[[i]], filename = file.path(temp_loc, paste0("2000_tile_", i, ".tif")))
+    writeRaster(dummies[[i]], filename = file.path(temp_loc, paste0("2001_tile_", i, ".tif")))
+  })
+
+  files <- list.files(temp_loc, full.names = TRUE)
+  footprints <- .make_footprints(files)
+  x <- st_bbox(dummy) %>% st_as_sfc() %>% st_as_sf()
+  extent <- c(-180, 180, -90, 90)
+  names(extent) <- c("xmin", "xmax", "ymin", "ymax")
+
+  tiled_temporal <- .read_raster_source(x, footprints)
+  expect_equal(names(tiled_temporal), c("2000_tile_1", "2001_tile_1"))
+  expect_equal(as.vector(ext(tiled_temporal)), extent)
+
+  tiled <- .read_raster_source(x, footprints[grep("2001", footprints$location), ])
+  expect_equal(names(tiled), "2001_tile_1")
+  expect_equal(as.vector(ext(tiled)), extent)
+
+  temporal <- .read_raster_source(x, footprints[grep("tile_12.tif", footprints$location), ])
+  extent[c(1:4)] <- c(90, 180, -45, 0)
+  expect_equal(names(temporal), c("2000_tile_12", "2001_tile_12"))
+  expect_equal(as.vector(ext(temporal)), extent)
+
+  single <- .read_raster_source(x, footprints[grep("2000_tile_10.tif", footprints$location), ])
+  extent[c(1:4)] <- c(-90, 0, -45, 0)
+  expect_equal(names(single), "2000_tile_10")
+  expect_equal(as.vector(ext(single)), extent)
+
+  expect_error(.read_raster_source(x, footprints[1:24, ]))
+
+})