Merge pull request #1144 from M3nin0/feature/dem-regularization

Support for DEM from MPC

Author: gilbertocamara

NAMESPACE

@@ -126,6 +126,7 @@ S3method(.raster_xres,terra)
 S3method(.raster_ymax,terra)
 S3method(.raster_ymin,terra)
 S3method(.raster_yres,terra)
+S3method(.reg_s2tile_convert,dem_cube)
 S3method(.reg_s2tile_convert,grd_cube)
 S3method(.reg_s2tile_convert,rtc_cube)
 S3method(.slice_dfr,numeric)
@@ -137,6 +138,7 @@ S3method(.source_collection_access_test,usgs_cube)
 S3method(.source_cube,stac_cube)
 S3method(.source_filter_tiles,"cdse_cube_sentinel-1-rtc")
 S3method(.source_filter_tiles,"deafrica_cube_sentinel-1-rtc")
+S3method(.source_filter_tiles,"mpc_cube_cop-dem-glo-30")
 S3method(.source_filter_tiles,"mpc_cube_sentinel-1-grd")
 S3method(.source_filter_tiles,stac_cube)
 S3method(.source_item_get_bands,stac_cube)
@@ -157,6 +159,7 @@ S3method(.source_items_fid,stac_cube)
 S3method(.source_items_new,"aws_cube_landsat-c2-l2")
 S3method(.source_items_new,"deafrica_cube_sentinel-1-rtc")
 S3method(.source_items_new,"deafrica_cube_sentinel-2-l2a")
+S3method(.source_items_new,"mpc_cube_cop-dem-glo-30")
 S3method(.source_items_new,"mpc_cube_landsat-c2-l2")
 S3method(.source_items_new,"mpc_cube_sentinel-1-grd")
 S3method(.source_items_new,"mpc_cube_sentinel-1-rtc")
@@ -172,6 +175,7 @@ S3method(.source_items_tile,"aws_cube_landsat-c2-l2")
 S3method(.source_items_tile,"cdse_cube_sentinel-1-rtc")
 S3method(.source_items_tile,"deafrica_cube_rainfall-chirps-daily")
 S3method(.source_items_tile,"deafrica_cube_rainfall-chirps-monthly")
+S3method(.source_items_tile,"mpc_cube_cop-dem-glo-30")
 S3method(.source_items_tile,"mpc_cube_landsat-c2-l2")
 S3method(.source_items_tile,"mpc_cube_sentinel-1-grd")
 S3method(.source_items_tile,"mpc_cube_sentinel-1-rtc")
@@ -187,6 +191,7 @@ S3method(.source_roi_tiles,"mpc_cube_landsat-c2-l2")
 S3method(.source_roi_tiles,sdc_cube)
 S3method(.source_roi_tiles,stac_cube)
 S3method(.source_tile_get_bbox,"cdse_cube_sentinel-1-rtc")
+S3method(.source_tile_get_bbox,"mpc_cube_cop-dem-glo-30")
 S3method(.source_tile_get_bbox,"mpc_cube_sentinel-1-grd")
 S3method(.source_tile_get_bbox,"mpc_cube_sentinel-1-rtc")
 S3method(.source_tile_get_bbox,stac_cube)
@@ -373,6 +378,7 @@ S3method(sits_reclassify,default)
 S3method(sits_reduce,raster_cube)
 S3method(sits_reduce,sits)
 S3method(sits_regularize,default)
+S3method(sits_regularize,dem_cube)
 S3method(sits_regularize,derived_cube)
 S3method(sits_regularize,raster_cube)
 S3method(sits_regularize,sar_cube)

R/api_cube.R

@@ -350,12 +350,25 @@ NULL
     },
     .default = FALSE
     )
+
+    dem_cube <- .try({
+        .conf("sources", source, "collections", collection, "dem_cube")
+    },
+    .default = FALSE
+    )
+
     if (sar_cube) {
         if (grepl("rtc", col_class, fixed = TRUE))
             unique(c(col_class, "rtc_cube", "sar_cube", s3_class, class(cube)))
         else
             unique(c(col_class, "grd_cube", "sar_cube", s3_class, class(cube)))
-    } else {
+    }
+
+    else if (dem_cube) {
+        unique(c(col_class, "dem_cube", s3_class, class(cube)))
+    }
+
+    else {
         unique(c(col_class, s3_class, class(cube)))
     }
 }

R/api_plot_raster.R

@@ -310,7 +310,7 @@
     # verifies if stars package is installed
     .check_require_packages("stars")
     # verifies if tmap package is installed
-    .check_require_packages("plot")
+    .check_require_packages("tmap")
 
     # deal with color palette
     .check_palette(palette)
@@ -394,7 +394,7 @@
     # verifies if stars package is installed
     .check_require_packages("stars")
     # verifies if tmap package is installed
-    .check_require_packages("plot")
+    .check_require_packages("tmap")
     # precondition - check color palette
     .check_palette(palette)
     # revert the palette

R/api_regularize.R

@@ -236,3 +236,51 @@
     cube_class <- c(cube_class[[1]], "sar_cube", cube_class[-1])
     .cube_set_class(cube, cube_class)
 }
+#' @noRd
+#' @export
+#'
+.reg_s2tile_convert.dem_cube<- function(cube, roi = NULL, tiles = NULL) {
+    # generate Sentinel-2 tiles and intersects it with doi
+    tiles_mgrs <- .s2tile_open(roi, tiles)
+
+    # create a new cube according to Sentinel-2 MGRS
+    cube_class <- .cube_s3class(cube)
+
+    cube <- tiles_mgrs |>
+        dplyr::rowwise() |>
+        dplyr::group_map(~{
+            # prepare a sf object representing the bbox of each image in
+            # file_info
+            cube_crs <- dplyr::filter(cube, .data[["crs"]] == .x[["crs"]])
+            if (nrow(cube_crs) == 0) {
+                cube_crs <- cube
+            }
+            fi_bbox <- .bbox_as_sf(.bbox(
+                x = .fi(cube_crs),
+                default_crs = .crs(cube_crs),
+                by_feature = TRUE
+            ))
+            file_info <- .fi(cube_crs)[.intersects({{fi_bbox}}, .x), ]
+            .cube_create(
+                source = .tile_source(cube_crs),
+                collection = .tile_collection(cube_crs),
+                satellite = .tile_satellite(cube_crs),
+                sensor = .tile_sensor(cube_crs),
+                tile = .x[["tile_id"]],
+                xmin = .xmin(.x),
+                xmax = .xmax(.x),
+                ymin = .ymin(.x),
+                ymax = .ymax(.x),
+                crs = paste0("EPSG:", .x[["epsg"]]),
+                file_info = file_info
+            )
+        }) |>
+        dplyr::bind_rows()
+
+    # Filter non-empty file info
+    cube <- .cube_filter_nonempty(cube)
+
+    # Finalize customizing cube class
+    cube_class <- c(cube_class[[1]], "dem_cube", cube_class[-1])
+    .cube_set_class(cube, cube_class)
+}

R/api_source.R

@@ -65,9 +65,18 @@ NULL
             },
             .default = FALSE
         )
+        # is this a collection of DEM data ?
+        dem_cube <- .try({
+            .conf("sources", source, "collections", collection, "dem_cube")
+            },
+            .default = FALSE
+        )
         # if this is a SAR collection, add "sar_cube" to the class
         if (sar_cube)
             class(source) <- c("sar_cube", class(source))
+        # if this is a DEM collection, add "dem_cube" to the class
+        if (dem_cube)
+            class(source) <- c("dem_cube", class(source))
         # add a class combining source and collection
         class_source_col <- paste(classes[[1]], tolower(collection), sep = "_")
         class(source) <- unique(c(class_source_col, class(source)))

R/api_source_mpc.R

@@ -233,6 +233,35 @@
         collection = collection
     )
 }
+#' @title Get bbox from file info for COP-DEM-GLO-30
+#' @keywords internal
+#' @noRd
+#' @param source     Data source
+#' @param file_info  File info
+#' @param ...        Additional parameters.
+#' @param collection Image collection
+#' @return vector (xmin, ymin, xmax, ymax).
+#' @export
+`.source_tile_get_bbox.mpc_cube_cop-dem-glo-30` <- function(source,
+                                                            file_info, ...,
+                                                            collection = NULL) {
+    .check_set_caller(".source_tile_get_bbox_mpc_dem_30")
+
+    # pre-condition
+    .check_num(nrow(file_info), min = 1)
+
+    # get bbox based on file_info
+    xmin <- min(file_info[["xmin"]])
+    ymin <- min(file_info[["ymin"]])
+    xmax <- max(file_info[["xmax"]])
+    ymax <- max(file_info[["ymax"]])
+
+    # post-condition
+    .check_that(xmin < xmax && ymin < ymax)
+    # create a bbox
+    bbox <- c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax)
+    return(bbox)
+}
 #' @keywords internal
 #' @noRd
 #' @export
@@ -429,6 +458,53 @@
 #' @keywords internal
 #' @noRd
 #' @export
+`.source_items_new.mpc_cube_cop-dem-glo-30` <- function(source,
+                                                        collection,
+                                                        stac_query, ...,
+                                                        tiles = NULL) {
+    .check_set_caller(".source_items_new_mpc_cube_cop-dem-glo-30")
+
+    # COP-DEM-GLO-30 does not support tiles - convert to ROI
+    if (!is.null(tiles)) {
+        roi <- .s2_mgrs_to_roi(tiles)
+        stac_query[["params"]][["intersects"]] <- NULL
+        stac_query[["params"]][["bbox"]] <- c(roi[["lon_min"]],
+                                              roi[["lat_min"]],
+                                              roi[["lon_max"]],
+                                              roi[["lat_max"]]
+        )
+    }
+
+    # Fix temporal interval (All data available in the same date)
+    stac_query[["params"]][["datetime"]] <- "2021-04-21/2021-04-23"
+
+    # Search content
+    items_info <- rstac::post_request(q = stac_query, ...)
+    .check_stac_items(items_info)
+    # fetching all the metadata
+    items_info <- suppressWarnings(
+        rstac::items_fetch(items = items_info, progress = FALSE)
+    )
+
+    # assign href
+    access_key <- Sys.getenv("MPC_TOKEN")
+    if (!nzchar(access_key)) {
+        access_key <- NULL
+    }
+    # Clean old tokens cached in rstac
+    .mpc_clean_token_cache()
+    items_info <- suppressWarnings(
+        rstac::items_sign(
+            items_info, sign_fn = rstac::sign_planetary_computer(
+                httr::add_headers("Ocp-Apim-Subscription-Key" = access_key)
+            )
+        )
+    )
+    return(items_info)
+}
+#' @keywords internal
+#' @noRd
+#' @export
 `.source_items_tile.mpc_cube_sentinel-1-grd` <- function(source,
                                                          items, ...,
                                                          collection = NULL) {
@@ -485,6 +561,20 @@
     })
     rstac::items_reap(items, field = c("properties", "tile"))
 }
+#' @title Organizes items for MPC COP-DEM-GLO-30 collections
+#' @param source     Name of the STAC provider.
+#' @param items      \code{STACItemcollection} object from rstac package.
+#' @param ...        Other parameters to be passed for specific types.
+#' @param collection Collection to be searched in the data source.
+#' @return A list of items.
+#' @keywords internal
+#' @noRd
+#' @export
+`.source_items_tile.mpc_cube_cop-dem-glo-30` <- function(source,
+                                                        items, ...,
+                                                        collection = NULL) {
+    rep("NoTilingSystem", rstac::items_length(items))
+}
 #' @title Filter S1 GRD tiles
 #' @noRd
 #' @param source  Data source
@@ -509,6 +599,19 @@
         tiles = tiles)
 
 }
+#' @title Filter COP-DEM-GLO-30 tiles
+#' @noRd
+#' @param source  Data source
+#' @param cube    Cube to be filtered
+#' @param tiles   Tiles to be selected
+#' @return Filtered cube
+#' @export
+`.source_filter_tiles.mpc_cube_cop-dem-glo-30` <- function(source,
+                                                           collection,
+                                                           cube,
+                                                           tiles) {
+    return(cube)
+}
 #' @title Check if roi or tiles are provided
 #' @param source        Data source
 #' @param roi           Region of interest

R/sits_merge.R

@@ -207,16 +207,13 @@ sits_merge.raster_cube <- function(data1, data2, ...) {
     # Get data1 timeline.
     d1_tl <- unique(as.Date(.cube_timeline(data1)[[1]]))
     # Create new `file_info` using dates from `data1` timeline.
-    fi_new <- purrr::map(seq_len(nrow(data2)), function(row) {
-        data_row <- data2[row,]
-
-        fi <- .fi(data_row)
-        fi[["date"]] <- as.Date(d1_tl[1:nrow(data_row)])
-
-        return(fi)
+    fi_new <- purrr::map(sits_timeline(data1), function(date_row) {
+        fi <- .fi(data2)
+        fi[["date"]] <- as.Date(date_row)
+        fi
     })
     # Assign the new `file_into` into `data2`
-    data2[["file_info"]] <- fi_new
+    data2[["file_info"]] <- list(dplyr::bind_rows(fi_new))
     # Merge cubes and return
     .cube_merge(data1, data2)
 }

R/sits_regularize.R

@@ -205,6 +205,52 @@ sits_regularize.sar_cube <- function(cube, ...,
     )
     return(cube)
 }
+
+#' @rdname sits_regularize
+#' @export
+sits_regularize.dem_cube <- function(cube, ...,
+                                     res,
+                                     output_dir,
+                                     roi = NULL,
+                                     tiles = NULL,
+                                     multicores = 2L,
+                                     progress = TRUE) {
+    # Preconditions
+    .check_raster_cube_files(cube)
+    .check_num_parameter(res, exclusive_min = 0)
+    output_dir <- .file_normalize(output_dir)
+    .check_output_dir(output_dir)
+    .check_num_parameter(multicores, min = 1, max = 2048)
+    .check_progress(progress)
+    # Check for ROI and tiles
+    .check_roi_tiles(roi, tiles)
+    # Display warning message in case STAC cube
+    # Prepare parallel processing
+    .parallel_start(workers = multicores)
+    on.exit(.parallel_stop(), add = TRUE)
+    # Convert input sentinel1 cube to sentinel2 grid
+    cube <- .reg_s2tile_convert(cube = cube, roi = roi, tiles = tiles)
+    .check_that(nrow(cube) > 0,
+                msg = .conf("messages", "sits_regularize_roi")
+    )
+    # Filter tiles
+    if (is.character(tiles)) {
+        cube <- .cube_filter_tiles(cube, tiles)
+    }
+    # DEMs don't have the temporal dimension, so the period is fixed in 1 day.
+    period <- "P1D"
+    # Call regularize in parallel
+    cube <- .reg_cube(
+        cube = cube,
+        res = res,
+        roi = roi,
+        period = period,
+        output_dir = output_dir,
+        progress = progress
+    )
+    return(cube)
+}
+
 #' @rdname sits_regularize
 #' @export
 sits_regularize.derived_cube <- function(cube, ...) {

inst/extdata/config_messages.yml

@@ -278,6 +278,7 @@
 .source_url: "invalid URL for requested source provider"
 .source_tile_get_bbox: "unable to retrieve images given a bounding box"
 .source_tile_get_bbox_cdse_s1_rtc: "unable to retrieve file information for S1 RTC cubes"
+.source_tile_get_bbox_mpc_dem_30: "unable to retrieve file information for COP DEM GLO 30m"
 .source_tile_get_bbox_mpc_s1_grd: "unable to retrieve file information for S1 GRD cubes"
 .stac_format_platform: "platform name should be unique (e.g. Landsat-8 or Sentinel-1A)"
 .stac_select_bands: "some bands for this product are not pre-configured in sits\n please include them in you user configuration file."