Skip to content

add google elevation source? #28

Open

Description

@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"])

Activity

  1. temospena commented on Apr 25, 2021

    @temospena
    CollaboratorAuthor

    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"])

    image

  2. Robinlovelace commented on Apr 25, 2021

    @Robinlovelace
    Member

    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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions