How to access multiple elevation points at the same time with get_elev_point

425 Views Asked by At

I hava another geodata question. I am trying to access the elevation data for specific points with the get_elev_point ()

I have about 160 points that I would like to get the elevation data for, so I would like to not do it one by one but I do not know how to retrieve them all at once.

I have made a data.frame of all my points where the first column x is longitude and the second column y is latitude and these are the only data in the data.frame.

I am using the elevatr package

Reproducible example:

location_1 <- data.frame (x=-7.37,y=5.775)

location_1_elev <-get_elev_point(locations = location_1, units="meters", prj = ll_prj, src = "aws")

When I do this, everything is fine and I get an elevation point back, but when I try to access several points at the same time I run into errors.

I took the earthquake data from R and transformed it into a data.frame that has only the longitude and latitiude. Then tried to access the elevation points via get_elev_points and got the error message:

data(quakes)

head(quakes)

locations <- data.frame(x = c(quakes$long, 1000), y = c(quakes$lat, 1000))

quakes_elev <-get_elev_point(locations = locations, units="meters", prj = ll_prj, src = "aws")

Error: API did not return tif

Do you have any tips how to make this happen, to be able to access multiple elevation points?

Thank you! ps.: Sorry for my clumsy asking, I am only learning now.

2

There are 2 best solutions below

3
On BEST ANSWER

It's likely the error is caused because the points are in the sea so there is no land elevation tif file to get an elevation from. For example, here is an example that tries to find the elevation of a point in the sea and then a point on land.

    library(elevatr)
    ll_prj <- "EPSG:4326"
            
    sea <- data.frame(x=181.62, y=-20.42)
    # This errors
    sea_elev <- get_elev_point(locations = sea, units='meters', prj=ll_prj, src='aws')
    # Error: This url: https://s3.amazonaws.com/elevation-tiles-prod/geotiff/5/32/18.tif did not return a tif
    
    land <- data.frame(x=-71.3036, y=44.2700)
    # This works
    land_elev <- get_elev_point(locations = land, units='meters', prj=ll_prj, src='aws')
    land_elev$elevation
    # [1] 1478

Using the epqs option returns NA without error so I suppose you could replace that with 0 if you want to use sea level as the elevation.

0
On

I finally had time to look into this.

The root causes of the problem are 1) data near 180 degrees, and 2) longitude using 0 to 360.

The first issue was a problem as I was grabbing the next higher tile for a given longitude. The end result is an x/y/z tile that didn't exist and thus would return an error. I have a fix that only grabs the tile that corresponds with a given longitude and not the next higher. This is fixed and pushed to https://github/jhollist/elevatr.

The second issue is a problem as elevatr assumes longitude between -180 and 180. So if you convert the quakes dataset to this and use the newer version of elevatr on GitHub, @kamilla-choni-pléh code will work.

install.packages("remotes")
remotes::install_github("jhollist/elevatr")
library(elevatr)
ll_prj <- "EPSG:4326"
data(quakes)
locations <- data.frame(x = c(quakes$long), y = c(quakes$lat))
locations$x <- ifelse(locations$x >= 180, locations$x - 360, locations$x)
quakes_elev <-get_elev_point(locations = locations, units="meters", prj = ll_prj, src = "aws")
quakes_elev$elevation

If using a higher zoom level, though it may complain as it is currently creating a raster that spans -180 to 180 (i.e. the globe). I am still mulling over how I should deal with that in elevatr.