r elevatr AWS and USGS elevation are completely different

Question

I am using the elevatr package in R and realized that some of elevation values I get when calling get_aws_points() are completely different than when calling get_elev_point(). In some cases they are thousands of feet apart. I have have checked the docs from AWS here: https://github.com/tilezen/joerd/blob/master/docs/data-sources.md and see that California should be using USGS 3DEP (either 3 or 10 meter resolution in California). The USGS point query service (https://nationalmap.gov/epqs/) uses a 1/3 arc-second dem (so 10m). These should essentially be the same dataset. I have creatted the following reproducible example here. Does anyone have a clue what is going on here? None of the returned elevations are identical. Both functions call elevation in meters and the z argument in get_aws_points() seems to have no effect on the returned elevation value.

library(tidyverse)
library(elevatr)
library(USAboundaries)
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

# Compare get_aws_points and get_elev_points

# Import California Polygon and project to California Albers
ca.poly <- us_states()%>%
  filter(state_name == "California")%>%
  st_transform(3310)

# Create Sample Points
pts <- ca.poly%>%
  st_make_grid(10000)%>%
  st_centroid()%>%
  st_sf()%>%
  filter(st_intersects(., ca.poly, sparse = FALSE))

# Get elevation from AWS
e1 <- as.data.frame(get_aws_points(pts), z = 1)
#> Warning: PROJ support is provided by the sf and terra packages among others

#> Warning: PROJ support is provided by the sf and terra packages among others
#> Mosaicing & Projecting
#> Note: Elevation units are in meters.
pts$AWS_Meters <- e1$elevation

# Get Elevation with higher zoom from AWS
e2 <- as.data.frame(get_aws_points(pts), z = 14)
#> Warning: PROJ support is provided by the sf and terra packages among others

#> Warning: PROJ support is provided by the sf and terra packages among others
#> Mosaicing & Projecting
#> Note: Elevation units are in meters.
pts$AWS_Meters_Z14 <- e2$elevation

# Get elevation from USGS
e3 <- as.data.frame(get_elev_point(pts))
#> Warning: PROJ support is provided by the sf and terra packages among others
#> Downloading point elevations:
#> Note: Elevation units are in meters
pts$USGS_Meters <- e3$elevation

# Calculate difference
pts$elevDif <- pts$AWS_Meters - pts$USGS_Meters

summary(pts$elevDif)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -599.0500  -26.4975    0.6650   -0.2794   31.6025  486.2700

# Plot difference as map
ggplot(pts)+
  geom_sf(aes(color = elevDif), size = 1)+
  scale_colour_gradientn(colours=rainbow(6))+
  labs(title = "AWS Elevation - USGS Elevation [meters]",
       color = "Meters")

# Plot difference as points
ggplot(pts)+
  geom_abline(slope = 1, intercept = 0, color = "red", size = 1)+
  geom_point(aes(x = AWS_Meters, y = USGS_Meters), alpha = 0.2)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.

^{Created on 2022-11-22 with reprex v2.0.2}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.2 (2022-10-31 ucrt)
#>  os       Windows 10 x64 (build 19042)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language (EN)
#>  collate  English_United States.utf8
#>  ctype    English_United States.utf8
#>  tz       America/New_York
#>  date     2022-11-22
#>  pandoc   2.19.2 @ C:/Program Files/RStudio/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package           * version date (UTC) lib source
#>  assertthat          0.2.1   2019-03-21 [1] CRAN (R 4.2.2)
#>  backports           1.4.1   2021-12-13 [1] CRAN (R 4.2.0)
#>  broom               1.0.1   2022-08-29 [1] CRAN (R 4.2.2)
#>  cellranger          1.1.0   2016-07-27 [1] CRAN (R 4.2.2)
#>  class               7.3-20  2022-01-16 [2] CRAN (R 4.2.2)
#>  classInt            0.4-8   2022-09-29 [1] CRAN (R 4.2.2)
#>  cli                 3.4.1   2022-09-23 [1] CRAN (R 4.2.2)
#>  codetools           0.2-18  2020-11-04 [2] CRAN (R 4.2.2)
#>  colorspace          2.0-3   2022-02-21 [1] CRAN (R 4.2.2)
#>  crayon              1.5.2   2022-09-29 [1] CRAN (R 4.2.2)
#>  curl                4.3.3   2022-10-06 [1] CRAN (R 4.2.2)
#>  DBI                 1.1.3   2022-06-18 [1] CRAN (R 4.2.2)
#>  dbplyr              2.2.1   2022-06-27 [1] CRAN (R 4.2.2)
#>  digest              0.6.30  2022-10-18 [1] CRAN (R 4.2.2)
#>  dplyr             * 1.0.10  2022-09-01 [1] CRAN (R 4.2.2)
#>  e1071               1.7-12  2022-10-24 [1] CRAN (R 4.2.2)
#>  elevatr           * 0.4.2   2022-01-07 [1] CRAN (R 4.2.2)
#>  ellipsis            0.3.2   2021-04-29 [1] CRAN (R 4.2.2)
#>  evaluate            0.18    2022-11-07 [1] CRAN (R 4.2.2)
#>  fansi               1.0.3   2022-03-24 [1] CRAN (R 4.2.2)
#>  farver              2.1.1   2022-07-06 [1] CRAN (R 4.2.2)
#>  fastmap             1.1.0   2021-01-25 [1] CRAN (R 4.2.2)
#>  forcats           * 0.5.2   2022-08-19 [1] CRAN (R 4.2.2)
#>  fs                  1.5.2   2021-12-08 [1] CRAN (R 4.2.2)
#>  furrr               0.3.1   2022-08-15 [1] CRAN (R 4.2.2)
#>  future              1.29.0  2022-11-06 [1] CRAN (R 4.2.2)
#>  gargle              1.2.1   2022-09-08 [1] CRAN (R 4.2.2)
#>  generics            0.1.3   2022-07-05 [1] CRAN (R 4.2.2)
#>  ggplot2           * 3.4.0   2022-11-04 [1] CRAN (R 4.2.2)
#>  globals             0.16.1  2022-08-28 [1] CRAN (R 4.2.1)
#>  glue                1.6.2   2022-02-24 [1] CRAN (R 4.2.2)
#>  googledrive         2.0.0   2021-07-08 [1] CRAN (R 4.2.2)
#>  googlesheets4       1.0.1   2022-08-13 [1] CRAN (R 4.2.2)
#>  gtable              0.3.1   2022-09-01 [1] CRAN (R 4.2.2)
#>  haven               2.5.1   2022-08-22 [1] CRAN (R 4.2.2)
#>  highr               0.9     2021-04-16 [1] CRAN (R 4.2.2)
#>  hms                 1.1.2   2022-08-19 [1] CRAN (R 4.2.2)
#>  htmltools           0.5.3   2022-07-18 [1] CRAN (R 4.2.2)
#>  httr                1.4.4   2022-08-17 [1] CRAN (R 4.2.2)
#>  jsonlite            1.8.3   2022-10-21 [1] CRAN (R 4.2.2)
#>  KernSmooth          2.23-20 2021-05-03 [2] CRAN (R 4.2.2)
#>  knitr               1.40    2022-08-24 [1] CRAN (R 4.2.2)
#>  labeling            0.4.2   2020-10-20 [1] CRAN (R 4.2.0)
#>  lattice             0.20-45 2021-09-22 [2] CRAN (R 4.2.2)
#>  lifecycle           1.0.3   2022-10-07 [1] CRAN (R 4.2.2)
#>  listenv             0.8.0   2019-12-05 [1] CRAN (R 4.2.2)
#>  lubridate           1.9.0   2022-11-06 [1] CRAN (R 4.2.2)
#>  magrittr            2.0.3   2022-03-30 [1] CRAN (R 4.2.2)
#>  mime                0.12    2021-09-28 [1] CRAN (R 4.2.0)
#>  modelr              0.1.9   2022-08-19 [1] CRAN (R 4.2.2)
#>  munsell             0.5.0   2018-06-12 [1] CRAN (R 4.2.2)
#>  parallelly          1.32.1  2022-07-21 [1] CRAN (R 4.2.1)
#>  pillar              1.8.1   2022-08-19 [1] CRAN (R 4.2.2)
#>  pkgconfig           2.0.3   2019-09-22 [1] CRAN (R 4.2.2)
#>  prettyunits         1.1.1   2020-01-24 [1] CRAN (R 4.2.2)
#>  progress            1.2.2   2019-05-16 [1] CRAN (R 4.2.2)
#>  progressr           0.11.0  2022-09-02 [1] CRAN (R 4.2.2)
#>  proxy               0.4-27  2022-06-09 [1] CRAN (R 4.2.2)
#>  purrr             * 0.3.5   2022-10-06 [1] CRAN (R 4.2.2)
#>  R6                  2.5.1   2021-08-19 [1] CRAN (R 4.2.2)
#>  raster              3.6-3   2022-09-18 [1] CRAN (R 4.2.2)
#>  Rcpp                1.0.9   2022-07-08 [1] CRAN (R 4.2.2)
#>  readr             * 2.1.3   2022-10-01 [1] CRAN (R 4.2.2)
#>  readxl              1.4.1   2022-08-17 [1] CRAN (R 4.2.2)
#>  reprex              2.0.2   2022-08-17 [1] CRAN (R 4.2.2)
#>  rgdal               1.6-2   2022-11-09 [1] CRAN (R 4.2.2)
#>  rlang               1.0.6   2022-09-24 [1] CRAN (R 4.2.2)
#>  rmarkdown           2.18    2022-11-09 [1] CRAN (R 4.2.2)
#>  rstudioapi          0.14    2022-08-22 [1] CRAN (R 4.2.2)
#>  rvest               1.0.3   2022-08-19 [1] CRAN (R 4.2.2)
#>  scales              1.2.1   2022-08-20 [1] CRAN (R 4.2.2)
#>  sessioninfo         1.2.2   2021-12-06 [1] CRAN (R 4.2.2)
#>  sf                * 1.0-9   2022-11-08 [1] CRAN (R 4.2.2)
#>  slippymath          0.3.1   2019-06-28 [1] CRAN (R 4.2.2)
#>  sp                  1.5-1   2022-11-07 [1] CRAN (R 4.2.2)
#>  stringi             1.7.8   2022-07-11 [1] CRAN (R 4.2.1)
#>  stringr           * 1.4.1   2022-08-20 [1] CRAN (R 4.2.2)
#>  terra               1.6-17  2022-09-10 [1] CRAN (R 4.2.1)
#>  tibble            * 3.1.8   2022-07-22 [1] CRAN (R 4.2.2)
#>  tidyr             * 1.2.1   2022-09-08 [1] CRAN (R 4.2.2)
#>  tidyselect          1.2.0   2022-10-10 [1] CRAN (R 4.2.2)
#>  tidyverse         * 1.3.2   2022-07-18 [1] CRAN (R 4.2.2)
#>  timechange          0.1.1   2022-11-04 [1] CRAN (R 4.2.2)
#>  tzdb                0.3.0   2022-03-28 [1] CRAN (R 4.2.2)
#>  units               0.8-0   2022-02-05 [1] CRAN (R 4.2.2)
#>  USAboundaries     * 0.4.0   2021-10-12 [1] CRAN (R 4.2.2)
#>  USAboundariesData   0.4.0   2022-11-23 [1] https://ropensci.r-universe.dev (R 4.2.2)
#>  utf8                1.2.2   2021-07-24 [1] CRAN (R 4.2.2)
#>  vctrs               0.5.0   2022-10-22 [1] CRAN (R 4.2.2)
#>  withr               2.5.0   2022-03-03 [1] CRAN (R 4.2.2)
#>  xfun                0.34    2022-10-18 [1] CRAN (R 4.2.2)
#>  xml2                1.3.3   2021-11-30 [1] CRAN (R 4.2.2)
#>  yaml                2.3.6   2022-10-18 [1] CRAN (R 4.2.1)
#> 
#>  [1] C:/Users/AMURRA02/RPackages
#>  [2] C:/Program Files/R/R-4.2.2/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Answer 1

Thanks for the question and detailed reprex.

The problem is in the two lines:

e1 <- as.data.frame(get_aws_points(pts), z = 1)
e2 <- as.data.frame(get_aws_points(pts), z = 14)

The z argument is an argument for get_aws_points. You are passing it as an argument for as.data.frame. As such, your elevations on e1 and e2 are both using the default z of 5 which returns a large pixel (several km). I wouldn't expect the USGS and the AWS with z of 5 elevations to very closely match.

If you change those two lines to:

e1 <- as.data.frame(get_aws_points(pts, z = 1))
e2 <- as.data.frame(get_aws_points(pts, z = 14))

Then you will get what you want. The e1 elevations will be very different than USGS, but the e2 elevations should be very close. Big caveat here... z=14 will return a raster for approximately all of California with an approximate pixel size of 7 meters. This will be large! I would scale this back. I tested it with a z = 10 and the difference between the two sources ranged from -81 to 52 meters. This is expected as the terrain tiles scale based on the zoom so the elevation at a given zoom level will be resampled. The USGS point service is pulling directly from a the 10m NED and is not resampled.

Also, since your point elevations are all in the US, I would use the USGS point elevation service for these data and not mess with the AWS source. The AWS point elevations are, admittedly, a bit of a hack to have a point elevations for global applications. There is no global analogue to the USGS service that I know of.