R/get_site_ODS.R

In oggioniale/ReLTER: An Interface for the eLTER Community

#' Acquire various raster layers from
#' \href{https://maps.opendatascience.eu/}{ODS Europe}
#' and crops to an eLTER site boundary.
#' @description `r lifecycle::badge("stable")`
#' Download and return a SpatRaster object containing the requested
#' dataset from \href{https://maps.opendatascience.eu/}{ODS},
#' cropped to an eLTER site boundary, which is obtained from the DEIMS-SDR API.
#' @param deimsid  A `character`. The DEIMS ID of the site from
#' DEIMS-SDR website. DEIMS ID information
#' \href{https://deims.org/docs/deimsid.html}{here}.
#' @param dataset A `character`. The requested dataset. One of:
#' "landcover", "clc2018", "osm_buildings", "natura2000",
#' "ndvi_spring", "ndvi_summer", "ndvi_autumn", "ndvi_winter", "ndvi_trend",
#' "forest_broadleaf", "forest_mixed", "forest_coniferous".
#' Default is "landcover".
#' @details Supported datasets from the ODS repository include:
#' Landcover: Land-cover class according to the highest probability,
#'    generated by a spatiotemporal ensemble-ML model. 30 m. resolution
#' CLC2018: Corine land cover rasterized to 100m spatial resolution
#'    and provided by Copernicus Land Monitoring Service.
#' OSM buildings: Buildings according to OSM polygons
#'    and the Copernicus impervious build-up layer (2018),
#'    aggregated and rasterized first to 10m spatial resolution
#'    and after downsampled to 30m by spatial average.
#' Natura2000: Protected areas rasterized from NATURA 2000
#'    (A, B and C site categories)
#'    and OSM (IUCN Ia, IUCN Ib, IUCN 2, IUCN 3, IUCN 4, IUCN 5, IUCN 6
#'    and others categories),
#'    first to 10m spatial resolution and after downsampled
#'    to 30m by spatial average.
#'    The overlap areas are indicated in a new category.
#' NDVI:  NDVI time-series,
#'    derived from the Landsat quarterly temporal composites
#'    NDVI Trend from 2000 - 2019 as OLS regression
#' Forests: Broadleaf, coniferous or mixed forests
#' All datasets are georeferenced to the
#' EPSG:3035 coordinate reference system.
#' and all except clc2018 have 30 meters resolution
#' @return The function returns a SpatRaster object (from the `terra` package)
#' of the requested dataset, cropped to the site boundaries
#' The user should save the raster to disk, if necessary.
#' i.e. writeRaster(ds_site, "site_dataset.tif")
#' @author Micha Silver, phD (2020) \email{silverm@@post.bgu.ac.il}
#' @author Alessandro Oggioni, phD (2020) \email{oggioni.a@@irea.cnr.it}
#' @importFrom dplyr case_when
#' @importFrom sf st_transform
#' @importFrom terra mask crop vect rast crs plot
#' @importFrom Rdpack reprompt
#' @references
#'   \insertRef{dplyrR}{ReLTER}
#'
#'   \insertRef{sfR}{ReLTER}
#'
#'   \insertRef{terraR}{ReLTER}
#' @export
#' @examples
#'  \dontrun{
#' # Landcover for Angelo Mosso
#' siteLandcover <- get_site_ODS(
#'   deimsid = "https://deims.org/17210eba-d832-4759-89fa-9ff127cbdf6e",
#'   dataset = "landcover"
#' )
#' siteLandcover
#' terra::plot(siteLandcover)
#'
#' # NDVI for Eisenwurzen
#' siteNDVI <- get_site_ODS(
#'   deimsid = "https://deims.org/d0a8da18-0881-4ebe-bccf-bc4cb4e25701",
#'   dataset = "ndvi_summer"
#' )
#' siteNDVI
#' terra::plot(siteNDVI)
#' }
#'
#' @section The function output:
#' \figure{get_site_ods_fig.png}{NDVI for Eisenwurzen}
#'
### function get_site_ODS
get_site_ODS <- function(deimsid, dataset = "landcover") {
  # Base URL for OpenDataScience Europe
  ods_url <- "https://s3.eu-central-1.wasabisys.com/eumap/"
  # Define viscurl virtual dataset
  ods_url <- paste0("/vsicurl/", ods_url)
  # Links to individual datasets
  landcover <- paste0("lcv/",
                      "lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019",
                      "_eumap_epsg3035_v0.1.tif")
  clc2018 <- paste0("lcv/",
                    "lcv_landcover.clc_corine_c_100m_0..0cm_2018",
                    "_eumap_epsg3035_v2020.tif")
  osm_buildings <- paste0("lcv/",
                    "lcv_building_copernicus.osm_p_30m_0..0cm_2018..2021",
                    "_eumap_epsg3035_v0.1.tif")
  natura2000 <- paste0("adm/",
                "adm_protected.area_natura2000.osm_p_30m_0..0cm_2019..2021",
                "_eumap_epsg3035_v0.1.tif")
  ndvi_spring <- paste0("lcv/lcv_ndvi_landsat.glad.ard_p50_30m_0..0cm_201903",
              "_eumap_epsg3035_v1.0.tif")
  ndvi_summer <- paste0("lcv/lcv_ndvi_landsat.glad.ard_p50_30m_0..0cm_201906",
                        "_eumap_epsg3035_v1.0.tif")
  ndvi_autumn <- paste0("lcv/",
                        "lcv_ndvi_landsat.glad.ard_p50_30m_0..0cm_201909",
                        "_eumap_epsg3035_v1.0.tif")
  ndvi_winter <- paste0("lcv/lcv_ndvi_landsat.glad.ard_p50_30m_0..0cm_201912",
                        "_eumap_epsg3035_v1.0.tif")
  ndvi_trend <- paste0("lcv/lcv_ndvi_landsat.glad.trend.slope_p50_30m",
                       "_0..0cm_2000..2019_eumap_epsg3035_v1.0.tif")
  forest_broadleaf <- paste0("lcv/lcv_landcover.311_lucas.corine.eml_p_30m",
                             "_0..0cm_2019_eumap_epsg3035_v0.2.tif")
  forest_coniferous <- paste0("lcv/lcv_landcover.312_lucas.corine.eml_p_30m",
                              "_0..0cm_2019_eumap_epsg3035_v0.2.tif")
  forest_mixed <- paste0("/lcv/lcv_landcover.313_lucas.corine.eml_p_30m",
                         "_0..0cm_2019_eumap_epsg3035_v0.2.tif")
  full_url <- dplyr::case_when(
                dataset == "landcover" ~ paste0(ods_url, landcover),
                dataset == "clc2018" ~ paste0(ods_url, clc2018),
                dataset == "osm_buildings" ~ paste0(ods_url, osm_buildings),
                dataset == "natura2000" ~ paste0(ods_url, natura2000),
                dataset == "ndvi_spring" ~ paste0(ods_url, ndvi_spring),
                dataset == "ndvi_summer" ~ paste0(ods_url, ndvi_summer),
                dataset == "ndvi_autumn" ~ paste0(ods_url, ndvi_autumn),
                dataset == "ndvi_winter" ~ paste0(ods_url, ndvi_winter),
                dataset == "ndvi_trend" ~ paste0(ods_url, ndvi_trend),
                dataset == "forest_broadleaf" ~ paste0(ods_url, forest_broadleaf),
                dataset == "forest_coniferous" ~ paste0(ods_url, forest_coniferous),
                dataset == "forest_mixed" ~ paste0(ods_url, forest_mixed),
                TRUE ~ paste("Dataset:", dataset, "unavailable")
              )
  # Make sure `dataset` is among the list of implemented datasets
  if (grepl(pattern = "unavailable", x = full_url, fixed = TRUE)) {
    print(full_url)
    return(NULL)
  }
  # First check that site has a boundary
  boundary <- ReLTER::get_site_info(
    deimsid,
    category = "Boundaries"
  )
  if (is.null(boundary) || !inherits(boundary, "sf")) {
    print("No boundary for requested DEIMS site.")
    return(NULL)
  }
  # terra::rast can address a virtual dataset *without* downloading
  ds <- terra::rast(full_url)
  if (is.null(ds) || !inherits(ds, "SpatRaster")) {
    print("No raster dataset downloaded")
    return(NULL)
  }
  # Crop and mask the raster dataset to the boundary polygon
  # The boundary must be transformed first
  # to the European CRS (EPSG:3035) used by ODS
  boundary <- sf::st_transform(boundary, terra::crs(ds))
  boundary <- terra::vect(boundary)
  ds_site <- terra::mask(terra::crop(ds, boundary), boundary)
  # If this is NDVI, rescale back to (-1.0,1.0) range
  if (length(grep(pattern = "ndvi", x = dataset, fixed = TRUE)) > 0) {
    ds_site <- (ds_site - 100) / 100.0
  }
  return(ds_site)
}