-
-
Save Pakillo/c0bd8f6e96e87625e715d3870522653f to your computer and use it in GitHub Desktop.
| library(sf) | |
| library(terra) | |
| ## Define area | |
| coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3), | |
| y = c(36.7, 36.8, 36.7, 36.8)) | |
| coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326) | |
| ## Download elevation data | |
| elev.ras <- elevatr::get_elev_raster(coords.sf, z = 13) | |
| elev <- rast(elev.ras) | |
| ## Calculate hillshade | |
| slopes <- terrain(elev, "slope", unit = "radians") | |
| aspect <- terrain(elev, "aspect", unit = "radians") | |
| hs <- shade(slopes, aspect) | |
| ## Plot hillshading as basemap | |
| # (here using terra::plot, but could use tmap) | |
| plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE, | |
| xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82)) | |
| # overlay with elevation | |
| plot(elev, col = terrain.colors(25), alpha = 0.5, legend = FALSE, | |
| axes = FALSE, add = TRUE) | |
| # add contour lines | |
| contour(elev, col = "grey40", add = TRUE) | |
| #### Overlaying orthoimage #### | |
| ## Could use own image, but here downloading w/ maptiles | |
| library(maptiles) | |
| ortho <- get_tiles(ext(-5.50, -5.30, 36.7, 36.8), | |
| provider = "Esri.WorldImagery", zoom = 13) | |
| ## Plot | |
| plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE, | |
| xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82)) | |
| # overlay with elevation | |
| plot(ortho, alpha = 150, add = TRUE) | |
| #### 3-D map with rayshader and rayvista #### | |
| library(rayshader) | |
| library(rayvista) | |
| graz.3D <- plot_3d_vista(dem = elev) | |
| render_snapshot(filename = "rayvista.png") |
Alright. I can't reproduce the problem in my computer so I'm afraid I can't try proposed solutions.
But I'm guessing the below should work? Here using sp rather than sf, because {elevatr} still uses the former.
library(sp)
library(terra)
## Define area
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
y = c(36.7, 36.8, 36.7, 36.8))
coordinates(coords) <- c("x", "y")
proj4string(coords) <- CRS("EPSG:4326")
## Download elevation data
elev.ras <- elevatr::get_elev_raster(coords, z = 13)
elev <- rast(elev.ras)An alternative would be to download elevation data using geodata::elevation_3s(), see https://github.com/rspatial/geodata
@Pakillo oddly, I still get an error. Not sure what is going on? Any advice for troubleshooting?
After running the proj4string() line, I get:
Error in CRS("EPSG:4326") : NA
The problem is with CRS handling in elevatr/sp (see USEPA/elevatr#48). I'm afraid I don't know how to fix it easily.
You can use {geodata} to get elevation data, as an alternative to {elevatr}:
library(sf)
library(terra)
#### Define area ####
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
y = c(36.7, 36.8, 36.7, 36.8))
coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)
#### Download elevation data ####
## Using elevatr
# elev.ras <- elevatr::get_elev_raster(coords, z = 10)
# elev <- rast(elev.ras)
## Using geodata
elev <- geodata::elevation_3s(lon = -5.4, lat = 36.75, path = getwd())
elev <- crop(elev, coords.sf)
#### Calculate hillshade ####
slopes <- terrain(elev, "slope", unit = "radians")
aspect <- terrain(elev, "aspect", unit = "radians")
hs <- shade(slopes, aspect)
#### Make map ####
## Plot hillshading as basemap
# (here using terra::plot, but could use tmap)
plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE,
xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82))
# overlay with elevation
plot(elev, col = terrain.colors(25), alpha = 0.5, legend = FALSE,
axes = FALSE, add = TRUE)
# add contour lines
contour(elev, col = "grey40", add = TRUE)Hope it helps
I'm getting the same CRS/PROJ errors with elevatr as @BrentPease1 describes, but the geodata alternative worked.
Thanks for sharing this!
@Pakillo Sorry, I meant to mention that I also tried that yesterday with the same result:
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': NA