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

3

7

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
Magdamagdaia answered 6/2, 2014 at 4:14 Comment(2)
I wouldn't have known about the Google Elevation API if it weren't for your question! Nice one!Daimyo
me either, that is a very useful source of infoDeck
D
10

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.

Daimyo answered 6/2, 2014 at 4:54 Comment(2)
I modified this to also spit out the names of the locations if you feed it named data.frame by changing (doesn't let me format properly in a comment so I added a \n for the crucial new line): ans <- t(sapply(res[[1]], function(x) { c(x[['location']]['lat'], x[['location']]['lng'], x['elevation'], x['resolution']) })) \n rownames(ans) <- rownames(locs) \n return(ans)Deck
I just wanted to note, this method has a query limit of about 130K (as of this commenting). It also appears google would prefer you to use their api key (see console.developers.google.com/apis/credentials), and preferably their js api (see here: developers.google.com/maps/documentation/javascript/elevation).Adverbial
M
7

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)

enter image description here

Mazard answered 6/2, 2014 at 7:20 Comment(0)
P
2

I've written the googleway package to make @jbaums' answer easier to implement. For eleveation data you can use the google_elevation() function

library(googleway)

## you need an API key
key <- "your_api_key"

## data:
df <- data.frame(lat = c(-36.0736, -36.0491),
                 lon = c(146.9442, 146.4622))

google_elevation(df_locations = df, 
                 location_type = "individual",
                 key = key)

# $results
# elevation location.lat location.lng resolution
# 1  163.0000     -36.0736     146.9442   152.7032
# 2  171.7301     -36.0491     146.4622   152.7032

And, just for kicks you can plot the locations on a map too

key <- "your_api_key"

google_map(data = df, key = key) %>%
    add_markers()

enter image description here

Penicillate answered 16/10, 2016 at 21:22 Comment(2)
This was really helpful. What error will I get if I exceed the requests limit? And any way to handle those?Carduaceous
@JavierM88 If you exceed the limit Google will simply not supply you any more data until the following day. Google have paid options where you can increase the limit.Penicillate

© 2022 - 2024 — McMap. All rights reserved.