Making progress on switch to sf/terra.

jhollist · Jun 6, 2023 · 135d824 · 135d824
1 parent 24644a4
commit 135d824
Show file tree

Hide file tree

Showing 11 changed files with 92 additions and 54 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: elevatr
 Title: Access Elevation Data from Various APIs
-Version: 0.4.5.9999
+Version: 1.0.0.9999
 Authors@R: c(person("Jeffrey", "Hollister", email = "hollister.jeff@epa.gov", 
                   role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9254-9740")),
              person("Tarak","Shah", role = "ctb"),
@@ -11,28 +11,28 @@ URL: https://github.com/jhollist/elevatr/
 BugReports: https://github.com/jhollist/elevatr/issues/
 Maintainer: Jeffrey Hollister <hollister.jeff@epa.gov>
 Description: Several web services are available that provide access to elevation
-             data. This package provides access to several of those services and 
-             returns elevation data either as a SpatialPointsDataFrame from 
-             point elevation services or as a raster object from raster 
+             data. This package provides access to many of those services and 
+             returns elevation data either as a simple features POINT/MULTIPOINT from 
+             point elevation services or as a terra SpatRaster object from raster 
              elevation services.  Currently, the package supports access to the
              Amazon Web Services Terrain Tiles <https://registry.opendata.aws/terrain-tiles/>, 
              the Open Topography Global Datasets API <https://opentopography.org/developers/>, 
              and the USGS Elevation Point Query Service <https://apps.nationalmap.gov/epqs/>.
 Depends: R (>= 3.5.0)
 Imports:
-    raster,
     httr,
     jsonlite,
     progressr,
     sf,
+    terra,
     future,
     furrr,
     purrr,
     methods,
     units,
     slippymath,
     curl
-License: CC0
+License: MIT
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.2.3

diff --git a/R/get_elev_point.R b/R/get_elev_point.R
@@ -37,26 +37,31 @@
 #'            increase.  
 #'            Read \url{https://github.com/tilezen/joerd/blob/master/docs/data-sources.md#what-is-the-ground-resolution} 
 #'            for details.  
-#' @return Function returns a \code{SpatialPointsDataFrame} or \code{sf} object 
-#'         in the projection specified by the \code{prj} argument.
+#' @return Function returns an \code{sf} object in the projection specified by 
+#'         the \code{prj} argument.
 #' @export
 #' @examples
 #' \dontrun{
 #' library(elevatr)
 #' library(sf)
-#' library(raster)
+#' library(terra)
 #' 
 #' mts <- data.frame(x = c(-71.3036, -72.8145), 
 #'                   y = c(44.2700, 44.5438), 
 #'                   names = c("Mt. Washington", "Mt. Mansfield"))
 #' ll_prj <- 4326
 #' mts_sf <- st_as_sf(x = mts, coords = c("x", "y"), crs = ll_prj)
-#' mts_raster <- raster(mts_sf, ncol = 2, nrow = 2)
+#' #Empty Raster
+#' mts_raster <- rast(mts_sf, nrow = 5, ncol = 5)
+#' # Raster with cells for each location
+#' mts_raster_loc <- terra::rasterize(x, rast(x, nrow = 10, ncol = 10))
 #' 
 #' get_elev_point(locations = mts, prj = ll_prj)
 #' get_elev_point(locations = mts, units="feet", prj = ll_prj)
 #' get_elev_point(locations = mts_sf)
 #' get_elev_point(locations = mts_raster)
+#' get_elev_point(locations = mts_raster_loc)
+#' 
 #' 
 #' # Code to split into a loop and grab point at a time.
 #' # This is usually faster for points that are spread apart 
@@ -65,9 +70,8 @@
 #' 
 #' elev <- vector("numeric", length = nrow(mts))
 #' for(i in seq_along(mts)){
-#' elev[i]<-suppressMessages(get_elev_point(locations = mts[i,], prj = ll_prj, 
-#'                                         src = "aws", z = 14)$elevation)
-#'                                         }
+#' elev[i]<-get_elev_point(locations = mts[i,], prj = ll_prj, src = "aws", 
+#'                         z = 10)$elevation}
 #' mts_elev <- cbind(mts, elev)
 #' mts_elev
 #' }
@@ -231,13 +235,13 @@ get_epqs <- function(locations, units = c("meters","feet"),
     round(as.numeric(resp$value), 2)
   }
 
-
   coords_df <- split(data.frame(sf::st_coordinates(locations)), 
                      seq_along(locations$elevation))   
 
   #elev_column_name <- make.unique(c(names(locations), "elevation"))
   #elev_column_name <- elev_column_name[!elev_column_name %in% names(locations)]
   elev_column_name <- "elevation"
+  browser()
   message("Downloading point elevations:")
 
   progressr::handlers(

diff --git a/R/get_elev_raster.R b/R/get_elev_raster.R
@@ -49,8 +49,9 @@
 #' @param ... Extra arguments to pass to \code{httr::GET} via a named vector, 
 #'            \code{config}.   See
 #'            \code{\link{get_aws_terrain}} for more details. 
-#' @return Function returns a \code{RasterLayer} in the projection 
-#'         specified by the \code{prj} argument.
+#' @return Function returns a \code{SpatRaster} in the projection 
+#'         specified by the \code{prj} argument or in the projection of the 
+#'         provided locations.
 #' @details Currently, the \code{get_elev_raster} function utilizes the 
 #'          Amazon Web Services 
 #'          (\url{https://registry.opendata.aws/terrain-tiles/}) terrain 
@@ -66,7 +67,7 @@
 #' @examples 
 #' \dontrun{
 #' data(lake)
-#' 
+#' lake_buff  <- st_buffer(lake, 1000)
 #' loc_df <- data.frame(x = runif(6,min=sf::st_bbox(lake)$xmin, 
 #'                                max=sf::st_bbox(lake)$xmax),
 #'                      y = runif(6,min=sf::st_bbox(lake)$ymin, 
@@ -75,6 +76,7 @@
 #' x <- get_elev_raster(locations = loc_df, prj = st_crs(lake) , z=10)
 #' x <- get_elev_raster(lake, z = 12)
 #' x <- get_elev_raster(lake, src = "gl3", expand = 5000)
+#' x <- get_elev_raster(lake_buff, z = 10, clip = "locations")
 #' }
 
 get_elev_raster <- function(locations, z, prj = NULL, 

diff --git a/R/internal.R b/R/internal.R
@@ -46,7 +46,7 @@ loc_check <- function(locations, prj = NULL){
   } 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.")
@@ -64,7 +64,17 @@ loc_check <- function(locations, prj = NULL){
       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(as.points(locations), 
+                              coords = terra::crds(locations, df = TRUE),
+                              crs = st_as_sf)
+    locations$elevation <- vector("numeric", nrow(locations))
+    if(is.null(prj) |is.null(sf_crs) | is.na(sf_crs)){
+      stop("Please supply a valid sf crs via locations or prj.")
+    }
+  }
 
   #check for long>180
   if(is.null(prj)){

diff --git a/man/get_elev_point.Rd b/man/get_elev_point.Rd
diff --git a/man/get_elev_raster.Rd b/man/get_elev_raster.Rd
diff --git a/man/loc_check.Rd b/man/loc_check.Rd
diff --git a/tests/testthat/test-get_elev_point.R b/tests/testthat/test-get_elev_point.R
@@ -1,6 +1,7 @@
 context("get_elev_point")
 library(sf)
 library(elevatr)
+library(terra)
 data("pt_df")
 data("sf_big")
 skip_on_cran()
@@ -12,6 +13,8 @@ sf_sm <- st_as_sf(pt_df, coords = c("x", "y"), crs = ll_prj)
 sf_sm_prj <- st_transform(sf_sm, crs = aea_prj) 
 bad_sf <- st_as_sf(data.frame(x = 1000, y = 1000), coords = c("x", "y"), 
                    crs = ll_prj)
+blank_raster <- rast(sf_sm, nrow = 5, ncol = 5, vals = 1)
+sf_sm_raster <- rasterize(sf_sm, rast(sf_sm, nrow = 10, ncol = 10))
 
 test_that("get_elev_point returns correctly", {
 
@@ -25,13 +28,17 @@ test_that("get_elev_point returns correctly", {
   epqs_sf_aws_z <- get_elev_point(locations = sf_sm, src = "aws", z = 4)
   epqs_sf_aws <- get_elev_point(locations = sf_sm, src = "aws")
   epqs_ft_aws <- get_elev_point(locations = sf_sm, src = "aws", units = "feet")
+  epqs_rast <- get_elev_point(locations = blank_raster)
+  epqs_sf_raster <- get_elev_point(locations = sf_sm_raster)
 
 
 
   # class
   expect_is(epqs_df, "sf")
   expect_is(epqs_sf, "sf")
   expect_is(epqs_sf_prj, "sf")
+  expect_is(epqs_rast, "sf")
+  expect_is(epqs_sf_raster, "sf")
 
   # crs
   expect_equal(st_crs(sf_sm)$wkt,st_crs(epqs_sf)$wkt)
@@ -45,4 +52,7 @@ test_that("get_elev_point returns correctly", {
   expect_equal(epqs_ft_aws$elev_units[1],"feet")
   expect_equal(epqs_sf_aws$elev_units[1],"meters")
 
+  # num points from SpatRaster
+  expect_equal(nrow(epqs_rast), 25)
+  expect_equal(nrow(epqs_sf_raster), 5)
 })
diff --git a/tests/testthat/test-get_elev_raster.R b/tests/testthat/test-get_elev_raster.R
@@ -1,6 +1,7 @@
 context("get_elev_raster")
 library(sf)
 library(elevatr)
+library(terra)
 data("pt_df")
 data("sf_big")
 data("lake")
@@ -13,15 +14,21 @@ sf_sm <- st_as_sf(pt_df, coords = c("x", "y"), crs = ll_prj)
 sf_sm_prj <- st_transform(sf_sm, crs = aea_prj) 
 bad_sf <- st_as_sf(data.frame(x = 1000, y = 1000), coords = c("x", "y"), 
                    crs = ll_prj)
+blank_raster <- rast(sf_sm, nrow = 5, ncol = 5, vals = 1)
+sf_sm_raster <- rasterize(sf_sm, rast(sf_sm, nrow = 10, ncol = 10))
 
 test_that("get_elev_raster returns correctly", {
 
   aws <- get_elev_raster(locations = sf_sm, z = 6, src = "aws")
   aws_prj <- get_elev_raster(locations = sf_sm_prj, z = 6, src = "aws")
+  aws_blnk_raster <- get_elev_raster(locations = blank_raster, z = 6, src = "aws")
+  aws_sf_raster <- get_elev_raster(locations = sf_sm_raster, z = 6, src = "aws")
 
   #class
-  expect_is(aws,"RasterLayer")
-  expect_is(aws_prj,"RasterLayer")
+  expect_is(aws,"SpatRaster")
+  expect_is(aws_prj,"SpatRaster")
+  expect_is(aws_blnk_raster, "SpatRaster")
+  expect_is(aws_sf_raster, "SpatRaster")
 
   #project
   #expect_equal(st_crs(aws)$wkt,st_crs(ll_prj)$wkt)
@@ -34,16 +41,26 @@ test_that("get_elev_raster clip argument works", {
   default_clip <- get_elev_raster(lake, z = 5, clip = "tile")
   bbox_clip <- get_elev_raster(lake, z = 5, clip = "bbox")
   locations_clip <- get_elev_raster(lake, z = 5, clip = "locations")
+  spat_rast_tile <- get_elev_raster(locations = sf_sm_raster, z = 5, 
+                                    src = "aws", clip = "tile")
+  spat_rast_loc <- get_elev_raster(locations = sf_sm_raster, z = 5, 
+                                    src = "aws", clip = "locations")
 
   default_values <- terra::values(default_clip)
   num_cell_default <- length(default_values[!is.na(default_values)])
   bbox_values <- terra::values(bbox_clip)
   num_cell_bbox <- length(bbox_values[!is.na(bbox_values)])
   locations_values <- terra::values(locations_clip)
   num_cell_locations <- length(locations_values[!is.na(locations_values)])
+  default_spat_rast <- terra::values(spat_rast_tile)
+  num_cell_default_spat_rast <- length(default_spat_rast[!is.na(default_spat_rast)])
+  loc_spat_rast <- terra::values(spat_rast_loc)
+  num_cell_loc_spat_rast <- length(loc_spat_rast[!is.na(loc_spat_rast)])
+
 
   expect_true(num_cell_default > num_cell_bbox)
   expect_true(num_cell_bbox > num_cell_locations)
+  expect_true(num_cell_default_spat_rast > num_cell_loc_spat_rast)
 })
 
 test_that("get_elev_raster returns correctly from opentopo", {
@@ -54,8 +71,8 @@ test_that("get_elev_raster returns correctly from opentopo", {
                              clip = "bbox")
 
   #class
-  expect_is(gl1,"RasterLayer")
-  expect_is(gl1_prj,"RasterLayer")
+  expect_is(gl1,"SpatRaster")
+  expect_is(gl1_prj,"SpatRaster")
 
   #project
   #expect_equal(st_crs(gl1)$wkt,st_crs(ll_prj)$wkt)