2

I'm trying to calculate elevation data for hiking routes, using open data (avoiding licensing constraints like Google).

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"

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 able to calculate elevation gain, negative slope, min/max altitude...

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

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.

Any help will be greatly appreciated!

Cheers, MB

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..."

www
  • 38,575
  • 12
  • 48
  • 84
mbranco
  • 101
  • 2
  • 8
  • One solution could be GraphHopper where elevation data already will be fetched (either CGIAR or SRTM) and you can calculate hiking routes easily. – Karussell Nov 21 '14 at 08:37
  • thamk you Karussell, but SRTM data has resolution between 30 and 90 meters: too much for mountain routes... – mbranco Nov 21 '14 at 10:11
  • Well, sure you can do pre- and postprocessing but you'll get a ready and scalable solution. And if you have better data: just a simple snippet to get this into GraphHopper and available for real routes – Karussell Nov 21 '14 at 20:38

3 Answers3

4

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

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

Then look at using sp::over or raster::extract to get the values of the DEM at point locations.

Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184
1

Thank you very much Simon, your answser addressed me to solve my problem!

Just a punctuation: trying to use your code, I got:

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

I thought your statement worked, because the help for "coordinates" says: "set spatial coordinates to create a Spatial object"

However I used:

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

For newbies like me: after this you get a nice data frame with the elevation data for all points simply with

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

P.S. Sorry Simon, this is my first question on stackoverflow, so when I tried to vote your answer, I got "vote up requires 15 of reputation" ! Now I'm 11, as soon as I gain other 4 points, I came back :-)

P.S.2 I'd like to know if "over" take the value of the nearest point of the grid, or compute a weighted average of the nearest points...

P.S.3 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..."

mbranco
  • 101
  • 2
  • 8
  • P.S.2 `sp::over` takes the grid cell the point falls in - if you conceive grid cells as "points", it takes the nearest point, as long as we are in a grid cell (and NA if we're outside the grid, or if the grid cell has a NA cell value). – Edzer Pebesma Feb 23 '15 at 20:17
1

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

 library(raster)
 my.pts.elevation <- extract(my.grid,my.point)
delaye
  • 1,357
  • 27
  • 45
  • Thank you @delaye, I tried raster too, it gave me the same results as sp::over (even if I had a little problem, I describe it in my second question [link](http://stackoverflow.com/questions/27060307/incorrect-values-returned-by-spover-function) – mbranco Nov 25 '14 at 09:34