R/get_elevdata_from_bbox.R

In zappingseb/rayshaderanimate: A rayshader baded animation of GPX files

#' Generate an elevation matrix from boundary box
#' 
#' @param bbox Boundaries box as generated by \link{get_bbox_from_gpx_table}
#' @param type Whether SRTM or EUDEM data should be used for elevation profiles. For EUDEM data
#'   you need to download the right tif file at: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1
#' @param file In case of EUDEM data provide a link to the GeoTIFF file
#' 
#' 
#' @return A matrix with the latitude in columns and the longitude in rows and the elevation
#'   as values (in m)
#' @export
#' @importFrom magrittr %>%
get_elevdata_from_bbox <- function(bbox, type = c("SRTM", "EUDEM"), file = NULL) {
  
  if (is.data.frame(bbox)) {
    bbox <- as.matrix.data.frame(bbox)
  }
  
  if (length(type) > 1) {
    type <- type[1]
  }
  
  # Download raw data from SRTM
  if (type == "EUDEM" && !is.null(file)) {
    data_raw <- .get_raw_elevdata_from_bbox_and_file(bbox, file)
  } else {
    if (type == "EUDEM") {
      warning("No EUDEM file was provided. Using SRTM.")
    }
    data_raw <- .get_raw_elevdata_from_bbox(bbox)
  }
  
  elevation_matrix <- .get_elmatdata_matrix_from_raw(data_raw)
  
  elevation_matrix <- t(elevation_matrix)
  
  if (type == "EUDEM") {
    elevation_labels <- .get_elmat_labels_bbox_size(bbox, dim(elevation_matrix))
  }else {
    elevation_labels <- .get_elmat_labels_from_raw(data_raw)
  }

  elmat_df <- as.data.frame(elevation_matrix)
  elmat_df <- cbind(elmat_df, elevation_labels$deg_elmat_lat) %>% as.data.frame
  elmat_df <- rbind(elmat_df, elevation_labels$deg_elmat_lon) %>% as.data.frame
  
  # Remove the data and make it colnames
  colnames(elmat_df)[1:(dim(elmat_df)[2]-1)] <- elmat_df[
    dim(elmat_df)[1],
    1:(dim(elmat_df)[2]-1)
    ]
  colnames(elmat_df)[length(colnames(elmat_df))] <- "deg_elmat_lat"
  return(elmat_df)
}

#' @noRd
#' @importFrom raster getData
.get_raw_elevdata_from_bbox <- function(bbox) {
  # Get elevation data
  raster_data <- raster::getData("SRTM", lon = bbox[1, 1], lat = bbox[2, 2])
  return(raster::crop(raster_data, raster::extent(bbox)))
}

#' @import sp
#' @importFrom raster raster crop
.get_raw_elevdata_from_bbox_and_file <- function(bbox, file) {
  
  d <- data.frame(lon = bbox[1, ], lat = bbox[2, ])
  coordinates(d) <- c("lon", "lat")
  proj4string(d) <- CRS("+init=epsg:4326") # WGS 84
  # https://spatialreference.org/ref/epsg/3035/proj4/
  d.3035 <- spTransform(d, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs "))
  
  eu_dem_data <- raster::raster(file)
  eu_dem_data_crop <- raster::crop(eu_dem_data, raster::extent(d.3035@bbox))
  
  return(eu_dem_data_crop)
}

#' @noRd
#' @importFrom raster extract extent
.get_elmatdata_matrix_from_raw <- function(srtm_data) {
  return(matrix(
    raster::extract(srtm_data, raster::extent(srtm_data), buffer = 1000),
    nrow = ncol(srtm_data), ncol = nrow(srtm_data)))
}

#' @importFrom raster res
.get_elmat_labels_from_raw <- function(srtm_data) {
  # Create a vector with the long/lat coordinates of the elmat
  return(list(
    deg_elmat_lon = seq(from = srtm_data[[1]]@extent@xmin,
                        to = srtm_data[[1]]@extent@xmax - raster::res(srtm_data)[1],
                        by = raster::res(srtm_data)[1]),
    deg_elmat_lat = seq(from = srtm_data[[1]]@extent@ymax - raster::res(srtm_data)[2],
                        to = srtm_data[[1]]@extent@ymin,
                        by = -raster::res(srtm_data)[2])
  ))
}
.get_elmat_labels_bbox_size <- function(bbox, size) {
  # Create a vector with the long/lat coordinates of the elmat
  return(list(
    deg_elmat_lon = seq(from = bbox[1, 1],
                        to = bbox[1, 2] - (bbox[1, 2] - bbox[1, 1]) / size[2],
                        by = (bbox[1, 2] - bbox[1, 1]) / size[2]),
    deg_elmat_lat = seq(from = bbox[2, 2] - (bbox[2, 2] - bbox[2, 1]) / size[1],
                        to = bbox[2, 1], by = -(bbox[2, 2] - bbox[2, 1]) / size[1])
  ))
}