## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 7,
  fig.ext = "png"
)
eval_chunks <- Sys.getenv("COPERNICUS_DATA_VIGNETTE") != ""

## ----setup, message=FALSE-----------------------------------------------------
library(CopernicusDataspace)
library(dplyr)
library(jsonlite)
library(sf)
library(stars)

## ----data-download, eval=eval_chunks------------------------------------------
# if (dse_has_s3_secret()) {
#   ## Define the area of interest:
#   bbox <- st_bbox(c(xmin =  4.66,
#                     ymin = 52.842,
#                     xmax =  7.15,
#                     ymax = 53.62), crs = 4326)
# 
#   ## Assume that we have already located the asset below within our
#   ## area of interest:
#   id <- "S1A_IW_GRDH_1SDV_20241125T055820_20241125T055845_056707_06F55C_12F9_COG"
# 
#   ## And it's from this collection:
#   cllctn <- "sentinel-1-grd"
# 
#   ## Get the URI for the manifest.safe file. This is the preferred entry
#   ## point into the data. It contains all bands and all meta-information
#   ## (like accurate Ground Control Points):
#   uri_safe <- dse_stac_get_uri(id, "safe_manifest", cllctn)
# 
#   ## Translate the URI into a virtual system indicator (VSI), such that
#   ## the stars package can interact with it:
#   vsi_safe <- dse_s3_uri_to_vsi(uri_safe)
# 
#   ## Set up the S3 credentials such that the stars package (and its
#   ## underpinning GDAL library) will recognise them.
#   dse_s3_set_gdal_options()
# 
#   ## Create a stars proxy object from the VSI. No data is downloaded
#   ## yet. But you can set up lazy data manipulations with tidyverse.
#   proxy <- read_stars(vsi_safe, proxy = TRUE)
# }

## ----meta-info, eval=eval_chunks----------------------------------------------
# ## Use `gdal_info` to retrieve the raster information
# json_info <-
#   gdal_utils("info", vsi_safe, options = c("-json"), quiet = TRUE) |>
#   fromJSON()
# 
# ## Get the ground control points for the specific raster.
# ## You will need it for subsetting the proxy.
# gcp <-
#   st_as_sf(json_info$gcps$gcpList, coords= c("x", "y"), crs =
#              st_crs(json_info$gcps$coordinateSystem$wkt))
# 
# ## A plot showing where the raster's GCPs are, and where the bbox lies
# par(mar = c(1, 1, 1, 1))
# plot(gcp["pixel"] |> st_geometry(), main = "Ground Control Points", col = NA)
# plot(gcp["pixel"], add = TRUE)
# plot(bbox |> st_as_sfc() |> st_geometry(), border = "red", add = TRUE)
# 
# ## Determine which raster indices are within the bbox
# ## of our area of interest:
# gcp_crop <-
#   st_crop(gcp,
#           st_transform(bbox |>
#                          st_as_sfc() |>
#                          st_buffer(units::as_units(1, "km")),
#                        st_crs(json_info$gcps$coordinateSystem$wkt)))
# 
# dim_x <- st_get_dimension_values(proxy, "x") |> sort()
# dim_x <- which(dim_x >= min(gcp_crop$pixel) & dim_x <= max(gcp_crop$pixel))
# dim_y <- st_get_dimension_values(proxy, "y") |> sort()
# dim_y <- which(dim_y >= min(gcp_crop$line) & dim_y <= max(gcp_crop$line))

## ----echo=FALSE---------------------------------------------------------------
knitr::include_graphics("../man/figures/read-data-gcp.png")

## ----data-mask, eval=eval_chunks----------------------------------------------
# masked_proxy <-
#   proxy[,dim_x, dim_y,] |>
#   st_apply(c("x", "y"), \(x) {
#     ## From json_info$bands we know that band 1=VH and 2=VV
#     ## We want VH but only when VV >= 30:
#     ifelse(x[2] >= 30, x[1], NA)
#   })
# 
# ## Store the masked file:
# masked_file <- tempfile(fileext = ".tif")
# write_stars(masked_proxy, masked_file, chunk_size = c(1024, 1024),
#             NA_value = -9999)

## ----data-read-masked, eval=eval_chunks---------------------------------------
# masked <- read_stars(masked_file)
# plot(masked, main = "masked")

## ----echo=FALSE---------------------------------------------------------------
knitr::include_graphics("../man/figures/read-data-masked.png")

## ----data-add-gcp, eval=eval_chunks-------------------------------------------
# ## The masked file has lost original GCP information
# ## So create a reference file where we explicitly add the GCP info
# ref_file <- tempfile(fileext = ".tif")
# 
# ## Update GCP info for selected raster indices:
# gcp <-
#   json_info$gcps$gcpList |>
#   mutate(pixel = pixel - min(dim_x) - 1,
#          line = line - min(dim_y) - 1) |>
#   filter(line >= 0 & pixel >= 0)
# opts <-
#   cbind("-gcp", as.matrix(gcp[, c("pixel", "line", "x", "y", "x")])) |>
#   apply(1, c) |> c()
# 
# ## Create reference file with GCP info
# gdal_utils(
#   "translate",
#   source = masked_file,
#   destination = ref_file,
#   options = c(opts,
#               "-a_srs", json_info$gcps$coordinateSystem$wkt,
#               "-a_nodata", "-9999") ## missing values
# )

## ----data-warp, eval=eval_chunks----------------------------------------------
# ## Warp the masked file to a fixed raster with meaningful CRS
# ## (In this case UTM31 with 10x10m cells)
# warped_file <- tempfile(fileext = ".tif")
# 
# gdal_utils(
#   "warp",
#   source = ref_file,
#   destination = warped_file,
#   options = c(
#     "-t_srs", "EPSG:32631",
#     "-tps",
#     "-tr", c(10, 10), "-tap", "-r", "bilinear",
#     "-srcnodata", "-9999", ## missing values
#     "-dstnodata", "-9999", ## missing values
#     "-overwrite"))
# warped <- read_stars(warped_file)
# 
# plot(warped, reset = FALSE, axes = TRUE, main = "warped")

## ----echo=FALSE---------------------------------------------------------------
knitr::include_graphics("../man/figures/read-data-warped.png")