Merge e45d411 into 90db6ad

DESCRIPTION

@@ -27,7 +27,7 @@ Imports:
 License: CC0
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.0
 Suggests: 
     testthat,
     knitr,

NAMESPACE

@@ -5,3 +5,4 @@ export(get_aws_terrain)
 export(get_elev_point)
 export(get_elev_raster)
 export(get_epqs)
+export(merge_rasters)

R/get_elev_raster.R

@@ -60,36 +60,37 @@
 #' data(lake)
 #' x <- get_elev_raster(lake, z = 12)
 #' }
-#' 
-get_elev_raster <- function(locations, z, prj = NULL,src = c("aws"),
+
+get_elev_raster <- function(locations, z, prj = NULL, src = c("aws"),
                            expand = NULL, clip = c("tile", "bbox", "locations"), 
                            verbose = TRUE, ...){
-  src <- match.arg(src)
+
+  src  <- match.arg(src)
   clip <- match.arg(clip) 
+
   # Check location type and if sp, set prj.  If no prj (for either) then error
   locations <- loc_check(locations,prj)
-  prj <- sp::proj4string(locations)
+  prj       <- sp::proj4string(locations)
 
   # Pass of locations to APIs to get data as raster
   if(src == "aws") {
     raster_elev <- get_aws_terrain(locations, z, prj = prj, 
-                                   expand = expand, ...)
+                                   expand = expand)
   }
-  # Re-project from webmerc back to original and return
+
   if(clip != "tile"){
     message(paste("Clipping DEM to", clip))
     raster_elev <- clip_it(raster_elev, locations, expand, clip)
   }
-  message(paste("Reprojecting DEM to original projection"))
-  raster_elev <- raster::projectRaster(raster_elev, crs = sp::CRS(prj))
+
   if(verbose){
     message(paste("Note: Elevation units are in meters.\nNote: The coordinate reference system is:\n", prj))
   }
+
   raster_elev
+
 }
 
-
-
 #' Get a digital elevation model from the AWS Terrain Tiles
 #' 
 #' This function uses the AWS Terrain Tile service to retrieve an elevation
@@ -118,39 +119,80 @@ get_elev_raster <- function(locations, z, prj = NULL,src = c("aws"),
 #'            vector.              
 #' @export
 #' @keywords internal
+
 get_aws_terrain <- function(locations, z, prj, expand=NULL, ...){
   # Expand (if needed) and re-project bbx to dd
+
   bbx <- proj_expand(sp::bbox(locations),prj,expand)
-  base_url <- "https://s3.amazonaws.com/elevation-tiles-prod/geotiff/"
+
+  base_url <- "https://s3.amazonaws.com/elevation-tiles-prod/geotiff"
+
   tiles <- get_tilexy(bbx,z)
-  dem_list<-vector("list",length = nrow(tiles))
+
+  urls  <-  sprintf("%s/%s/%s/%s.tif", base_url, z, tiles[,1], tiles[,2])
+
+  dem_list <- vector("list",length = nrow(tiles))
+
+  dir <- tempdir()
+
   pb <- progress::progress_bar$new(
     format = "Downloading DEMs [:bar] :percent eta: :eta",
-    total = nrow(tiles), clear = FALSE, width= 60)
-  for(i in seq_along(tiles[,1])){
+    total = length(urls), clear = FALSE, width= 60)
+
+  for(i in seq_along(urls)){
     pb$tick()
     Sys.sleep(1/100)
-    tmpfile <- tempfile()
-    url <- paste0(base_url,z,"/",tiles[i,1],"/",tiles[i,2],".tif")
-    resp <- httr::GET(url,httr::write_disk(tmpfile,overwrite=TRUE), ...)
+    tmpfile <- tempfile(fileext = ".tif")
+
+    resp <- httr::GET(urls[i], httr::write_disk(tmpfile,overwrite=TRUE))
+
     if (!grepl("image/tif", httr::http_type(resp))) {
-      stop(paste("This url:", url,"did not return a tif"), call. = FALSE)
+      stop(paste("This url:", urls[i],"did not return a tif"), call. = FALSE)
     } 
-    dem_list[[i]] <- raster::raster(tmpfile)
-    raster::projection(dem_list[[i]]) <- web_merc
+
+    dem_list[[i]] <- tmpfile
+  }
+
+  merged_elevation_grid <- merge_rasters(dem_list, target_prj = prj)
+
+  return(merged_elevation_grid)
+}
+
+#' Merge Rasters
+#' 
+#' Merge multiple downloaded raster files into a single file. The input `target_prj` 
+#' describes the projection for the new grid.
+#' 
+#' @param raster_list a list of raster file paths to be mosaiced
+#' @param target_prj the target projection of the output raster
+#' @param method the method for resampling/reprojecting. Default is 'bilinear'. 
+#' Options can be found [here](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r)
+#' @param returnRaster if TRUE, return a raster object (default), else, return the file path to the object
+#' @export
+#' @keywords internal
+
+merge_rasters = function(raster_list,  target_prj, method = "bilinear", returnRaster = TRUE){
+
+  message(paste("Mosaicing & Projecting"))
+
+  destfile <- tempfile(fileext = ".tif")
+  files    <- unlist(raster_list)
+
+  if(is.null(target_prj)){
+    r <- raster::raster(files[1])
+    target_prj <- raster::crs(r)
   }
-  origins<-t(data.frame(lapply(dem_list,raster::origin)))
-  min_origin<-c(min(origins[,1]),min(origins[,2]))
-  change_origins <- function(x,y){
-    raster::origin(x)<-y
-    x
+
+  sf::gdal_utils(util = "warp", 
+                 source = files, 
+                 destination = destfile,
+                 options = c("-t_srs", as.character(target_prj),
+                             "-r", method)
+             )
+
+  if(returnRaster){
+    raster::raster(destfile)
+  } else {
+    destfile
   }
-  dem_list <- lapply(dem_list, function(x,y) change_origins(x,min_origin))
-
-  if(length(dem_list) == 1){ 
-    return(dem_list[[1]])
-  } else if (length(dem_list) > 1){
-    message("Merging DEMs")
-    return(do.call(raster::merge, dem_list))
-  } 
 }

R/internal.R

@@ -122,12 +122,12 @@ proj_expand <- function(bbx,prj,expand){
 #' function to clip the DEM
 #' @keywords internal
 clip_it <- function(rast, loc, expand, clip){
-  loc_wm <- sp::spTransform(loc, sp::CRS(web_merc))
+  loc_wm <- sp::spTransform(loc, raster::crs(rast))
   if(clip == "locations" & !grepl("Points", class(loc_wm))){
     dem <- raster::mask(raster::crop(rast,loc_wm), loc_wm)
   } else if(clip == "bbox" | grepl("Points", class(loc_wm))){
-    bbx <- proj_expand(sp::bbox(loc_wm), web_merc, expand)
-    bbx_sp <- sp::spTransform(bbox_to_sp(bbx), sp::CRS(web_merc))
+    bbx <- proj_expand(sp::bbox(loc_wm), as.character(raster::crs(rast)), expand)
+    bbx_sp <- sp::spTransform(bbox_to_sp(bbx), raster::crs(rast))
     dem <- raster::mask(raster::crop(rast,bbx_sp), bbx_sp)
   }
   dem

tests/testthat/test-get_elev_raster.R

@@ -3,10 +3,12 @@ data("pt_df")
 data("sp_big")
 data("lake")
 library(sp)
-ll_prj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
-aea_prj <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
+#ll_prj  <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
+ ll_prj  <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
+#aea_prj <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
+ aea_prj <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
 
-sp_sm <- SpatialPoints(coordinates(pt_df),CRS(ll_prj))
+sp_sm     <- SpatialPoints(coordinates(pt_df),CRS(ll_prj))
 sp_sm_prj <- spTransform(sp_sm,CRS(aea_prj))
 
 test_that("get_elev_raster returns correctly", {