#' function to convert lat long to xyz tile with decimals
#' rounding to tile occurs in \code{get_tilexy}
#' @keywords internal
latlong_to_tilexy <- function(lon_deg, lat_deg, zoom){
# Code assumes lon is 180 to 180, so converts to that
lon_deg <- ifelse(lon_deg > 180, lon_deg - 360, lon_deg)
#Code from https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Coordinates_to_tile_numbers_2
lat_rad <- lat_deg * pi /180
n <- 2.0 ^ zoom
xtile <- (lon_deg + 180.0) / 360.0 * n
ytile <- (1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n
return(c(xtile, ytile))
}
#' function to get a data.frame of all xyz tiles to download
#' @keywords internal
get_tilexy <- function(bbx,z){
#Convert to -180 - +180
bbx[c("xmin","xmax")] <- ifelse( bbx[c("xmin","xmax")] > 180, bbx[c("xmin","xmax")] - 360, bbx[c("xmin","xmax")])
min_tile <- unlist(slippymath::lonlat_to_tilenum(bbx["xmin"],bbx["ymin"],z))
max_tile <- unlist(slippymath::lonlat_to_tilenum(bbx["xmax"],bbx["ymax"],z))
x_all <- seq(from = floor(min_tile[1]), to = floor(max_tile[1]))
y_all <- seq(from = floor(min_tile[2]), to = floor(max_tile[2]))
if(z == 1){
x_all <- x_all[x_all<2]
y_all <- y_all[y_all<2]
} else if(z == 0){
x_all <- x_all[x_all<1]
y_all <- y_all[y_all<1]
}
return(expand.grid(x_all,y_all))
}
#' function to check input type and projection. All input types convert to a
#' SpatialPointsDataFrame for point elevation and bbx for raster.
#' @keywords internal
loc_check <- function(locations, prj = NULL){
if(is.null(nrow(locations))){
nfeature <- length(locations)
} else {
nfeature <- nrow(locations)
}
if(all(class(locations)=="data.frame")){
if(is.null(prj) & !any(class(locations) %in% c("sf", "sfc", "sfg"))){
stop("Please supply a valid sf crs via locations or prj.")
}
locations <- sf::st_as_sf(x = locations, coords = c("x", "y"), crs = prj)
locations$elevation <- vector("numeric", nfeature)
} else if(any(class(locations) %in% c("sf", "sfc", "sfg"))){
sf_crs <- sf::st_crs(locations)
locations$elevation <- vector("numeric", nfeature)
if((is.null(sf_crs) | is.na(sf_crs)) & is.null(prj)){
stop("Please supply an sf object with a valid crs.")
}
} else if(any(class(locations) %in% c("SpatRaster", "SpatVector"))){
sf_crs <- sf::st_crs(locations)
locations <- sf::st_as_sf(terra::as.points(locations, values = FALSE),
coords = terra::crds(locations, df = TRUE),
crs = sf_crs)
locations$elevation <- vector("numeric", nrow(locations))
if((is.null(sf_crs) | is.na(sf_crs)) & is.null(prj)){
stop("Please supply a valid sf crs via locations or prj.")
}
}
#check for long>180
if(is.null(prj)){
prj_test <- sf::st_crs(locations)
} else {
prj_test <- prj
}
lll <- sf::st_is_longlat(prj_test)
if(lll){
if(any(sf::st_coordinates(locations)[,1]>180)){
stop("The elevatr package requires longitude in a range from -180 to 180.")
}
}
locations
}
#' function to project bounding box and if needed expand it
#' @keywords internal
proj_expand <- function(locations,prj,expand){
lll <- sf::st_is_longlat(prj)
#any(grepl("\\bGEOGCRS\\b",sf::st_crs(prj)) |
# grepl("\\bGEODCRS\\b", sf::st_crs(prj)) |
# grepl("\\bGEODETICCRS\\b", sf::st_crs(prj)) |
# grepl("\\bGEOGRAPHICCRS\\b", sf::st_crs(prj)) |
# grepl("\\blonglat\\b", sf::st_crs(prj)) |
# grepl("\\blatlong\\b", sf::st_crs(prj)))
if(is.null(nrow(locations))){
nfeature <- length(locations)
} else {
nfeature <- nrow(locations)
}
if(any(sf::st_bbox(locations)[c("ymin","ymax")] == 0) & lll & is.null(expand)){
# Edge case for lat exactly at the equator - was returning NA
expand <- 0.01
} else if(nfeature == 1 & lll & is.null(expand)){
# Edge case for single point and lat long
expand <- 0.01
} else if(nfeature == 1 & is.null(expand)){
# Edge case for single point and projected
# set to 1000 meters
unit <- sf::st_crs(sf::st_as_sf(locations), parameters = TRUE)$ud_unit
expand <- units::set_units(units::set_units(1000, "m"), unit,
mode = "standard")
expand <- as.numeric(expand)
}
if(!is.null(expand)){
bbx <- sf::st_bbox(locations) + c(-expand, -expand, expand, expand)
} else {
bbx <- sf::st_bbox(locations)
}
bbx <- bbox_to_sf(bbx, prj = prj)
bbx <- sf::st_bbox(sf::st_transform(bbx, crs = ll_geo))
bbx_coord_check <- as.numeric(bbx)
if(any(!bbx_coord_check >= -180 & bbx_coord_check <= 360)){
stop("The elevatr package requires longitude in a range from -180 to 180.")
}
bbx
#sf expand - save for later
#loc_sf <- sf::st_as_sf(locations)
#loc_bbx <- sf::st_bbox(loc_sf)
#bbx_sf <- loc_bbx + c(-expand,-expand,expand,expand)
#names(bbx_sf) <- c("xmin", "ymin", "xmax", "ymax")
#attr(bbx_sf, "class") <- "bbox"
#bbx_sf
}
#' function to clip the DEM
#' @keywords internal
clip_it <- function(rast, loc, expand, clip){
loc_wm <- sf::st_transform(loc, crs = terra::crs(rast))
if(clip == "locations" & !grepl("sfc_POINT", class(sf::st_geometry(loc_wm))[1])){
dem <- terra::mask(terra::crop(rast,loc_wm), loc_wm)
} else if(clip == "bbox" | grepl("sfc_POINT", class(sf::st_geometry(loc_wm))[1])){
bbx <- proj_expand(loc_wm, as.character(terra::crs(rast)), expand)
bbx_sf <- sf::st_transform(bbox_to_sf(bbx), crs = terra::crs(rast))
dem <- terra::mask(terra::crop(rast,bbx_sf), bbx_sf)
}
dem
}
#' Assumes geographic projection
#' sf bbox to poly
#' @param bbox an sf bbox object
#' @param prj defaults to "EPSG:4326"
#' @keywords internal
bbox_to_sf <- function(bbox, prj = 4326) {
sf_bbx <- sf::st_as_sf(sf::st_as_sfc(bbox))
sf_bbx <- sf::st_transform(sf_bbx, crs = prj)
sf_bbx
}
#' Estimate download size of DEMs
#' @param locations the locations
#' @param prj prj string as set earlier by get_elev_point or get_elev_raster
#' @param src the src
#' @param z zoom level if source is aws
#' @keywords internal
estimate_raster_size <- function(locations, prj, src, z = NULL){
locations <- bbox_to_sf(sf::st_bbox(locations),
prj = prj)
locations <- sf::st_transform(locations, crs = 4326)
# Estimated cell size (at equator) from zoom level source
# https://github.com/tilezen/joerd/blob/master/docs/data-sources.md#sources-native-resolution
# Each degree at equator = 111319.9 meters
# Convert ground res to dd
# zoom level 0 = 156543 meters 156543/111319.9
# old resolution (no idea how I calculated these...)
# c(0.54905236, 0.27452618, 0.15455633, 0.07145545, 0.03719130, 0.01901903,
# 0.00962056, 0.00483847, 0.00241219, 0.00120434, 0.00060173, 0.00030075,
# 0.00015035, 0.00007517, 0.00003758)
m_at_equator <- c(156543.0, 78271.5, 39135.8, 19567.9, 9783.9, 4892.0, 2446.0,
1223.0, 611.5, 305.7, 152.9, 76.4, 38.2, 19.1, 9.6, 4.8,
2.4)
z_res <- data.frame(z = 0:16, res_dd = m_at_equator/111319.9)
bits <- switch(src,
aws = 32,
gl3 = 32,
gl1 = 32,
alos = 32,
srtm15plus = 32)
if(src == "aws"){
res <- z_res[z_res$z == z,]$res_dd
} else{
res <- switch(src,
gl3 = 0.0008333,
gl1 = 0.0002778,
alos = 0.0002778,
srtm15plus = 0.004165)
}
num_rows <- (sf::st_bbox(locations)$xmax - sf::st_bbox(locations)$xmin)/res
num_cols <- (sf::st_bbox(locations)$ymax - sf::st_bbox(locations)$ymin)/res
num_megabytes <- (num_rows * num_cols * bits)/8388608
num_megabytes
}
#' OpenTopo Key
#'
#' The OpenTopography API now requires an API Key. This function will grab your
#' key from an .Renviron file
#'
#' @keywords internal
get_opentopo_key <- function(){
if(Sys.getenv("OPENTOPO_KEY")==""){
stop("You have not set your OpenTopography API Key.
Please use elevatr::set_opentopo_key().")
}
Sys.getenv("OPENTOPO_KEY")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.