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 infosessioninfo::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
#>
#> ──────────────────────────────────────────────────────────────────────────────
CodePudding user response：
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.