add google elevation source?

[temospena]

Maybe to add as an alternative to ceramic/mapbox? It has a very good resolution depending on the place on earth.

Using googleway to retrieve google elevation:

library(tidyverse)
library(sf)
library(slopes)
library(googleway)
set_key(Sys.getenv("GOOGLE_KEY"))

#use small sample od 10 linestrings from slopes package
Roads = st_transform(lisbon_road_segments, 4326) %>% slice(20:30)
mapview::mapview(Roads)

Roads_coords = data.frame(st_coordinates(Roads$geom))
names(Roads_coords) = c("lon", "lat", "l")

google = google_elevation(Roads_coords)
Roads_coords$elev = google$results$elevation
# Roads_coords$resol = google$results$resolution

Roads_xyz = Roads
Roads_xyz$geom= NULL
Roads_xyz$bbox= NULL

l=list()
for (i in 1:nrow(Roads)){ 
  line = Roads_coords[Roads_coords$l == i,]

  l[i] = st_as_sf(line, coords = c("lon","lat","elev")) %>% st_combine() %>%  st_cast("LINESTRING")

  Roads_xyz$geom[i] = l[i] 
}
#> Warning: Unknown or uninitialised column: `geom`.

Roads_xyz = st_as_sf(Roads_xyz, dim = "XYZ", crs=st_crs(Roads))

slope_xyz(Roads_xyz)
#> Loading required namespace: geodist
#>          1          2          3          4          5          6          7 
#> 0.14818945 0.08901145 0.19950000 0.07891966 0.10043988 0.07550628 0.11870625 
#>          8          9         10         11 
#> 0.09655327 0.09932231 0.05408169 0.13421425
Roads$slope = slope_xyz(Roads_xyz)

mapview::mapview(Roads["slope"])

[temospena]

Due the api request limits, it needs another loop to make one request at a time per linestring.

So, for the full lisbon_road_segments dataset, we can do:

library(tidyverse)
library(sf)
library(slopes)
library(googleway)
set_key(Sys.getenv("GOOGLE_KEY"))

Roads = st_transform(lisbon_road_segments, 4326)
# mapview::mapview(Roads)

#a df with only corddinates in WGS84 (google does not like others)
Roads_coords = data.frame(st_coordinates(Roads$geom))
names(Roads_coords) = c("lon", "lat", "l")

#a df without geometries
Roads_xyz = Roads
Roads_xyz$geom= NULL 
Roads_xyz$bbox= NULL
Roads_xyz$geom= NA

#get elevarion per point, and transform each line in a xyz linestring again
for (i in 1:nrow(Roads)){ 
  google = google_elevation(Roads_coords[Roads_coords$l == i,])
  Roads_coords$elev[Roads_coords$l == i] = google$results$elevation
  Roads_coords$resol[Roads_coords$l == i] = google$results$resolution

  line = Roads_coords[Roads_coords$l == i,]
  Roads_xyz$geom[i] = st_as_sf(line, coords = c("lon","lat","elev")) %>% st_combine() %>%  st_cast("LINESTRING")
}

Roads_xyz = st_as_sf(Roads_xyz, dim = "XYZ", crs=st_crs(Roads))

Roads$slope = slope_xyz(Roads_xyz) #to store the values in the original dataset

summary(Roads$slope)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> 0.0004856 0.0064922 0.0279617 0.0455289 0.0780702 0.1995000 

summary(Roads_coords$resol) #check average resolution
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   1.193   1.193   1.193   1.193   1.193   1.193 

mapview::mapview(Roads["slope"])

[Robinlovelace]

Fantastic, I think that could be the basis of a new function to get slopes from Google's API. Hypothesis: will be faster and more accurate than the current approach that downloads raster DEMs covering the entire area.