slopes-eg.md

author	date	output	title
ateucher	2021-06-08	reprex::reprex_document	ash-crane_reprex.R

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(raster)
#> Loading required package: sp
library(slopes)
library(dplyr)

bb <- getbb("Galiano Island, British Columbia, Canada")

gi_roads <- opq(bbox = bb) %>%
  add_osm_feature(key = "name", value = c("Sturdies Bay Ferry Terminal", "Sturdies Bay Road", "Georgeson Bay Road", "Montague Road", "Montague Park Road")) %>%
  osmdata_sf()

route <- gi_roads$osm_lines
exclude <- c(618412545, 618412549, 151400263, 151400266, 151395296, 151395293, 151395294, 151395295, 151400336)
route <- route[!route$osm_id %in% exclude, ]
plot(route["osm_id"])

elevations <- elevations_get(route)
#> Loading required namespace: ceramic
#> Preparing to download: 16 tiles at zoom = 13 from 
#> https://api.mapbox.com/v4/mapbox.terrain-rgb/
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps unknown in Proj4 definition: +proj=merc +lon_0=0
#> +k=1 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum unknown in Proj4 definition
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded ellps unknown in Proj4 definition: +proj=merc +lon_0=0
#> +k=1 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs +type=crs
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum unknown in Proj4 definition
route <- st_transform(route, st_crs(elevations))

route_3d <- slope_3d(route, elevations)

plot_slope(route_3d)

# The line segments don't appear to be in sequential order, so 
# I merged it into a single line, and chopped up the route 
# into smaller (25m) segments (which were then in sequential order):

# Make it into one line:
route2 <- st_line_merge(st_union(route))

# Function to extract segments of a linestring into 
# individual linestrings. From:
# https://dieghernan.github.io/201905_Cast-to-subsegments/
stdh_cast_substring <- function(x, to = "MULTILINESTRING") {
  ggg <- st_geometry(x)
  
  if (!unique(st_geometry_type(ggg)) %in% c("POLYGON", "LINESTRING")) {
    stop("Input should be  LINESTRING or POLYGON")
  }
  for (k in 1:length(st_geometry(ggg))) {
    sub <- ggg[k]
    geom <- lapply(
      1:(length(st_coordinates(sub)[, 1]) - 1),
      function(i) {
        rbind(
          as.numeric(st_coordinates(sub)[i, 1:2]),
          as.numeric(st_coordinates(sub)[i + 1, 1:2])
        )
      }
    ) %>%
      st_multilinestring() %>%
      st_sfc()
    
    if (k == 1) {
      endgeom <- geom
    }
    else {
      endgeom <- rbind(endgeom, geom)
    }
  }
  endgeom <- endgeom %>% st_sfc(crs = st_crs(x))
  if (class(x)[1] == "sf") {
    endgeom <- st_set_geometry(x, endgeom)
  }
  
  if (to == "LINESTRING") {
    endgeom <- endgeom %>% st_cast("LINESTRING")
  }
  return(endgeom)
}

# use st_segmentize to add vertices every 25m then chop up
# into indivdual linestrings using above function:
route2 <- st_segmentize(route2, 25) %>% 
  stdh_cast_substring("LINESTRING") %>%
  st_sf() %>% 
  mutate(id = seq_len(nrow(.)))

plot(elevations)
plot(route2, add = TRUE)

route_3d2 <- slope_3d(route2, elevations)

plot_slope(route_3d2)

^{Created on 2021-06-08 by the reprex package (v2.0.0)}