2

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

我能够使用 readGDAL（来自 RGDAL 包）读取我国的公共 DEM（分辨率为 10 米），并且 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

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

4

3 回答 3

3

不要转换 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 值。

于 2014-11-19T10:35:53.237 回答
1

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

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

> 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)

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

PS2我想知道“over”是否取网格最近点的值，或者计算最近点的加权平均值......

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

于 2014-11-19T17:34:30.640 回答
1

As Simon said , if you are in the same projection you can use : 

 library(raster)
 my.pts.elevation <- extract(my.grid,my.point)
于 2014-11-20T08:34:39.230 回答