In ITSLeeds/slopes: Calculate Slopes of Roads, Rivers and Trajectories

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(slopes)

The aim of this vignette is to document issues that have been identified and (hopefully, eventually) fixed. It does not run by default. To run it, change FALSE to TRUE in the code chunk below.

knitr::opts_chunk$set(eval = FALSE)

Error in weighted mean (#12)

This issue showed that the function slope_raster() errors in the final line:

library(sf)
library(raster)
library(slopes)

#import datasets
network_original = st_read("https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true")
network = st_zm(network_original, drop = T) # make sure it has no Z values stored
network
u = "https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true"
download.file(u, "LisboaNASA_clip.tif")
dem = raster("LisboaNASA_clip.tif")
dem

# do they overlap?
raster::plot(dem)
plot(sf::st_geometry(network), add = TRUE)
# why this error?
network$slopeNASA = slope_raster(network, e = dem)
# In addition: Warning messages:
# 1: In max(d, na.rm = TRUE) :
#   no non-missing arguments to max; returning -Inf
# 2: In max(d, na.rm = TRUE) :
#   no non-missing arguments to max; returning -Inf

Does it error on a subset of the network?

network_subset = network[1:9, ]
network$slopesNASA = slope_raster(network_subset, e = dem) # also fails with the same error
dem_subset = crop(dem, network)
network_subset$slopesNASA = slope_raster(network_subset, e = dem_subset) # also fails with the same error
unique(st_geometry_type(network_subset))
unique(st_geometry_type(lisbon_road_segment))
mapview::mapview(network_subset)
network_subset_ls = sf::st_cast(network_subset, to = "LINESTRING")
mapview::mapview(network_subset_ls)
network_subset_ls$slopesNASA = slope_raster(network_subset_ls, e = dem_subset)

Testing on the full network:

# install latest version
remotes::install_github("ropensci/slopes")
library(sf)
library(raster)
library(slopes)

#import datasets
network_original = st_read("https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true")
network = st_zm(network_original, drop = T) # make sure it has no Z values stored
network
u = "https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true"
download.file(u, "LisboaNASA_clip.tif")
dem = raster("LisboaNASA_clip.tif")
dem

# do they overlap?
raster::plot(dem)
plot(sf::st_geometry(network), add = TRUE)
# why this error?
network$slopeNASA = slope_raster(network, e = dem)
network_ls = sf::st_cast(network, "LINESTRING")
network_ls$slopeNASA = slope_raster(network_ls, e = dem)
plot(network_ls["slopeNASA"], lwd = 5)

Issues with terra

remotes::install_github("ropensci/slopes")
library(slopes)
class(lisbon_road_network)
raster::plot(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_network), add = TRUE)
demterra = terra::rast(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_network), add = TRUE)
lisbon_road_network$sloperaster = slope_raster(lisbon_road_network, e = dem_lisbon_raster)
library(terra)
plot(demterra)
lisbon_road_network$slopeterra = slope_raster(lisbon_road_network, e = demterra, terra = T)
summary(lisbon_road_network$sloperaster)
summary(lisbon_road_network$slopeterra)

# Looking at the code, this seems to be the culprit:

# res = as.numeric(terra::extract(e, m[, 1:2], method = method))
# vs
# res = as.numeric(raster::extract(e, m[, 1:2], method = method))
# lets try it...
method = "bilinear"
m = sf::st_coordinates(lisbon_road_segment)[1:2, ]
e = demterra
as.numeric(terra::extract(e, m[, 1:2], method = method))
as.numeric(terra::extract(dem_lisbon_raster, m[, 1:2], method = method))
res_terra = terra::extract(e, m[, 1:2], method = method)
class(res_terra)
res_terra
as.numeric(res_terra[, "r1"])