问题描述
我正在尝试使用开放数据(避免像Google这样的许可限制)来计算远足路线的海拔数据.
I'm trying to calculate elevation data for hiking routes, using open data (avoiding licensing constraints like Google).
我能够读取我所在国家/地区的公开DEM(分辨率为10米)使用readGDAL(来自RGDAL软件包)和proj4string(mygrid)给我:
I was able to read a public DEM of my country (with a 10-metres resolution)using readGDAL (from package RGDAL), and proj4string(mygrid) gives me:
"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
.asc文件的开头是:
The beginning of the .asc file is:
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 ....
.....
我需要做的就是从此网格中获取路线特定节点的高程数据,以便能够计算高程增益,负斜率,最小/最大高度...
All that I need is to pick up from this grid the elevation data for the specific nodes of a route, so to be ableto calculate elevation gain, negative slope, min/max altitude...
我使用不错的OSMAR包从OpenStreetMap带来了路线数据,所以我路线的数据表是这样的:
I bring the routes data from OpenStreetMap, using the nice package OSMAR, so the data table of my route is something like this:
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坐标参考系统中变换纬度/经度坐标,然后如何带来相应的网格值(对最近点进行平均?)
But I've no idea how to transform the latitude/longitude coordinates in the DEM coordinate reference system, and then how to bring the corresponding grid values (doing a sort of average of the nearest points?)
我发现所有使用谷歌搜索的文档都将渲染网格图,而不是从中提取值.
All the documentation I've found googling it is about to render grid maps, not to extract values from them.
任何帮助将不胜感激!
干杯,MB
P.S.第二个问题应该是:具有多个网格拼贴,如果一条路线跨越两个或多个拼贴,该怎么办?将它们合并,同时引用两者..."
P.S. The second question should be:"Having several grid-tiles, what could I do if a route is across two or more tiles? Merge them, reference both..."
推荐答案
请勿转换DEM.改变你的观点.您的DEM投影在常规网格上.如果您将其重新投影到另一个项目中,则可能不会.而是从OSM数据中得出空间点:
Don't transform the DEM. Transform your points. Your DEM is projected over a regular grid. If you re-project it into another it may not be. Instead, make spatial points from your OSM data:
require(sp); coordinates(my.OSM.points) <- long + lat
然后将它们转换为DEM的坐标参考系统(请注意,您可能需要首先设置点的CRS.因此,您可以这样做:
Then transform them to the coordinate reference system of the DEM (note you may need to set the CRS of the points first. So you can do this:
# 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值.
Then look at using sp::over
or raster::extract
to get the values of the DEM at point locations.
这篇关于如何从DEM(数字高程模型)中提取特定值?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!