diff --git a/NEWS.md b/NEWS.md index 6580bb2a..9583e30f 100644 --- a/NEWS.md +++ b/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 diff --git a/R/calc_indicators.R b/R/calc_indicators.R index d0d816cf..3d00487d 100644 --- a/R/calc_indicators.R +++ b/R/calc_indicators.R @@ -123,58 +123,32 @@ calc_indicators <- function(x, indicators, ...) { .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 @@ calc_indicators <- function(x, indicators, ...) { 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 diff --git a/R/get_resources.R b/R/get_resources.R index 822badd0..17440f31 100644 --- a/R/get_resources.R +++ b/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) +} diff --git a/tests/testthat/test-calc_indicator.R b/tests/testthat/test-calc_indicator.R index e48df74b..dcb14099 100644 --- a/tests/testthat/test-calc_indicator.R +++ b/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, ])) + +})