Is there a way to quickly get elevation/altitude given longitude and latitude?

General Usage
1

I have a data frame with longitude and latitude, and I would like to obtain the elevation/altitude for each point. At the moment, I use the R package geonames, which is a functional API for geonames.org.

For example:

using RCall

lon = 10
lat = -5 

@rput lon lat
R"""
library(geonames)
options(geonamesUsername = "username")
elev = unlist(GNsrtm3(lat, lon)[1])
"""
@rget elev

However, this method is incredibly time-consuming. I have data-frames with thousands of locations, and each GNsrtm3() takes 50+ seconds.

Do you know of a way that I can obtain the elevation/altitude quicker - preferably within Julia?

2

The Geonames.org web interface is really simple for elevation data. You could roll your own query in a couple of lines of code using HTTP.jl.

If you choose to have the data returned in JSON format, use JSON3.jl to parse it out.

2 Likes
3

GMT is very fast in doing that but one need to have the grid with the altitudes stored locally.
From the name GNsrtm3 I assume the data is from the SRTM 3 arc seconds (~90 meters resolution) grid.

The GMT site stores several grids at different resolutions that can be downloaded automatically. But because high res data is very big we should download only an area of interest. The command

using GMT
gmt("docs data")

opens the site with the documentation about the available data sets and instructions on how to access with plain GMT (shell commands) instructions, but it’s actually very easy to do the same with the GMT.jl interface. For example, for your example case one will do

julia> G = grdcut("@earth_relief_03s", limits=(9,11,-6,-4));
GMT [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
GMT [WARNING]: gmt_get_dataset_tiles: No earth_relief_03s_g tiles available for your region.
grdblend [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

grdblend [NOTICE]: SRTM15 Earth Relief original at 15x15 arc seconds [Tozer et al., 2019].
grdblend [NOTICE]:   -> Download 10x10 degree grid tile (earth_relief_15s_p): S10E000

The download will be done only once and the data is stored in your computer for later use. Now to get the altitudes we use the grdtrack module that accepts a Mx2 matrix with lon,lat (I think DataFrames will work as well)

But mind you if you are on Windows. There is a strange bug that only shows up on Windows, when the limits selected includes both land and ocean. In that case the best is to use a 15 arc seconds resolution grid (there are no 3 arc seconds global grids for the oceans).

julia> grdtrack([10 -5], G)
1×3 GMTdataset{Float64, 2}
 Row │     Lon      Lat         Z
     │ Float64  Float64   Float64
─────┼────────────────────────────
   1 │    10.0     -5.0  -2955.12

On the second run, that is after compiling,

julia> @time grdtrack([10 -5], G);
  0.067492 seconds (1.11 k allocations: 70.727 KiB)

and if we had passed a thousands points the time would have been only slightly longer.

4 Likes
4

This works perfectly! Thank you so much for your answer.

1 Like