Add gsw

DESCRIPTION

@@ -51,6 +51,7 @@ Collate:
     'calc_gsw_occurrence.R'
     'calc_gsw_recurrence.R'
     'calc_gsw_seasonality.R'
+    'calc_gsw_time_series.R'
     'calc_gsw_transitions.R'
     'calc_hfp.R'
     'get_resources.R'
@@ -84,6 +85,7 @@ Collate:
     'get_gfw_treecover.R'
     'get_gmw.R'
     'get_gsw.R'
+    'get_gsw_time_series.R'
     'get_hfp.R'
     'get_ipbes_biomes.R'
     'get_key_biodiversity_areas.R'

NAMESPACE

@@ -15,6 +15,7 @@ export(calc_gsw_change)
 export(calc_gsw_occurrence)
 export(calc_gsw_recurrence)
 export(calc_gsw_seasonality)
+export(calc_gsw_time_series)
 export(calc_gsw_transitions)
 export(calc_humanfootprint)
 export(calc_indicators)
@@ -55,6 +56,7 @@ export(get_global_surface_water_recurrence)
 export(get_global_surface_water_seasonality)
 export(get_global_surface_water_transitions)
 export(get_gmw)
+export(get_gsw_time_series)
 export(get_humanfootprint)
 export(get_ipbes_biomes)
 export(get_irr_carbon)

NEWS.md

@@ -18,6 +18,7 @@
   to write a GTiff file and errors if unsuccessful (#335)
 
 - new resources: 
+  - `get_gsw_time_series()`
   - `get_biodiversity_intactness_index()`
   - `get_chelsa()`
   - `get_ipbes_biomes()`
@@ -32,6 +33,8 @@
   - `calc_biodiversity_intactness_index ()`
   - `calc_exposed_population()`
   - `calc_humanfootprint()`
+  - `calc_exposed_population()`
+  - `calc_gsw_time_series()` 
   - `calc_precipitation_chelsa()`
   - `calc_vul_carbon()`, `calc_man_carbon()`, and `calc_irr_carbon()`
 

R/calc_gsw_time_series.R

@@ -0,0 +1,110 @@
+#' Calculate Global Surface Water Time Series
+#'
+#' This function calculates the total area of the global surface water time
+#' series data, separated by the following classes:
+#'
+#' - No Observation: It was not possible to determine whether a pixel was water
+#' (this may be the case for frozen areas or during the polar night in extreme
+#' latitudes).
+#' - Permanent Water: Water was detected in twelve months per year or in a
+#' combination of permanent and no observation.
+#' - Seasonal Water: Water and no water was detected.
+#' - No Water: No Water was detected.
+#'
+#' The required resources for this indicator are:
+#'  - [gsw_time_series_resource]
+#'
+#' @name gsw_time_series_indicator
+#' @keywords indicator
+#' @format A function returning a tibble with time series of global surface
+#' water data classes.
+#' @include register.R
+#' @export
+#' @examples
+#' \dontshow{
+#' mapme.biodiversity:::.copy_resource_dir(file.path(tempdir(), "mapme-data"))
+#' }
+#' \dontrun{
+#' library(sf)
+#' library(mapme.biodiversity)
+#'
+#' outdir <- file.path(tempdir(), "mapme-data")
+#' dir.create(outdir, showWarnings = FALSE)
+#'
+#' mapme_options(
+#'   outdir = outdir,
+#'   verbose = FALSE
+#' )
+#'
+#' aoi <- read_sf(
+#'   system.file("extdata", "shell_beach_protected_area_41057_B.gpkg",
+#'               package = "mapme.biodiversity"
+#' ))
+#' aoi <- get_resources(aoi, get_gsw_time_series (years = 2000:2001))
+#' aoi <- calc_indicators(calc_gsw_time_series())
+#' aoi <- portfolio_long(aoi)
+#'
+#' aoi
+#' }
+calc_gsw_time_series <- function() {
+  check_namespace("exactextractr")
+
+  function(x = NULL,
+           gsw_time_series,
+           name = "gsw_timeseries",
+           mode = "asset",
+           aggregation = "sum",
+           verbose = mapme_options()[["verbose"]]) {
+
+    if (is.null(gsw_time_series)) {
+      return(NULL)
+    }
+
+    gsw_time_series <- mask(
+      gsw_time_series,
+      x
+    )
+
+    gsw_time_series <- clamp(
+      gsw_time_series,
+      lower = 0,
+      upper = 3,
+      values = FALSE
+    )
+
+    idx_yr <- as.numeric(gregexpr("yearlyClassification", names(gsw_time_series))) + 20
+    years_resource <- substr(names(gsw_time_series), idx_yr, idx_yr + 3)
+
+    names(gsw_time_series) <- years_resource
+    coverage_fractions <- exactextractr::exact_extract(gsw_time_series, x, "frac", coverage_area = TRUE)
+    if(nlyr(gsw_time_series) == 1) {
+      names(coverage_fractions) <- paste(names(coverage_fractions), years_resource, sep = ".")
+    }
+    x_total_area <- as.numeric(st_area(x)) / 10000
+
+    results <- purrr::map_df(names(coverage_fractions), function(colname) {
+      year <- substr(colname, nchar(colname) - 3, nchar(colname))
+      variable <- switch(
+        substr(colname, 6, 6),
+        "0" = "no_observations",
+        "1" = "not_water",
+        "2" = "seasonal_water",
+        "3" = "permanent_water"
+      )
+      tibble::tibble(
+        datetime = as.POSIXct(paste0(year, "-01-01T00:00:00Z")),
+        variable = variable,
+        unit = "ha",
+        value = coverage_fractions[[colname]] * x_total_area
+      )
+    })
+
+    return(results)
+  }
+}
+
+register_indicator(
+  name = "gsw_time_series",
+  description = "Global Surface Water - Yearly Time Series area estimation of water classes.",
+  resources = "gsw_time_series"
+)

R/get_gsw_time_series.R

@@ -0,0 +1,111 @@
+#' Helper function to download Global Surface Water (GSW) yearly time series
+#' data
+#'
+#' This function constructs  the necessary data URLs for a given data set,
+#' version and polygon and downloads them for further processing with the
+#' mapme.biodiversity package.
+#'
+#' @details
+#' The available surface water classes for a given pixel are the following:
+#'
+#' - No Observation: It was not possible to determine whether a pixel was water
+#' (this may be the case for frozen areas or during the polar night in extreme
+#' latitudes).
+#' - Permanent Water: Water was detected in twelve months per year or in a
+#' combination of permanent and no observation.
+#' - Seasonal Water: Water and no water was detected.
+#' - No Water: No Water was detected.
+#'
+#' @name gsw_time_series_resource
+#' @keywords resource
+#' @returns A function that returns a character vector of file paths.
+#' @references
+#' * Global Surface Water Explorer: \url{https://global-surface-water.appspot.com/}
+#' * Data Users Guide: \url{https://storage.cloud.google.com/global-surface-water/downloads_ancillary/DataUsersGuidev2021.pdf}
+#' * Research Article: \url{https://www.nature.com/articles/nature20584}
+#' @source Raw Data: \url{https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GSWE/YearlyClassification/LATEST/tiles/}
+#' @param years Numeric vector of years to process between 1984 and 2021.
+#' Default: `1984:2021`.
+#' @param version Version of the data set to process. Available options are
+#' (`VER1-0`, `VER2-0`, `VER3-0`, `VER4-0`, `VER5-0` and `LATEST`) Default:
+#' `LATEST`. Choosing `LATEST` will result in the latest available version.
+#' @export
+get_gsw_time_series <- function(years, version = "LATEST") {
+  version <- unique(version)
+  available_versions = c("VER1-0", "VER2-0", "VER3-0", "VER4-0", "VER5-0",
+                         "LATEST")
+  stopifnot(version %in% available_versions)
+  if(version == "LATEST") {
+    version <- "VER5-0"
+  }
+
+  years <- unique(years)
+  years <- years[!is.na(as.numeric(years))]
+  available_years <- 1984:2021
+  years <- check_available_years(years, available_years, "gsw_time_series")
+
+  function(x,
+           name = "gsw_time_series",
+           type = "raster",
+           outdir = mapme_options()[["outdir"]],
+           verbose = mapme_options()[["verbose"]]) {
+    # make the gsw grid and construct urls for intersecting tiles
+    baseurl <- sprintf(
+      "/vsicurl/https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GSWE/YearlyClassification/%s/tiles/",
+      version
+    )
+    grid_gsw <- make_global_grid(
+      xmin = -180, xmax = 180, dx = 10,
+      ymin = -60, ymax = 80, dy = 10
+    )
+    tile_ids <- unique(unlist(sf::st_intersects(x, grid_gsw)))
+    if (length(tile_ids) == 0) {
+      stop("The extent of the portfolio does not intersect with the GSW grid.",
+           call. = FALSE
+      )
+    }
+    ids <- sapply(tile_ids, function(n) .get_gsw_ts_tile_id(grid_gsw[n, ]))
+
+    source_combinations <- expand.grid(years = years, ids = ids)
+    urls <- purrr::map2_chr(source_combinations$years, source_combinations$ids, \(year, tile_url_id) {
+      # Warning: file naming system is different for 2021
+      separator <- ifelse(year == 2021, "_", "-")
+      sprintf(
+        "%syearlyClassification%s/yearlyClassification%s%s%s.tif",
+        baseurl, year, year, separator, tile_url_id
+      )
+    })
+
+    fps <- grid_gsw[expand.grid(years, tile_ids = tile_ids)$tile_ids, ]
+    fps[["source"]] <- urls
+    make_footprints(fps,
+                    filenames = paste0(version, "_", basename(fps$source)),
+                    what = "raster",
+                    co = c("-co", "INTERLEAVE=BAND", "-co", "COMPRESS=LZW", "-ot", "Byte")
+    )
+  }
+}
+
+.get_gsw_ts_tile_id <- function(tile) {
+  min_x <- st_bbox(tile)[1]
+  max_y <- st_bbox(tile)[4]
+
+  x_idx <- which(seq(-180, 170, by = 10) == min_x) - 1
+  y_idx <- which(seq(80, -50, by = -10) == max_y) - 1
+
+  tile_id <- sprintf(
+    "%010d-%010d",
+    y_idx * 40000,
+    x_idx * 40000
+  )
+
+  return(tile_id)
+}
+
+register_resource(
+  name = "gsw_time_series",
+  description = "Global Surface Water - Yearly Time Series",
+  licence = "https://global-surface-water.appspot.com/download",
+  source = "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GSWE/YearlyClassification/LATEST/tiles/",
+  type = "raster"
+)

README.md

@@ -72,6 +72,7 @@ Below is a list of the resources currently supported by
 | global_surface_water_seasonality | Global Surface Water - Seasonality of water occurrrence                                                                    | <https://www.copernicus.eu/en/access-data>                 |
 | global_surface_water_transitions | Global Surface Water - Transition classes                                                                                  | <https://www.copernicus.eu/en/access-data>                 |
 | gmw                              | Global Mangrove Watch - Vector data of mangrove extent                                                                     | CC BY 4.0                                                  |
+| gsw_time_series                  | Global Surface Water - Yearly Time Series                                                                                  | <https://global-surface-water.appspot.com/download>        |
 | humanfootprint                   | Time series on human pressures on natural ecosystems.                                                                      | CC BY 4.0                                                  |
 | ipbes_biomes                     | Global Assessment Report on Biodiversity and Ecosystem Services division of the earth’s surface into biomes and anthromes. | CC 4.0                                                     |
 | irr_carbon                       | Amount of carbon irrecoverably lost by a typical land use conversion event until mid-century.                              | CC NC 4.0                                                  |
@@ -107,6 +108,7 @@ Next, is a list of supported indicators.
 | gsw_occurrence                | Areal statistic of surface water based on occurrence threshold                 |
 | gsw_recurrence                | Areal statistic of surface water based on reccurence threshold                 |
 | gsw_seasonality               | Areal statistic of surface water by seasonality                                |
+| gsw_time_series               | Global Surface Water - Yearly Time Series area estimation of water classes.    |
 | gsw_transitions               | Areal statistics of surface water grouped by transition class                  |
 | humanfootprint                | Statistics of the human footprint data set per polygon.                        |
 | ipbes_biomes                  | Area distibution of IBPES biomes within a polygon.                             |

data-raw/gsw_time_series.R

@@ -0,0 +1,31 @@
+library(dplyr)
+library(mapme.indicators)
+library(terra)
+
+years <- 2000:2001
+outdir <- system.file("res", package = "mapme.biodiversity")
+mapme_options(outdir = outdir)
+
+x <- read_sf(
+  system.file("extdata", "shell_beach_protected_area_41057_B.gpkg",
+              package = "mapme.biodiversity"
+  )
+) %>%
+get_resources(get_gsw_time_series(years = years))
+
+gsw_time_series <- prep_resources(x)
+gsw_time_series <- gsw_time_series [[1]]
+
+outdir <- "inst/res/gsw_time_series"
+dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
+
+gsw_time_series <- mask(gsw_time_series, x)
+
+for(lyr_id in seq_len(nlyr(gsw_time_series))) {
+  lyr_name <- names(gsw_time_series) [lyr_id]
+  lyr_name <- sub("VER5-0_yearlyClassification", "v5_", lyr_name)
+  fname <- paste0(file.path(outdir, lyr_name), ".tif")
+  gsw_time_series %>%
+    subset(lyr_id) %>%
+    writeRaster(fname)
+}