R/topo.R

In AustralianAntarcticDivision/raadtools: Tools for Synoptic Environmental Spatial Data

#' Functions to provide topographic (bathymetry and/or topography) data.
#'
#' Use \code{readtopo} (or its alias \code{readbathy}) to read data
#' from the chosen data set. The function \code{topofile} is used to
#' find the full file name.
#'
#' `readtopo()` and `readbathy()`  accept `terra::ext` or
#' `raster::extent` or just numeric objects for 'xylim'. Alternatively these can be
#' a terra `SpatRaster` or a raster `BasicRaster` (RasterLayer, RasterBrick, or
#' RasterStack) as a template target raster for crop and resize, or reprojection to
#' new raster grid.  The 'resample' argument controls the kind of sampling when
#' regridded or warped. (i.e. if xylim is extent, you get crop(), if it's a grid,
#' it is remodelled to the grid. extent doesn't need crs, grid does).

#' The following data sets are available using the argument \code{topo}.
#' \describe{
#' \item{gebco_08}{The GEBCO_08 Grid, a global 30 arc-second grid largely generated by combining quality-controlled ship depth soundings with interpolation between sounding points guided by satellite-derived gravity data. \url{http://www.gebco.net/data_and_products/gridded_bathymetry_data/}}
#' \item{ibcso}{IBCSO bathymetry data, resolution 1min, use argument \code{polar = TRUE} to return the projected version (polar stereographic with true scale at 71S, WGS84), 500m resolution. \url{http://www.ibcso.org/data.html}. Default is Ice Surface 'ibcso_is', use 'ibcso_bed' for Bedrock.}
#' \item{etopo1}{ETOPO1 is a 1 arc-minute global relief model of Earth's surface that integrates land topography and ocean bathymetry. \url{http://www.ngdc.noaa.gov/mgg/global/global.html}}
#' \item{etopo2}{Historic and deprecated prior version of ETOPO1. \url{http://www.ngdc.noaa.gov/mgg/global/etopo2.html}}
#' \item{kerguelen}{Kerguelen Plateau Bathymetric Grid,  \url{https://data.gov.au/dataset/ds-ga-a05f7893-007f-7506-e044-00144fdd4fa6/details?q=}}
#' \item{george_v_terre_adelie}{A bathymetric Digital Elevation Model (DEM) of the George V and Terre Adelie continental shelf and margin - 100, 250, and 500 metre resolution. \url{http://data.aad.gov.au/aadc/metadata/metadata_redirect.cfm?md=AMD/AU/GVdem_2008}}
#' \item{smith_sandwell}{Global seafloor topography from satellite altimetry and ship depth soundings. \url{http://topex.ucsd.edu/WWW_html/mar_topo.html}}
#' \item{cryosat2}{Antarctica CryoSat-2 Digital Elevation Model (DEM). \url{https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/cryosat}}
#' \item{lake_superior}{Bathymetry of Lake Superior \url{https://www.ngdc.noaa.gov/mgg/greatlakes/superior.html}}
#' \item{ramp}{Radarsat Antarctic digital elevation model V2 \url{https://github.com/AustralianAntarcticDivision/blueant#radarsat-antarctic-digital-elevation-model-v2}}
#' \item{ga_srtm}{Digital Elevation Model (DEM) of Australia with 1 Second Grid}
#' \item{rema_1km, rema_200m, rema_100m, rema_8m}{Reference Elevation Model of Antartica (REMA) for the peninsula, see `read_rema_tiles` for the index}
#' \item{gebco_19}{GEBCO_2019 Grid https://www.gebco.net/data_and_products/gridded_bathymetry_data/gebco_2019/gebco_2019_info.html}
#' }
#' @title Topography data
#' @name readtopo
#' @aliases readtopo topofile readbathy read_rema_tiles
#' @param topo Data source, see Details.
#' @param lon180 Flag for returning data in Atlantic [-180, 180] rather than Pacific [0, 360] view.
#' @param xylim spatial extents or raster grid to crop from source data, see Details
#' @param polar Flag for returning the polar version of the IBCSO data.
#' @param returnfiles Return just the relevant file name
#' @param ... reserved for future use, ignored currently
#' @param resample method to use for GDAL warper resampling (currently via terra), the 'method' argument to 'project()'
#' @return
#' \describe{
#' \item{}{\code{topofile} returns a character string of the full path to a file name}
#' \item{}{\code{readtopo} and \code{readbathy} return the requested data as a RasterLayer (these are aliases)}
#' }
#' @export
#' @importFrom terra crs rast project crop rotate  ext setValues xmin xmax ymin ymax
readtopo <- function(topo = c("gebco_08", "ibcso",
                              "etopo1", "etopo2",
                              "kerguelen",
                              "smith_sandwell", "gebco_14",
                              "macrie1100m", "macrie2100m",
                              "cryosat2",
                              "lake_superior",
                              "ramp", "ibcso_is", "ibcso_bed",

                              "rema_1km",
                              "rema_200m",
                              "rema_100m",
                              "rema_8m",
                              "gebco_19", "gebco_21"),
                     polar = FALSE,
                     lon180 = TRUE,
                     xylim = NULL,
                     returnfiles = FALSE,
                     ..., resample = "bilinear") {
  if (topo[1] == "ga_strm") stop("ga_srtm is not available")
  if (topo[1] == "george_v_terre_adelie") stop("george_v_terre_adelie is not available")
  if (topo[1] == "srtm") stop("srtm is not available")

  topo <- match.arg(topo)
  nowarning <- options(warn = -1)
  on.exit(options(nowarning), add  = TRUE)
  if (!lon180 & topo %in% c("gebco_08", "ibcso", "etopo2")) {
    tfile <- topofile(topo = topo, polar = FALSE, ...)
    if (returnfiles) return(tfile)
    if (is.null(xylim)) res <- .rotate(raster(tfile))
  } else {
    tfile <- topofile(topo = topo, polar = polar, lon180 = lon180, ...)
    if (returnfiles) return(tfile)
    res <- raster(tfile)
    if (is.na(raster::projection(res))) {
      if (raster::couldBeLonLat(res, warnings = FALSE)) {
        raster::projection(res) <- "OGC:CRS84"
      }
    }
    # if (topo[1] == "smith_sandwell") {
    #   raster::projection(res) <- "EPSG:3857"
    # }
  }

  ## sources with wrong extent
  ##if (topo == "etopo1" || topo == "gebco_08") res <- raster::setExtent(res, raster::extent(-180, 180, -90, 90))
  ## sources with missing CRS metadata
  ##llprojs <- c("etopo2", "kerguelen", "george_v_terre_adelie", "lake_superior")
  ##if (topo %in% llprojs)  projection(res) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  ## self-describing, if you always work in lon-lat ...
  ##if (topo == "ibcso_bed") projection(res) <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

  ##if (topo == "cryosat2") projection(res) <- "+proj=stere +lat_0=-90 +lat_ts=-71 +datum=WGS84 +x_0=0 +y_0=0"

  if (!is.null(xylim) && topo == "rema_8m") {
    if (inherits(xylim, "SpatRaster") || inherits(xylim, "BasicRaster")) {
      ## ok keep going
    } else {
      stop("'xylim' is not supported to crop rema_8m, please provide a template raster (as terra or raster)")
    }
  }
  if (!is.null(xylim)) {
    ## terra!
    ## so use might give
    ## c(xmin, xmax, ymin, ymax) or
    ## a terra rast or
    ## a raster raster or
    ## a terra ext or
    ## a raster extent or
    ## and every case is sane, so do the task
    if (inherits(xylim, "SpatExtent")) xylim <- raster::extent(xylim)
    if (inherits(xylim, "numeric") || inherits(xylim, "Extent")) {
      res <- crop(res, xylim)
    }
    if (inherits(xylim, "BasicRaster")) {
      if (is.na(raster::projection(xylim))) stop("cannot use xylim raster template that has no projection/crs")
      tr <- terra::rast(res)
      if (is.na(terra::crs(tr))) stop(sprintf("problem with topo crs: %s", topo[1]))
      res <- raster::raster(terra::project(tr, terra::rast(xylim), method = resample))
    }

    if (inherits(xylim, "SpatRaster")) {
      if (is.na(terra::crs(xylim))) stop("cannot use xylim raster template that has no projection/crs")
      res <- raster::raster(terra::project(terra::rast(res), xylim, method = resample))
    }

  }
  res
}


#' @rdname readtopo
#' @export
readbathy <- readtopo



#' @rdname readtopo
#' @export
#' @importFrom vapour vapour_vrt
topofile <- function(topo = c("gebco_08",  "ibcso",
                              "etopo1", "etopo2",
                              "kerguelen",
                              "smith_sandwell", "gebco_14", "macrie1100m", "macrie2100m", "cryosat2",
                              "lake_superior",
                              "ramp", "ibcso_is", "ibcso_bed",

                              "rema_1km",
                              "rema_200m",
                              "rema_100m",
                              "rema_8m",
                              "gebco_19", "gebco_21"),
                     polar = FALSE,
                     lon180 = TRUE,
                     ...) {
  topo <- match.arg(topo)
  if (topo == "ibcso") topo <- "ibcso_is" ## ??

  if (topo == "rema_8m") {
    r8m_files <- raadfiles::rema_8m_files()
    ##warning(sprintf("rema_8m is a very large **virtual** raster consisting of many (%i) files on disk,\n beware of making subsets that will pull a lot of data into memory", nrow(r8m_files)))
  }
  allf <- allfiles()

  topopath <-  switch(topo,
                      smith_sandwell = if (lon180) raadfiles::smith_sandwell_lon180_files()$fullname else raadfiles::smith_sandwell_files()$fullname,
                      gebco_08 = raadfiles::gebco08_files()$fullname,
                      gebco_14 = raadfiles::gebco14_files()$fullname,
                      gebco_19 = raadfiles::gebco19_files()$fullname,
                      gebco_21 = raadfiles::gebco21_files()$fullname,
                      ibcso_is = dplyr::filter(raadfiles::ibcso_files(all = FALSE))$fullname,
                      ibcso_bed = dplyr::filter(raadfiles::ibcso_bed_files(all = FALSE))$fullname,
                      etopo1 = raadfiles::etopo1_files()$fullname,
                      etopo2 = raadfiles::etopo2_files()$fullname,
                      kerguelen = raadfiles::kerguelen_files()$fullname,

                      macrie1100m = raadfiles::macquarie100m_57S_files()$fullname,
                      macrie2100m = raadfiles::macquarie100m_58S_files()$fullname,
                      cryosat2 = raadfiles::cryosat2_files()$fullname,
                      lake_superior = raadfiles::lakesuperior_files()$fullname,
                      ramp = raadfiles::ramp_files()$fullname,
                      rema_100m = raadfiles::rema_100m_files()$fullname[1L],
                      rema_200m = raadfiles::rema_200m_files()$fullname[1L],
                      rema_1km = raadfiles::rema_1km_files()$fullname[1L],
                      rema_8m =   file.path(dirname(dirname(r8m_files$fullname[1])), "rema_mosaic_8m_dem.vrt"))


  if (topo == "lake_superior") {
    topopath <- vapour::vapour_vrt(topopath, extent = c(-92.20042, -83.99958, 45.99958, 49.50042), projection = "OGC:CRS84")

  }
  ## no longer needed we're using the tif from v2
  # if (topo == "ibcso_bed") {
  #   topopath <- vapour::vapour_vrt(topopath, extent = c(  -3333500, 3334000, -3337000, 3333500),
  #                                  projection = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ")
  # 
  # }
  if (topo %in% c("gebco_08", "etopo1", "etopo2")) {
    topopath <- vapour::vapour_vrt(topopath, extent = c(-180, 180, -90, 90), projection = "OGC:CRS84")
  }
  if (topo == "smith_sandwell") {
##    topopath <- vapour::vapour_vrt(topopath, projection = "EPSG:3857")
  }
  if (topo == "cryosat2") {
    topopath <- vapour::vapour_vrt(topopath,  extent = c(-2821000, 2819000, rev(c(-2421000, 2419000)) ),
                                   projection = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
  }


  if (length(topopath) < 1 || is.na(topopath) || is.null(topopath)) stop(sprintf("cannot find %s", topo))
  topopath[1L]
}


#' @export
#' @name readtopo
read_rema_tiles <- function(...) {
  raster::shapefile(raadfiles::rema_tile_files(all = FALSE)$fullname[1])
}