Extracting elevation from website for lat/lon points in Australia, using R

Question

G'day Everyone,

I am trying to get some elevation data for about 700 points I have. I thought I might use the code provided to the same question (Conversion for latitude/longitude to altitude in R), unfortunately I get errors when using the geonames package, and the website the best answer provides does not have Australian elevation data available (Errors provided below FYI).

I found another website which provides very accurate elevation data for Australia, but I have no idea how I might extract the information from the webpage . I think it is using the google elevation API, but again I have no idea how to access that.

When I put 'lat, lon' coordinates into the 'search for location' box it gives the elevation data below the map. However, I can't seem to find this in the source page. The website is http://www.daftlogic.com/sandbox-google-maps-find-altitude.htm.

some example lon lat values that work:

-36.0736, 146.9442

-36.0491, 146.4622

I am wondering if anyone can help me query this site from R, and extract the elevation data? Or does it seem like to much of a hassle? I realise there is a batch function (up to 100 locations) on the website, but it would be cool to be able to do this from R.

Thanks everyone, sorry if this is extremely obvious.

Cheers, Adam

ERRORS

When using geonames:

elevation <- GNgtopo30(adult$lat, adult$lon)
Error in getJson("gtopo30JSON", list(lat = lat, lng = lng)) : 
  error code 10 from server: Please add a username to each call in order for geonames to    be able to identify the calling application and count the credits usage.
In addition: Warning message:
In readLines(u) :
  incomplete final line found on 'http://ws.geonames.org/gtopo30JSON?  lat=-36.0736&lng=146.9442'

When using query code:

library(RCurl)
library(XML)
url <- paste("http://earthtools.org/height", adult$lat, adult$lon, sep = '/')
page <- getURL(url)
ans <- xmlTreeParse(page, useInternalNodes = TRUE)
Space required after the Public Identifier
SystemLiteral " or ' expected
SYSTEM or PUBLIC, the URI is missing
Extra content at the end of the document
Error: 1: Space required after the Public Identifier
2: SystemLiteral " or ' expected
3: SYSTEM or PUBLIC, the URI is missing
4: Extra content at the end of the document

Answer 1

There's an Elevation API provided by Google, which returns either a JSON or XML response. Here's an example using a JSON response, parsed by fromJSON in the RJSONIO package.

googEl <- function(locs)  {
  require(RJSONIO)
  locstring <- paste(do.call(paste, list(locs[, 2], locs[, 1], sep=',')),
                     collapse='|')
  u <- sprintf('http://maps.googleapis.com/maps/api/elevation/json?locations=%s&sensor=false',
               locstring)
  res <- fromJSON(u)
  out <- t(sapply(res[[1]], function(x) {
    c(x[['location']]['lat'], x[['location']]['lng'], 
      x['elevation'], x['resolution']) 
  }))    
  rownames(out) <- rownames(locs)
  return(out)
}

m <- matrix(c(146.9442, 146.4622, -36.0736, -36.0491), nc=2)

googEl(m)

      lat     lng      elevation resolution
[1,] -36.0736 146.9442 163       152.7032  
[2,] -36.0491 146.4622 171.7301  152.7032  

The googEl function expects a matrix or data.frame of coordinates, with longitude in the first column and latitude in the second.

Answer 2

There is ?getData for SRTM elevation in the raster package.

For example:

library(raster)
m <- data.frame(lon = c(146.9442, 146.4622), lat = c(-36.0736, -36.0491))

x <- getData('alt', country = "AUS")

cbind(m, alt = extract(x, m))
       lon      lat alt
1 146.9442 -36.0736 164
2 146.4622 -36.0491 172

Use interpolation into the cell rather than nearest neighbour:

cbind(m, alt = extract(x, m, method = "bilinear"))
       lon      lat      alt
1 146.9442 -36.0736 164.9519
2 146.4622 -36.0491 172.1293

The download step simply finds the file that was saved previously in the working directory, so that only has to happen once.

The data object is a RasterLayer which you can visualize with plot etc.:

plot(x)
points(m)
x
class       : RasterLayer 
dimensions  : 5496, 5568, 30601728  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 112.8, 159.2, -54.9, -9.1  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : C:\temp\AUS_msk_alt.grd 
names       : AUS_msk_alt 
values      : -43, 2143  (min, max)