4 Digital Elevation Model (DEM) Mosaic

Here we create the best quality available Digital Elevation Model (DEM). Derived from the last return of the LiDAR survey, carried out in 2016, Copernicus GLO30 (Neteler, Haas, and Metz 2022) and NASADEM (Crippen et al. 2016).

The first step is to correct the elevation of the digital terrain model, derived from the LiDAR survey. The reported elevations form the survey are relative to the reference ellipsoid. Therefore these were corrected to orthometric elevations using the Earth Gravitational Model 2008 (EGM2008) geoid. This step is required because the LiDAR DTM does not have complete coverage of the site and we therefore need to merge it with global DTM products, the heights of which are reported relative to the geoid height.

Three DEM sources were used, prioritising the LiDAR data, followed by COP GLO30, and, where neither of these sources contained data, NASADEM was used. These sources were mosaicked using gdalwarp via the {ezwarp} R package (Graham 2024), with bilinear resampling and resolution of 30 m.

The finalised DEM is presented in the following Figure 4.1.

Figure 4.1: Composite DEM comprising sources from ALS, Copernicus GLO30 and NASADEM.

4.1 Code

Merge DEM Code R/mergeDEM_Sabah.R

#' Merge DEM sources
#'
#' Create a continuous surface elevation map for Kuamut.
#'
#' @param kuamut_bounds character file path for the Kuamut area boundary
#' @param Kuamut_src character file path for the Kuamut LiDAR DTM
#' @return a list with the path to the merged DEM
#'
#' @details
#' Download two GLO30 and NASADEM DEM sources and merge them with the Kuamut
#' LiDAR DTM to create a continuous surface elevation map of the Kuamut area.
merge_dem_srcs <- function(kuamut_bounds, Kuamut_src) {
  kuamut_terr <- rast(Kuamut_src)

  glo30_dem_src <- demsrc::mpc_dtm_src(kuamut_terr)
  nasa_dem_src <- demsrc::mpc_dtm_src(kuamut_terr, collection = "nasadem")
  # save file
  if (!dir.exists(set_path("RoutputProducts"))) {
    dir.create(set_path("RoutputProducts"))
  }

  out_ras_path <- set_path("RoutputProducts/dem30_filled_Sabah.tif")

  kb <- read_sf(kuamut_bounds) |>
    st_make_valid()

  warp_dtm <- function(.res, .out) {
    ezwarp(
      x = c(
        as.list(nasa_dem_src),
        as.list(glo30_dem_src),
        kuamut_terr
      ),
      y = kb,
      cutline = kb,
      res = .res,
      nodata = -999,
      filename = .out,
      engine = "sf"
    )
  }
  dtm1 <- warp_dtm(30, out_ras_path)


  return(out_ras_path)
}

Convert from ellipsoidal to orthometric datum R/geoid_tools.R

#' Download the geoid model gtx file
#' @param geoid character the geoid model to download one of: "egm_08_25"
#' or "egm_96_15".
#' @param force logical if TRUE, force download even if file exists
#' @return character file path for the downloaded gtx file
#' @details This function downloads the geoid model gtx file from the OSGeo
#' download server. The geoid model is used to correct the LiDAR DTM elevations
#' to be relative to the geoid.
download_gtx <- function(geoid = c("egm_08_25", "egm_96_15"), force = FALSE) {
  geoid <- rlang::arg_match(geoid)

  geoid_url <- switch(geoid,
    egm_08_25 = "https://download.osgeo.org/proj/vdatum/egm08_25/egm08_25.gtx",
    egm_96_15 = "https://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx"
  )

  if (!dir.exists(geoid)) {
    dir.create(geoid)
  }

  geoid_file <- file.path(geoid, paste0(geoid, ".gtx"))

  # if force is true and file exists
  if (!force && file.exists(geoid_file)) {
    cli::cli_alert_success("Geoid model gtx file already downloaded!")
    return(geoid_file)
  }

  cli::cli_alert_info("Downloading {geoid} model gtx file...")
  curl::curl_download(geoid_url, geoid_file, quiet = FALSE)
  return(geoid_file)
}


#' Make the Kuamut LiDAR DTM elevations relative to the geoid
#' @param Kuamut_src character file path for the kuamut liDAR DTM
#' @param destfile character file path for the corrected DTM
#' @return character file path for the corrected DTM
#' @details This function calls `terr_geoid` to get the geoid model for the
#' Kuamut area and then corrects the LiDAR DTM elevations, which are relative to
#' the ellipsoid, to be relative to the geoid. This function uses the `gdalwarp`
#' utility to perform the geoid correction with the appropriate gtx file.
geoid_correct_dtm_gdal <- function(
    Kuamut_src,
    geoid = c("egm_08_25", "egm_96_15"),
    destfile = set_path(
      "lidar_dtm_geoid_corr_gdal.tif",
      "data_out/RoutputProducts"
    )) {
  geoid <- rlang::arg_match(geoid)
  geoid_file <- download_gtx(geoid)
  # Define the target spatial reference system with the geoid grid
  t_srs <- paste(
    terra::crs(terra::rast(Kuamut_src), proj = TRUE),
    paste0("+geoidgrids=", geoid_file)
  )

  # Use gdal_utils to perform the warp operation
  gdal_utils(
    util = "warp",
    source = Kuamut_src,
    destination = destfile,
    options = c("-t_srs", t_srs),
    quiet = FALSE
  )
  return(destfile)
}