From 981fb035bb4b650eabc468d0a7711a183239c9d6 Mon Sep 17 00:00:00 2001 From: Giacomo Falchetta <36954873+giacfalk@users.noreply.github.com> Date: Thu, 5 Feb 2026 17:08:33 +0100 Subject: [PATCH] Switch raster handling to stars --- CITATION.cff | 35 +++------ DESCRIPTION | 3 +- R/allocation.R | 141 ++++++++++++++++------------------ R/allocation_plot.R | 4 +- R/data.R | 8 +- R/friction.R | 81 +++++++++---------- R/get_cache_data.R | 9 +-- R/mask_raster_to_polygon.R | 45 +++++------ R/traveltime.R | 8 +- R/traveltime_plot.R | 4 +- R/traveltime_stats.r | 25 +++--- R/utils-checks.R | 2 +- R/utils.R | 67 ++++++++++++---- README.qmd | 4 +- codemeta.json | 23 ++---- data-raw/data_load.R | 10 +-- man-roxygen/params-demand.R | 2 +- man/allocation.Rd | 4 +- man/allocation_discrete.Rd | 4 +- man/friction.Rd | 2 +- man/mask_raster_to_polygon.Rd | 12 +-- man/naples_hot_days.Rd | 4 +- man/naples_population.Rd | 4 +- man/traveltime.Rd | 2 +- man/traveltime_stats.Rd | 2 +- paper/paper.md | 7 +- tests/testthat.R | 2 +- 27 files changed, 252 insertions(+), 262 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 67434d4..e3ff96f 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -305,19 +305,19 @@ references: year: '2025' version: '>= 4.4' - type: software - title: raster - abstract: 'raster: Geographic Data Analysis and Modeling' + title: stars + abstract: 'stars: Spatiotemporal Arrays, Raster and Vector Data Cubes' notes: Imports - url: https://rspatial.org/raster - repository: https://CRAN.R-project.org/package=raster + url: https://r-spatial.github.io/stars/ + repository: https://CRAN.R-project.org/package=stars authors: - - family-names: Hijmans - given-names: Robert J. - email: r.hijmans@gmail.com - orcid: https://orcid.org/0000-0001-5872-2872 + - family-names: Pebesma + given-names: Edzer + email: edzer.pebesma@uni-muenster.de + orcid: https://orcid.org/0000-0001-8049-7069 year: '2025' - doi: 10.32614/CRAN.package.raster - version: '>= 3.6.32' + doi: 10.32614/CRAN.package.stars + version: '>= 0.6.8' - type: software title: sf abstract: 'sf: Simple Features for R' @@ -357,20 +357,6 @@ references: address: Vienna, Austria year: '2025' version: '>= 4.4' -- type: software - title: terra - abstract: 'terra: Spatial Data Analysis' - notes: Imports - url: https://rspatial.org/ - repository: https://CRAN.R-project.org/package=terra - authors: - - family-names: Hijmans - given-names: Robert J. - email: r.hijmans@gmail.com - orcid: https://orcid.org/0000-0001-5872-2872 - year: '2025' - doi: 10.32614/CRAN.package.terra - version: '>= 1.8.86' - type: software title: tools abstract: 'R: A Language and Environment for Statistical Computing' @@ -549,4 +535,3 @@ references: year: '2025' doi: 10.32614/CRAN.package.testthat version: '>= 3.3.0' - diff --git a/DESCRIPTION b/DESCRIPTION index 17cf874..f5a8d5f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,11 +41,10 @@ Imports: malariaAtlas (>= 1.6.4), osmdata (>= 0.3.0), parallel (>= 4.4), - raster (>= 3.6.32), sf (>= 1.0.23), spatialEco (>= 2.0.3), + stars (>= 0.6.8), stats (>= 4.4), - terra (>= 1.8.86), tools (>= 4.4), zip (>= 2.3.3) Suggests: diff --git a/R/allocation.R b/R/allocation.R index 496cc96..de2e875 100644 --- a/R/allocation.R +++ b/R/allocation.R @@ -22,7 +22,7 @@ #' @template params-bb-area #' @template params-facilities #' @template params-traveltime-b -#' @param weights (optional) A raster with the weights for the demand (default: +#' @param weights (optional) A `stars` raster with the weights for the demand (default: #' `NULL`). #' @template params-objectiveminutes #' @template params-objectiveshare-a @@ -73,7 +73,7 @@ #' - `objective_share`: The value of the `objectiveshare` parameter used. #' - `facilities`: A [`sf`][sf::sf()] object with the newly allocated #' facilities. -#' - `travel_time`: A [`raster`][raster::raster()] RasterLayer object +#' - `travel_time`: A [`stars`][stars::st_as_stars()] object #' representing the travel time map with the newly allocated facilities. #' #' @family location-allocation functions @@ -152,16 +152,24 @@ apply_demand_transformation <- function(demand, weights, approach, exp_demand, e # Helper function: Find next facility location based on heuristic find_next_facility <- function(demand, heur, all_prev = NULL) { if (heur == "kd") { - return(raster::which.max(raster::raster(spatialEco::sp.kde( - x = sf::st_as_sf(raster::rasterToPoints(demand, spatial = TRUE)), - bw = terra::res(demand)*10, - res = terra::res(demand), + kde <- spatialEco::sp.kde( + x = sf::st_as_sf(demand, as_points = TRUE), + bw = min(stars::st_dimensions(demand)$x$delta, + stars::st_dimensions(demand)$y$delta) * 10, + res = min(stars::st_dimensions(demand)$x$delta, + stars::st_dimensions(demand)$y$delta), standardize = TRUE, scale.factor = 10000 - )) %>% raster::projectRaster(demand))) } + ) + + kde <- as_stars_if_needed(kde) |> + stars::st_warp(dest = demand, method = "near") + + return(stars_max_cell(kde)) + } # heur == "max" - raster::which.max(demand) + stars_max_cell(demand) } # Helper function: Initialize or accumulate new facilities @@ -178,14 +186,11 @@ accumulate_facility <- function(new_facilities, pos) { # Helper function: Apply travel time mask to demand mask_demand_by_traveltime <- function(demand, traveltime, objectiveminutes) { - demand |> - raster::overlay( - traveltime, - fun = function(x, y) { - x[y <= objectiveminutes] <- NA - x - } - ) + demand_values <- stars::st_values(demand) + traveltime_values <- stars::st_values(traveltime) + demand_values[traveltime_values <= objectiveminutes] <- NA + stars::st_values(demand) <- demand_values + demand } allocation <- function( @@ -204,14 +209,22 @@ allocation <- function( exp_demand = 1, exp_weights = 1 ) { - checkmate::assert_class(demand, "RasterLayer") + demand <- as_stars_if_needed(demand) + demand <- as_stars_if_needed(demand) + assert_stars(demand) assert_bb_area(bb_area) assert_facilities(facilities) assert_traveltime(traveltime, null_ok = TRUE) checkmate::assert_choice(mode, choices = c("walk", "fastest")) checkmate::assert_choice(dowscaling_model_type, choices = c("lm", "rf")) checkmate::assert_count(res_output, positive = TRUE) - checkmate::assert_class(weights, "RasterLayer", null.ok = TRUE) + if (!is.null(weights)) { + weights <- as_stars_if_needed(weights) + } + if (!is.null(weights)) { + weights <- as_stars_if_needed(weights) + } + assert_stars(weights, null_ok = TRUE) checkmate::assert_number(objectiveminutes, lower = 0) checkmate::assert_number(objectiveshare, lower = 0, upper = 1) checkmate::assert_choice(heur, choices = c("max", "kd")) @@ -252,17 +265,14 @@ allocation <- function( traveltime_raster_outer[[1]] |> mask_raster_to_polygon(bb_area) - raster::crs(traveltime) <- - "+proj=longlat +datum=WGS84 +no_defs +type=crs" - - traveltime <- raster::projectRaster(traveltime, demand) - - raster::crs(traveltime) <- "+proj=longlat +datum=WGS84 +no_defs +type=crs" + sf::st_crs(traveltime) <- 4326 + traveltime <- stars::st_warp(traveltime, dest = demand, method = "near") + sf::st_crs(traveltime) <- 4326 # Apply demand transformation using helper demand <- apply_demand_transformation(demand, weights, approach, exp_demand, exp_weights, bb_area) - totalpopconstant <- demand |> raster::cellStats("sum", na.rm = TRUE) + totalpopconstant <- stars_sum(demand) demand <- mask_demand_by_traveltime(demand, traveltime, objectiveminutes) @@ -274,11 +284,7 @@ allocation <- function( iter <- iter + 1 all <- find_next_facility(demand, heur, all_prev = NULL) - - pos <- - demand |> - raster::xyFromCell(all) |> - as.data.frame() + pos <- all$xy new_facilities <- accumulate_facility(new_facilities, pos) @@ -300,21 +306,18 @@ allocation <- function( traveltime_raster_new <- traveltime_raster_outer[[2]][[3]] |> gdistance::accCost(xy_matrix) |> - raster::crop(raster::extent(demand)) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> - raster::projectRaster(demand) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> + as_stars_if_needed() |> + stars::st_crop(sf::st_bbox(demand)) |> + sf::st_set_crs(4326) |> + stars::st_warp(dest = demand, method = "near") |> + sf::st_set_crs(4326) |> mask_raster_to_polygon(bb_area) demand <- mask_demand_by_traveltime(demand, traveltime_raster_new, objectiveminutes) k <- demand |> - raster::cellStats("sum", na.rm = TRUE) |> + stars_sum() |> magrittr::divide_by(totalpopconstant) k_save[iter] <- k @@ -330,7 +333,9 @@ allocation <- function( if (k < (1 - objectiveshare)) { break } else if (k == k_save[iter - 1]) { - raster::values(demand)[all] <- NA + demand_values <- stars::st_values(demand) + demand_values[all$index] <- NA + stars::st_values(demand) <- demand_values } } @@ -382,28 +387,22 @@ apply_demand_transformation <- function(demand, weights, approach, exp_demand, e # Helper: Mask demand by travel time threshold mask_demand_by_traveltime <- function(demand, traveltime_raster, objectiveminutes) { - demand |> - raster::overlay( - traveltime_raster, - fun = function(x, y) { - x[y <= objectiveminutes] <- NA - x - } - ) + demand_values <- stars::st_values(demand) + traveltime_values <- stars::st_values(traveltime_raster) + demand_values[traveltime_values <= objectiveminutes] <- NA + stars::st_values(demand) <- demand_values + demand } # Helper: Calculate travel time raster from cost surface calculate_traveltime_raster <- function(traveltime_raster_outer, xy_matrix, demand, bb_area) { traveltime_raster_outer[[2]][[3]] |> gdistance::accCost(xy_matrix) |> - raster::crop(raster::extent(demand)) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> - raster::projectRaster(demand) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> + as_stars_if_needed() |> + stars::st_crop(sf::st_bbox(demand)) |> + sf::st_set_crs(4326) |> + stars::st_warp(dest = demand, method = "near") |> + sf::st_set_crs(4326) |> mask_raster_to_polygon(bb_area) } @@ -471,9 +470,9 @@ initialize_traveltime <- function(traveltime, demand, bb_area, mode, dowscaling_ return(list(traveltime_outer = traveltime, is_new = FALSE)) } - traveltime_new <- demand |> - raster::`values<-`(objectiveminutes + 1) |> - mask_raster_to_polygon(bb_area) + traveltime_new <- demand + stars::st_values(traveltime_new) <- objectiveminutes + 1 + traveltime_new <- traveltime_new |> mask_raster_to_polygon(bb_area) friction <- bb_area |> friction( @@ -504,7 +503,7 @@ allocation_discrete <- function( exp_weights = 1, par = FALSE ) { - checkmate::assert_class(demand, "RasterLayer") + assert_stars(demand) assert_bb_area(bb_area) checkmate::assert_multi_class(candidate, c("sf", "sfc")) assert_facilities(facilities, null_ok = TRUE) @@ -515,7 +514,7 @@ allocation_discrete <- function( checkmate::assert_choice(mode, choices = c("walk", "fastest")) checkmate::assert_choice(dowscaling_model_type, choices = c("lm", "rf")) checkmate::assert_count(res_output, positive = TRUE) - checkmate::assert_class(weights, "RasterLayer", null.ok = TRUE) + assert_stars(weights, null_ok = TRUE) checkmate::assert_number(objectiveminutes, lower = 0) checkmate::assert_number(objectiveshare, lower = 0, upper = 1, null.ok = TRUE) checkmate::assert_choice(approach, choices = c("norm", "absweights")) @@ -567,16 +566,12 @@ allocation_discrete <- function( ) traveltime_raster_outer <- traveltime_init$traveltime_outer - totalpopconstant <- demand |> raster::cellStats("sum", na.rm = TRUE) + totalpopconstant <- stars_sum(demand) traveltime_raster_outer[[1]] <- traveltime_raster_outer[[1]] |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> - raster::projectRaster(demand) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) + sf::st_set_crs(4326) |> + stars::st_warp(dest = demand, method = "near") |> + sf::st_set_crs(4326) demand <- mask_demand_by_traveltime(demand, traveltime_raster_outer[[1]], objectiveminutes) demand_raster_bk <- demand @@ -615,7 +610,7 @@ allocation_discrete <- function( demand_rasterio <- mask_demand_by_traveltime(demand_rasterio, traveltime_raster_new, objectiveminutes) demand_rasterio |> - raster::cellStats("sum", na.rm = TRUE) |> + stars_sum() |> magrittr::divide_by(totalpopconstant) } @@ -645,7 +640,7 @@ allocation_discrete <- function( demand <- mask_demand_by_traveltime(demand, traveltime_raster_new, objectiveminutes) k <- demand |> - raster::cellStats("sum", na.rm = TRUE) |> + stars_sum() |> magrittr::divide_by(totalpopconstant) out <- list( @@ -713,7 +708,7 @@ allocation_discrete <- function( demand_rasterio <- mask_demand_by_traveltime(demand_rasterio, traveltime_raster_new, objectiveminutes) demand_rasterio |> - raster::cellStats("sum", na.rm = TRUE) |> + stars_sum() |> magrittr::divide_by(totalpopconstant) } @@ -743,7 +738,7 @@ allocation_discrete <- function( demand <- mask_demand_by_traveltime(demand, traveltime_raster_new, objectiveminutes) k <- demand |> - raster::cellStats("sum", na.rm = TRUE) |> + stars_sum() |> magrittr::divide_by(totalpopconstant) cli::cli_alert_info( diff --git a/R/allocation_plot.R b/R/allocation_plot.R index 5a751ae..69605c2 100644 --- a/R/allocation_plot.R +++ b/R/allocation_plot.R @@ -108,7 +108,7 @@ allocation_plot <- function( max_limit <- allocation |> magrittr::extract2("travel_time") |> - raster::values() |> + stars::st_values() |> max(na.rm = TRUE) plot <- @@ -118,7 +118,7 @@ allocation_plot <- function( data = allocation |> #nolint magrittr::extract2("travel_time") |> mask_raster_to_polygon(bb_area) |> - raster::as.data.frame(xy = TRUE) |> + as.data.frame(xy = TRUE) |> stats::na.omit() ) + ggplot2::geom_sf( diff --git a/R/data.R b/R/data.R index 7cf73e3..215ce12 100644 --- a/R/data.R +++ b/R/data.R @@ -64,7 +64,7 @@ #' #' @description #' -#' A [`RasterLayer`][raster::raster()] object representing the population +#' A [`stars`][stars::st_as_stars()] object representing the population #' density in the city of Naples, Italy. #' #' The dataset is based on the Global Human @@ -73,7 +73,7 @@ #' ([GHS-POP]( #' https://human-settlement.emergency.copernicus.eu/download.php?ds=pop)). #' -#' @format A [`RasterLayer`][raster::raster()] object with 1 layer. +#' @format A [`stars`][stars::st_as_stars()] object with 1 layer. #' #' @source Global Human Settlement Layer #' ([GHSL](https://human-settlement.emergency.copernicus.eu)). @@ -85,7 +85,7 @@ #' #' @description #' -#' A [`RasterLayer`][raster::raster()] object representing the number of hot +#' A [`stars`][stars::st_as_stars()] object representing the number of hot #' days in the city of Naples, Italy. #' #' This 100-meter resolution heat hazard map shows the number of days with @@ -94,7 +94,7 @@ #' above 25 °C during 2008–2017, based on simulations from the #' [UrbClim](https://www.urban-climate.eu/model) model. #' -#' @format A [`RasterLayer`][raster::raster()] object with 1 layer. +#' @format A [`stars`][stars::st_as_stars()] object with 1 layer. #' #' @source [UrbClim](https://www.urban-climate.eu/model). #' diff --git a/R/friction.R b/R/friction.R index bd44536..6d0d1fc 100644 --- a/R/friction.R +++ b/R/friction.R @@ -36,7 +36,7 @@ #' #' @return An [invisible][base::invisible] [`list`][base::list] with the #' following elements: -#' - `friction_layer`: A [`RasterLayer`][raster::raster()] object with the +#' - `friction_layer`: A [`stars`][stars::st_as_stars()] object with the #' friction surface layer. #' - `transition_matrix`: A [`TransitionLayer`][gdistance::transition()] with #' the transition matrix for cost-distance calculations. @@ -91,16 +91,15 @@ friction <- function( checkmate::assert_file_exists(file) - friction_layer <- file |> terra::rast() + friction_layer <- file |> stars::read_stars() friction_layer <- friction_layer |> - terra::crop( + stars::st_crop( bb_area |> sf::st_transform(sf::st_crs(friction_layer)) |> - terra::ext() - ) |> - raster::raster() + sf::st_bbox() + ) } else { if (mode == "fastest") { dataset <- "Accessibility__201501_Global_Travel_Speed_Friction_Surface" @@ -120,7 +119,7 @@ friction <- function( malariaAtlas::getRaster( extent = matrix(sf::st_bbox(bb_area), ncol = 2), ) |> - raster::raster() + stars::st_as_stars() } } @@ -138,65 +137,51 @@ friction <- function( r <- bb_area |> sf::st_transform(crs = 3395) |> - raster::extent() |> - raster::raster(res = res_output, crs = sf::st_crs(3395)$proj4string) + sf::st_bbox() |> + stars::st_as_stars(dx = res_output, dy = res_output) |> + sf::st_set_crs(3395) streets <- x$osm_lines |> sf::st_transform(crs = 3395) |> sf::st_buffer(res_output) |> - terra::vect() |> - terra::rasterize( - r |> terra::rast(), - background = NA, - fun = "sum" - ) |> - raster::raster() + stars::st_rasterize( + template = r + ) d <- streets |> - terra::rast() |> - terra::project(y = sf::st_crs(4326)$proj4string) + stars::st_warp(crs = sf::st_crs(4326), method = "near") - terra::values(d) <- ifelse( - terra::values(d) < 0, - 0, - terra::values(d) - ) - - terra::values(d) <- ifelse( - is.na(terra::values(d)), - 0, - terra::values(d) - ) + d_values <- stars::st_values(d) + d_values <- ifelse(is.na(d_values) | d_values < 0, 0, d_values) + stars::st_values(d) <- d_values d <- d |> - raster::raster() |> - raster::crop(friction_layer) + stars::st_crop(sf::st_bbox(friction_layer)) d_2 <- d - d <- raster::stack(d, d_2) - - names(d) <- paste0("l", 1:raster::nlayers(d)) + d <- c(d, d_2) + names(d) <- paste0("l", seq_along(d)) min_iter <- 2 max_iter <- 10 p_train <- 0.5 - if (length(unique(raster::values(friction_layer))) == 1) { - raster::values(friction_layer) <- + if (length(unique(stars::st_values(friction_layer))) == 1) { + stars::st_values(friction_layer) <- stats::runif( - min = unique(raster::values(friction_layer)) * 0.9, - max = unique(raster::values(friction_layer)) * 1.1, - n = length(raster::values(friction_layer)) + min = unique(stars::st_values(friction_layer)) * 0.9, + max = unique(stars::st_values(friction_layer)) * 1.1, + n = length(stars::st_values(friction_layer)) ) } res_rf <- dissever::dissever( - coarse = friction_layer, # stack of fine resolution covariates - fine = d, # coarse resolution raster + coarse = as_raster_if_needed(friction_layer), # stack of fine resolution covariates + fine = as_raster_if_needed(d), # coarse resolution raster method = dowscaling_model_type, # regression method used for disseveration p = p_train, # proportion of pixels sampled for training regression model min_iter = min_iter, # minimum iterations @@ -204,13 +189,16 @@ friction <- function( verbose = TRUE ) - raster::values(res_rf$map) <- ifelse( - raster::values(res_rf$map) <= 0, - min(raster::values(friction_layer), na.rm = TRUE), - raster::values(res_rf$map) + res_rf_map <- as_stars_if_needed(res_rf$map) + res_rf_values <- stars::st_values(res_rf_map) + res_rf_values <- ifelse( + res_rf_values <= 0, + min(stars::st_values(friction_layer), na.rm = TRUE), + res_rf_values ) + stars::st_values(res_rf_map) <- res_rf_values - friction_layer <- res_rf$map + friction_layer <- res_rf_map } else { friction_layer <- friction_layer } @@ -219,6 +207,7 @@ friction <- function( transition_matrix <- friction_layer |> + as_raster_if_needed() |> # RAM intensive, can be very slow for large areas. gdistance::transition(\(x) 1 / mean(x), 16) diff --git a/R/get_cache_data.R b/R/get_cache_data.R index 4436373..243a6a2 100644 --- a/R/get_cache_data.R +++ b/R/get_cache_data.R @@ -25,15 +25,14 @@ get_cache_data <- function(dataset, bb_area) { if (!is.na(file)) { cli::cli_progress_step("Using cached surface friction data") - out <- file |> terra::rast() + out <- file |> stars::read_stars() out |> - terra::crop( + stars::st_crop( bb_area |> sf::st_transform(sf::st_crs(out)) |> - terra::ext() - ) |> - raster::raster() + sf::st_bbox() + ) } else { cli::cli_progress_step("Downloading and caching surface friction data") diff --git a/R/mask_raster_to_polygon.R b/R/mask_raster_to_polygon.R index f042218..355d314 100644 --- a/R/mask_raster_to_polygon.R +++ b/R/mask_raster_to_polygon.R @@ -1,19 +1,19 @@ -#' Mask a `RasterLayer` object to a `sf` object +#' Mask a `stars` object to a `sf` object #' #' @description #' -#' This function rapidly masks a [`RasterLayer`][raster::raster()] object to +#' This function rapidly masks a [`stars`][stars::st_as_stars()] object to #' a [`sf`][sf::st_as_sf()] object. #' -#' @param ras A [`RasterLayer`][raster::raster()] object. -#' @param mask A [`RasterLayer`][raster::raster()] or [`sf`][sf::st_as_sf()] -#' object. +#' @param ras A [`stars`][stars::st_as_stars()] object. +#' @param mask A [`stars`][stars::st_as_stars()] or [`sf`][sf::st_as_sf()] +#' object. #' @param inverse A [`logical`][base::logical()] flag. If `TRUE`, the mask #' is inverted (default: `FALSE`). -#' @param updatevalue The value to update the [`RasterLayer`][raster::raster()] +#' @param updatevalue The value to update the [`stars`][stars::st_as_stars()] #' object with (default: `NA`). #' -#' @return A masked [`RasterLayer`][raster::raster()] object. +#' @return A masked [`stars`][stars::st_as_stars()] object. #' #' @family utility functions #' @keywords masking @@ -27,8 +27,13 @@ mask_raster_to_polygon <- function( inverse = FALSE, updatevalue = NA ) { - checkmate::assert_class(ras, "RasterLayer") - checkmate::assert_multi_class(mask, c("Raster", "sf")) + ras <- as_stars_if_needed(ras) + if (inherits(mask, "Raster")) { + mask <- as_stars_if_needed(mask) + } + + checkmate::assert_class(ras, "stars") + checkmate::assert_multi_class(mask, c("stars", "sf")) checkmate::assert_flag(inverse) checkmate::assert_atomic(updatevalue, len = 1) @@ -51,21 +56,13 @@ mask_raster_to_polygon <- function( mask <- mask |> - sf::st_crop( - c( - xmin = raster::xmin(ras), - ymin = raster::ymin(ras), - xmax = raster::xmax(ras), - ymax = raster::ymax(ras) - ) - ) |> - sf::st_cast() |> - terra::vect() |> - suppressWarnings() + sf::st_crop(sf::st_bbox(ras)) |> + sf::st_cast() } - ras |> - terra::rast() |> - terra::mask(mask, inverse = inverse) |> - raster::raster() + if (inverse) { + stars::st_mask(ras, mask, invert = TRUE, update = updatevalue) + } else { + stars::st_mask(ras, mask, update = updatevalue) + } } diff --git a/R/traveltime.R b/R/traveltime.R index bdccaec..aff177d 100644 --- a/R/traveltime.R +++ b/R/traveltime.R @@ -13,7 +13,7 @@ #' #' @return An [invisible][base::invisible] [`list`][base::list] with the #' following elements: -#' - `travel_time`: A [`RasterLayer`][raster::raster()] object with the travel +#' - `travel_time`: A [`stars`][stars::st_as_stars()] object with the travel #' time map. #' - `friction`: A [`list`][base::list] with the outputs of the #' [`friction()`][friction] function. @@ -83,8 +83,10 @@ traveltime <- function( travel_time <- friction_data[[3]] |> gdistance::accCost(points) |> - mask_raster_to_polygon(bb_area) |> - raster::`crs<-`(value = "+proj=longlat +datum=WGS84 +no_defs +type=crs") + as_stars_if_needed() |> + mask_raster_to_polygon(bb_area) + + sf::st_crs(travel_time) <- 4326 list( travel_time = travel_time, diff --git a/R/traveltime_plot.R b/R/traveltime_plot.R index b9b050f..8c16919 100644 --- a/R/traveltime_plot.R +++ b/R/traveltime_plot.R @@ -67,7 +67,7 @@ traveltime_plot <- function( traveltime |> magrittr::extract2("travel_time") |> mask_raster_to_polygon(bb_area) |> - raster::as.data.frame(xy = TRUE) |> + as.data.frame(xy = TRUE) |> stats::na.omit() plot <- @@ -101,7 +101,7 @@ traveltime_plot <- function( max_travel_time <- traveltime |> magrittr::extract2("travel_time") |> - raster::values() |> + stars::st_values() |> max(na.rm = TRUE) if (max_travel_time > contour_traveltime) { diff --git a/R/traveltime_stats.r b/R/traveltime_stats.r index 5520eeb..9983b2c 100644 --- a/R/traveltime_stats.r +++ b/R/traveltime_stats.r @@ -57,7 +57,8 @@ traveltime_stats <- function( print = TRUE ) { assert_traveltime(traveltime) - checkmate::assert_class(demand, "RasterLayer") + demand <- as_stars_if_needed(demand) + assert_stars(demand) checkmate::assert_numeric(breaks, min.len = 1) checkmate::assert_number(objectiveminutes, lower = 1) checkmate::assert_flag(print) @@ -67,21 +68,19 @@ traveltime_stats <- function( traveltime_values <- demand_values <- P15_cumsum <- th <- NULL # nolint end - raster::crs(demand) <- "+proj=longlat +datum=WGS84 +no_defs +type=crs" + sf::st_crs(demand) <- 4326 - data_curve <- + traveltime_warped <- traveltime[[1]] |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" - ) |> - raster::projectRaster(demand) |> - raster::`crs<-`( - value = "+proj=longlat +datum=WGS84 +no_defs +type=crs" + sf::st_set_crs(4326) |> + stars::st_warp(dest = demand, method = "near") + + data_curve <- + data.frame( + traveltime_values = as.vector(stars::st_values(traveltime_warped)), + demand_values = as.vector(stars::st_values(demand)) ) |> - raster::values() |> - data.frame(raster::values(demand)) |> - stats::na.omit(data_curve) |> - magrittr::set_colnames(c("traveltime_values", "demand_values")) |> + stats::na.omit() |> dplyr::arrange(traveltime_values) |> dplyr::as_tibble() |> dplyr::mutate( diff --git a/R/utils-checks.R b/R/utils-checks.R index f64947e..5c3adc2 100644 --- a/R/utils-checks.R +++ b/R/utils-checks.R @@ -128,7 +128,7 @@ assert_minimal_coverage <- function( null_ok = FALSE ) { assert_traveltime(traveltime) - checkmate::assert_class(demand, "RasterLayer") + assert_stars(demand) checkmate::assert_number(objectiveminutes, lower = 0) checkmate::assert_number(objectiveshare, lower = 0, upper = 1, null.ok = TRUE) checkmate::assert_flag(null_ok) diff --git a/R/utils.R b/R/utils.R index fb83573..3e3ccdf 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,19 +1,57 @@ -#' This function is a workaround to ensure that the `raster` package -#' is loaded when plotting `RasterLayer` objects from within the -#' `locationallocation` package (e.g., `plot(naples_population)`). -#' -#' Simply importing `plot` with `importFrom()` does not work because -#' S4 method dispatch requires the `raster` package to be loaded for -#' RasterLayer objects to be plotted correctly. -#' #' @noRd #' @export plot <- function(x, ...) { - if (inherits(x, "RasterLayer")) { - raster::plot(x, ...) - } else { - graphics::plot(x, ...) + graphics::plot(x, ...) +} + +as_stars_if_needed <- function(x) { + if (inherits(x, "stars")) { + return(x) + } + + if (!requireNamespace("stars", quietly = TRUE)) { + cli::cli_abort("The {.pkg stars} package is required to handle raster data.") } + + stars::st_as_stars(x) +} + +as_raster_if_needed <- function(x) { + if (inherits(x, "Raster")) { + return(x) + } + + if (!requireNamespace("raster", quietly = TRUE)) { + cli::cli_abort( + "The {.pkg raster} package is required for cost-distance calculations." + ) + } + + methods::as(x, "Raster") +} + +assert_stars <- function(x, null_ok = FALSE) { + if (isTRUE(null_ok) && is.null(x)) { + return(invisible(TRUE)) + } + + checkmate::assert_class(x, "stars") +} + +stars_sum <- function(x) { + sum(stars::st_values(x), na.rm = TRUE) +} + +stars_max_cell <- function(x) { + value_name <- names(x)[1] + values <- as.vector(stars::st_values(x)) + values_no_na <- ifelse(is.na(values), -Inf, values) + max_index <- which.max(values_no_na) + + data <- as.data.frame(x, xy = TRUE) + max_row <- data[which.max(data[[value_name]]), c("x", "y"), drop = FALSE] + + list(index = max_index, xy = max_row) } facilities_coordinates <- function(facilities, bb_area = NULL) { @@ -47,8 +85,9 @@ points_to_matrix <- function(points, n = NULL) { } normalize_raster <- function(r) { - r_min <- raster::cellStats(r, stat = "min") - r_max <- raster::cellStats(r, stat = "max") + r_values <- stars::st_values(r) + r_min <- min(r_values, na.rm = TRUE) + r_max <- max(r_values, na.rm = TRUE) (r - r_min) / (r_max - r_min) } diff --git a/README.qmd b/README.qmd index b4046f7..73519ca 100644 --- a/README.qmd +++ b/README.qmd @@ -14,7 +14,7 @@ library(badger) library(ggplot2) library(here) library(magrittr) -library(raster) +library(stars) library(sf) library(tidyterra) @@ -93,7 +93,7 @@ These datasets include the coordinates of public drinking water fountains (blue #| echo: false library(ggplot2) -library(raster) +library(stars) library(sf) library(tidyterra) diff --git a/codemeta.json b/codemeta.json index ff977eb..f19612c 100644 --- a/codemeta.json +++ b/codemeta.json @@ -321,16 +321,16 @@ }, "15": { "@type": "SoftwareApplication", - "identifier": "raster", - "name": "raster", - "version": ">= 3.6.32", + "identifier": "stars", + "name": "stars", + "version": ">= 0.6.8", "provider": { "@id": "https://cran.r-project.org", "@type": "Organization", "name": "Comprehensive R Archive Network (CRAN)", "url": "https://cran.r-project.org" }, - "sameAs": "https://CRAN.R-project.org/package=raster" + "sameAs": "https://CRAN.R-project.org/package=stars" }, "16": { "@type": "SoftwareApplication", @@ -365,25 +365,12 @@ "version": ">= 4.4" }, "19": { - "@type": "SoftwareApplication", - "identifier": "terra", - "name": "terra", - "version": ">= 1.8.86", - "provider": { - "@id": "https://cran.r-project.org", - "@type": "Organization", - "name": "Comprehensive R Archive Network (CRAN)", - "url": "https://cran.r-project.org" - }, - "sameAs": "https://CRAN.R-project.org/package=terra" - }, - "20": { "@type": "SoftwareApplication", "identifier": "tools", "name": "tools", "version": ">= 4.4" }, - "21": { + "20": { "@type": "SoftwareApplication", "identifier": "zip", "name": "zip", diff --git a/data-raw/data_load.R b/data-raw/data_load.R index 100e1ca..8489479 100644 --- a/data-raw/data_load.R +++ b/data-raw/data_load.R @@ -1,5 +1,5 @@ library(here) -library(raster) +library(stars) library(sf) library(usethis) @@ -12,9 +12,9 @@ library(usethis) #' representing the administrative boundary of the city of Naples, Italy. #' - `naples_fountains`: A [`sf`][sf::st_as_sf()] point geometry object #' representing the water fountains in the city of Naples, Italy. -#' - `naples_population`: A [`Raster`][raster::raster()] object representing the +#' - `naples_population`: A [`stars`][stars::st_as_stars()] object representing the #' population density in the city of Naples, Italy. -#' - `naples_hot_days`: A [`Raster`][raster::raster()] object representing the +#' - `naples_hot_days`: A [`stars`][stars::st_as_stars()] object representing the #' number of hot days in the city of Naples, Italy. #' #' @return An invisible `NULL`. This function is used for its side effect. @@ -35,12 +35,12 @@ data_load <- function() { ), list( file = here::here("data-raw", "naples-population.tif"), - fun = \(x) raster::readAll(raster::raster(x)), + fun = stars::read_stars, name = "naples_population" ), list( file = here::here("data-raw", "naples-hot-days.tif"), - fun = \(x) raster::readAll(raster::raster(x)), + fun = stars::read_stars, name = "naples_hot_days" ) ) diff --git a/man-roxygen/params-demand.R b/man-roxygen/params-demand.R index 8faa23a..1ec7788 100644 --- a/man-roxygen/params-demand.R +++ b/man-roxygen/params-demand.R @@ -1,2 +1,2 @@ -#' @param demand A [`RasterLayer`][raster::raster()] object with the +#' @param demand A [`stars`][stars::st_as_stars()] object with the #' demand layer (e.g. population density). diff --git a/man/allocation.Rd b/man/allocation.Rd index f7f96c2..9666fdd 100644 --- a/man/allocation.Rd +++ b/man/allocation.Rd @@ -22,7 +22,7 @@ allocation( ) } \arguments{ -\item{demand}{A \code{\link[raster:raster]{RasterLayer}} object with the +\item{demand}{A \code{\link[stars:st_as_stars]{stars}} object with the demand layer (e.g. population density).} \item{bb_area}{A \code{\link[sf:st_as_sf]{sf}} boundary box object with the area @@ -107,7 +107,7 @@ of demand that remains unmet after allocating the new facilities. \item \code{objective_share}: The value of the \code{objectiveshare} parameter used. \item \code{facilities}: A \code{\link[sf:sf]{sf}} object with the newly allocated facilities. -\item \code{travel_time}: A \code{\link[raster:raster]{raster}} RasterLayer object +\item \code{travel_time}: A \code{\link[stars:st_as_stars]{stars}} object representing the travel time map with the newly allocated facilities. } } diff --git a/man/allocation_discrete.Rd b/man/allocation_discrete.Rd index 9b2765f..e3bb48d 100644 --- a/man/allocation_discrete.Rd +++ b/man/allocation_discrete.Rd @@ -25,7 +25,7 @@ allocation_discrete( ) } \arguments{ -\item{demand}{A \code{\link[raster:raster]{RasterLayer}} object with the +\item{demand}{A \code{\link[stars:st_as_stars]{stars}} object with the demand layer (e.g. population density).} \item{bb_area}{A \code{\link[sf:st_as_sf]{sf}} boundary box object with the area @@ -125,7 +125,7 @@ of demand that remains unmet after allocating the new facilities. \item \code{objective_share}: The value of the \code{objectiveshare} parameter used. \item \code{facilities}: A \code{\link[sf:sf]{sf}} object with the newly allocated facilities. -\item \code{travel_time}: A \code{\link[raster:raster]{raster}} RasterLayer object +\item \code{travel_time}: A \code{\link[stars:st_as_stars]{stars}} object representing the travel time map with the newly allocated facilities. } } diff --git a/man/friction.Rd b/man/friction.Rd index 0c6ac90..bde5e72 100644 --- a/man/friction.Rd +++ b/man/friction.Rd @@ -53,7 +53,7 @@ function will use this file instead of downloading the friction data.} An \link[base:invisible]{invisible} \code{\link[base:list]{list}} with the following elements: \itemize{ -\item \code{friction_layer}: A \code{\link[raster:raster]{RasterLayer}} object with the +\item \code{friction_layer}: A \code{\link[stars:st_as_stars]{stars}} object with the friction surface layer. \item \code{transition_matrix}: A \code{\link[gdistance:transition]{TransitionLayer}} with the transition matrix for cost-distance calculations. diff --git a/man/mask_raster_to_polygon.Rd b/man/mask_raster_to_polygon.Rd index bd94ace..7aa2dbe 100644 --- a/man/mask_raster_to_polygon.Rd +++ b/man/mask_raster_to_polygon.Rd @@ -2,27 +2,27 @@ % Please edit documentation in R/mask_raster_to_polygon.R \name{mask_raster_to_polygon} \alias{mask_raster_to_polygon} -\title{Mask a \code{RasterLayer} object to a \code{sf} object} +\title{Mask a \code{stars} object to a \code{sf} object} \usage{ mask_raster_to_polygon(ras, mask, inverse = FALSE, updatevalue = NA) } \arguments{ -\item{ras}{A \code{\link[raster:raster]{RasterLayer}} object.} +\item{ras}{A \code{\link[stars:st_as_stars]{stars}} object.} -\item{mask}{A \code{\link[raster:raster]{RasterLayer}} or \code{\link[sf:st_as_sf]{sf}} +\item{mask}{A \code{\link[stars:st_as_stars]{stars}} or \code{\link[sf:st_as_sf]{sf}} object.} \item{inverse}{A \code{\link[base:logical]{logical}} flag. If \code{TRUE}, the mask is inverted (default: \code{FALSE}).} -\item{updatevalue}{The value to update the \code{\link[raster:raster]{RasterLayer}} +\item{updatevalue}{The value to update the \code{\link[stars:st_as_stars]{stars}} object with (default: \code{NA}).} } \value{ -A masked \code{\link[raster:raster]{RasterLayer}} object. +A masked \code{\link[stars:st_as_stars]{stars}} object. } \description{ -This function rapidly masks a \code{\link[raster:raster]{RasterLayer}} object to +This function rapidly masks a \code{\link[stars:st_as_stars]{stars}} object to a \code{\link[sf:st_as_sf]{sf}} object. } \examples{ diff --git a/man/naples_hot_days.Rd b/man/naples_hot_days.Rd index 5417c0e..d70943d 100644 --- a/man/naples_hot_days.Rd +++ b/man/naples_hot_days.Rd @@ -5,7 +5,7 @@ \alias{naples_hot_days} \title{Hot days in the city of Naples, Italy} \format{ -A \code{\link[raster:raster]{RasterLayer}} object with 1 layer. +A \code{\link[stars:st_as_stars]{stars}} object with 1 layer. } \source{ \href{https://www.urban-climate.eu/model}{UrbClim}. @@ -14,7 +14,7 @@ A \code{\link[raster:raster]{RasterLayer}} object with 1 layer. naples_hot_days } \description{ -A \code{\link[raster:raster]{RasterLayer}} object representing the number of hot +A \code{\link[stars:st_as_stars]{stars}} object representing the number of hot days in the city of Naples, Italy. This 100-meter resolution heat hazard map shows the number of days with diff --git a/man/naples_population.Rd b/man/naples_population.Rd index 29139f0..203793f 100644 --- a/man/naples_population.Rd +++ b/man/naples_population.Rd @@ -5,7 +5,7 @@ \alias{naples_population} \title{Population density in the city of Naples, Italy} \format{ -A \code{\link[raster:raster]{RasterLayer}} object with 1 layer. +A \code{\link[stars:st_as_stars]{stars}} object with 1 layer. } \source{ Global Human Settlement Layer @@ -15,7 +15,7 @@ Global Human Settlement Layer naples_population } \description{ -A \code{\link[raster:raster]{RasterLayer}} object representing the population +A \code{\link[stars:st_as_stars]{stars}} object representing the population density in the city of Naples, Italy. The dataset is based on the Global Human diff --git a/man/traveltime.Rd b/man/traveltime.Rd index 593fcdd..a162af1 100644 --- a/man/traveltime.Rd +++ b/man/traveltime.Rd @@ -57,7 +57,7 @@ function will use this file instead of downloading the friction data.} An \link[base:invisible]{invisible} \code{\link[base:list]{list}} with the following elements: \itemize{ -\item \code{travel_time}: A \code{\link[raster:raster]{RasterLayer}} object with the travel +\item \code{travel_time}: A \code{\link[stars:st_as_stars]{stars}} object with the travel time map. \item \code{friction}: A \code{\link[base:list]{list}} with the outputs of the \code{\link[=friction]{friction()}} function. diff --git a/man/traveltime_stats.Rd b/man/traveltime_stats.Rd index 54635bc..36cc623 100644 --- a/man/traveltime_stats.Rd +++ b/man/traveltime_stats.Rd @@ -16,7 +16,7 @@ traveltime_stats( \item{traveltime}{A \code{\link[base:list]{list}} with the output of the \code{\link[=traveltime]{traveltime()}} function.} -\item{demand}{A \code{\link[raster:raster]{RasterLayer}} object with the +\item{demand}{A \code{\link[stars:st_as_stars]{stars}} object with the demand layer (e.g. population density).} \item{objectiveminutes}{(optional) A number indicating the target travel time diff --git a/paper/paper.md b/paper/paper.md index 0bc88e2..d54eb2c 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -235,8 +235,7 @@ perform location-allocation spatial optimization at a high spatial-resolution (particularly useful in urban-scale applications). A set of reporting functions and graphical outputs are pre-calculated as part of the package. The package further relies on *malariaAtlas* and -*gdistance* functions, as well as on *raster*, *terra*, and *sf* object -classes.\ +*gdistance* functions, as well as on *stars* and *sf* object classes.\ - **Friction layers:** Malaria Atlas friction surfaces: walking or fastest mode. These datasets provide 1-km resolution global gridded @@ -249,7 +248,7 @@ classes.\ - **Point locations for candidate facilities:** Optional, a point simple feature geometry object ($sf$). -- **Demand weights:** A raster. It can be, for instance, population +- **Demand weights:** A raster (handled as a *stars* object). It can be, for instance, population counts per location (grid cell) - optionally in a weighted form by specifying the $weights$ argument to another raster of a function of several rasters. This would allow using a risk framework where the @@ -585,4 +584,4 @@ is thankful to Ahmed T. Hammad for the previous joint work in is also grateful to his scientific affiliations CMCC and IIASA for their continuous support. -# References \ No newline at end of file +# References diff --git a/tests/testthat.R b/tests/testthat.R index 9760e03..07bf686 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -11,7 +11,7 @@ test_check("locationallocation") # - Reduce cyclomatic complexity. # - Add unit tests. # - Add test coverage (Codecov). -# - Move to the r-spatial framework or adopt `terra` +# - Move to the r-spatial framework or adopt `stars` # (the first option aligns better with tidyverse principles). # - Standardize parameter names (for example, `updatevalue` to `new_value`). # - Merge `traveltime` and `traveltime_discrete`.