如何从 DEM（数字高程模型）中提取特定值？答案

【问题标题】：How to extract specific values from a DEM (digital elevation model)?如何从 DEM（数字高程模型）中提取特定值？
【发布时间】：2015-01-16 19:28:11
【问题描述】：

我正在尝试使用开放数据计算远足路线的海拔数据（避免像 Google 这样的许可限制）。

我能够阅读我所在国家/地区的公共 DEM（分辨率为 10 米）使用 readGDAL（来自 RGDAL 包），proj4string(mygrid) 给了我：

"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

.asc 文件的开头是：

ncols         9000 
nrows         8884 
xllcorner     323256,181155 
yllcorner     4879269,74709 
cellsize      10 
NODATA_value  -9999 
978 998 1005 1008 1012 1016 1020 1025 ..... 
..... 
..... 400 Megabytes of elevation values .... 
.....

我只需要从这个网格中获取路线特定节点的高程数据，以便能够计算海拔增益、负斜率、最小/最大高度...

我使用漂亮的包 OSMAR 从 OpenStreetMap 带来路线数据，所以我的路线数据表是这样的：

    RouteId NodeId  lat         lon
1   -13828  -8754   45.36743    7.753664
2   -13828  -8756   45.36762    7.753878
3   -13828  -8758   45.36782    7.754344
4   -13828  -8760   45.36794    7.754541
....

但我不知道如何在 DEM 坐标参考系中转换纬度/经度坐标，然后如何带上相应的网格值（对最近点进行某种平均？）

我在谷歌上搜索到的所有文档都是关于渲染网格地图，而不是从中提取值。

任何帮助将不胜感激！

干杯，MB

附：第二个问题应该是： “有几个网格图块，如果一条路线跨越两个或多个图块，我该怎么办？合并它们，引用两者......”

【问题讨论】：

一种解决方案可能是 GraphHopper，其中已经获取了海拔数据（CGIAR 或 SRTM），您可以轻松计算远足路线。
感谢 Karussell，但 SRTM 数据的分辨率在 30 到 90 米之间：对于山区路线来说太多了...
好吧，您当然可以进行预处理和后处理，但您将获得一个现成且可扩展的解决方案。如果您有更好的数据：只需一个简单的 sn-p 即可将其输入 GraphHopper 并用于实际路线

标签： r openstreetmap rgdal osmar

【解决方案1】：

不要转换 DEM。转换你的观点。您的 DEM 投影在常规网格上。如果你将它重新投影到另一个它可能不是。相反，从您的 OSM 数据中创建空间点：

require(sp); coordinates(my.OSM.points) <- long + lat

然后将它们转换到 DEM 的坐标参考系（注意您可能需要先设置点的 CRS。所以您可以这样做：

#  Set pojection information for points
proj4string(my.OSM.points) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84")

#  Transform them
new.points <- spTransform( my.OSM.points , proj4string( mygrid ) )

然后看看使用sp::over 或raster::extract 来获取点位置的 DEM 值。

【讨论】：

【解决方案2】：

非常感谢西蒙，您的回答让我解决了我的问题！

只是一个标点符号：尝试使用您的代码，我得到了：

> coordinates(my.OSM.points) <- routesCoord[,lat,lon]
Error in coordinates(my.OSM.points) <- routesCoord[, lat, lon] : 
object 'my.OSM.points' not found

我认为您的陈述有效，因为“坐标”的帮助说： "设置空间坐标以创建一个空间对象"

但是我用过：

crsString <- "+proj=longlat +datum=WGS84 +ellps=WGS84"
my.OSM.points <- SpatialPoints( routesCoord[,lat,lon], 
                                proj4string=CRS(crsString), bbox = NULL)
new.points <- spTransform( my.OSM.points , CRS(proj4string(mygrid)) )

对于像我这样的新手：在此之后，您将获得一个漂亮的数据框，其中包含所有点的高程数据，只需使用

my.elev <- over(new.points, mygrid)

附：对不起西蒙，这是我在stackoverflow上的第一个问题，所以当我试图投票给你的答案时，我得到了“投票需要 15 的声誉”！现在我 11 岁了，只要我获得其他 4 分，我就回来了 :-)

P.S.2 我想知道“over”是取网格最近点的值，还是计算加权平均值最近的点...

P.S.3 第二个问题应该是：“有几个网格图块，如果路线跨越两个或多个图块，我该怎么办？合并它们，引用两者..."

【讨论】：

PS2 sp::over 采用该点所在的网格单元 - 如果您将网格单元视为“点”，只要我们在一个网格单元中，它就会采用最近的点（并且 NA如果我们在网格之外，或者网格单元格具有 NA 单元格值）。

【解决方案3】：

正如西蒙所说，如果你在同一个投影中，你可以使用：

 library(raster)
 my.pts.elevation <- extract(my.grid,my.point)

【讨论】：

谢谢@delaye，我也试过raster，它给了我和sp::over一样的结果（即使我有一点问题，我在第二个问题link中描述了它跨度>