从网站提取澳大利亚经纬度点的高程，使用R语言。

Question

大家好，

我正在尝试获取大约700个点的高程数据。我想使用提供给同样问题的代码（R中纬度/经度转高度的转换），但不幸的是，在使用geonames包时会出现错误，并且最佳答案提供的网站没有澳大利亚高程数据可用（下面提供了错误信息以供参考）。

我找到了另一个网站，该网站为澳大利亚提供了非常准确的高程数据，但我不知道如何从网页中提取信息。我认为它正在使用google高程API，但我不知道如何访问它。

当我将“纬度，经度”坐标放入“搜索位置”框中时，它会在地图下方提供高程数据。但是，我似乎找不到这个数据在源页面中的位置。该网站是http://www.daftlogic.com/sandbox-google-maps-find-altitude.htm。

一些示例经纬度值：

-36.0736, 146.9442

-36.0491, 146.4622

我想知道是否有人可以帮助我从R查询此网站，并提取高程数据？或者这似乎太麻烦了吗？我意识到该网站有一个批处理功能（最多100个位置），但能够从R中完成这项工作将是很不错的。

谢谢大家，如果这显然明显请原谅我。

谢谢， Adam

错误信息：

当使用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'

使用查询代码时：

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

- Adam

---

如果不是因为你的问题，我就不会知道有Google Elevation API这个东西！太棒了！ - jbaums

我也是，那是一个非常有用的信息来源。 - JeremyS

Answer 1

谷歌提供了一个海拔 API (Elevation API)，可以返回 JSON 或 XML 响应。以下是使用 JSON 响应的示例，通过 RJSONIO 包中的 fromJSON 进行解析。

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  

< p > googEl 函数需要一个包含经度在第一列，纬度在第二列的坐标矩阵或数据框。

- jbaums

---

非常好的答案！非常感谢，它完美地解决了我的问题。祝你好运，亚当。 - Adam

我修改了这个程序，如果你输入一个命名的数据框，它也会输出位置的名称。修改后的代码如下（在评论中无法正确格式化，因此我添加了一个\n来表示关键的新行）： 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) - JeremyS

我只是想指出，这种方法有一个查询限制，大约为130K（截至此评论）。Google似乎更喜欢您使用他们的API密钥（请参见https://console.developers.google.com/apis/credentials），最好使用他们的JS API（请参见此处：https://developers.google.com/maps/documentation/javascript/elevation）。 - EconomiCurtis

Answer 2

在raster包中有?getData函数用于获取SRTM高程数据。

例如：

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

使用插值法而不是最近邻法来填充单元格：

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

下载步骤只需在工作目录中查找之前保存的文件，因此只需执行一次即可。

数据对象是一个 `RasterLayer`，可以使用 `plot` 等工具进行可视化。

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)

- mdsumner

Answer 3

我编写了googleway包，以使@jbaums的回答更易于实现。对于高程数据，您可以使用google_elevation()函数。

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

而且，为了好玩，您还可以在地图上标注这些位置。

key <- "your_api_key"

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

- SymbolixAU

---

这真的很有帮助。如果我超过请求限制，会得到什么错误？有什么方法可以处理这些错误吗？ - GabrielMontenegro

@JavierM88，如果您超过限制，Google将不会在隔天之前提供任何更多数据。Google有付费选项，可以增加限制。 - SymbolixAU