hikereport

All the statistics and profiles of my hikes are obtained with a bit of computer work. What you need is two things: the data and the GIS software to work on that data.

The data

In my case I have obtained all the data from institutional sources (Regione Lombardia). They provide several maps and databases; they can send you cdroms for free, which I think is very cool… you don’t even pay for the shipping! In particular for my work I use the CTR (“Carta Tecnica Regionale”, 4 cdroms) at the highest resolution available (1:10000) and the Digital Elevation Model, that has a different resolution, with squares of 20×20 meters.

The GIS software

I have chosen GRASS because it is very mature, extensible, open source, and runs on Linux. In one word it is just great. Even if I am not at all a GIS expert I enjoy using it.

Preliminary steps

So, back to how to obtain the profiles. The first step is, of course, to create locations and mapsets and import the selected maps. It’s also good to patch some of the DEMs in order to have bigger maps. Once you have imported the needed maps, you can enjoy the processing

Processing

Actually you can draw a simple profile with GRASS’ command

d.profile

but that works only on a straight line, from point A to point B… but you don’t usually walk like that, do you?

No problem, GRASS has another command,

r.profile

that can accept a series of points (vertices of a polyline) as input, and outputs the profile (coordinates and height) along that polyline with a given resolution. Once we have that, the trick is nearly done… the part that’s missing is to use

r.profile

’s output to compute cumulative uphill, downhill, slopes and so on. That’s what I do with my

hikereport

script.

hikereport.py

hikereport

is a simple script that I have written. In the beginning it was a bash script, because I didn’t know python (you can find the bash version here). Unfortunately it was not portable, because it used arrays and arrays are not supported by

ash

. Moreover the bash version has dependencies (

awk, bc

) and runs very slow (because of all the piping with external tools like

bc

, for example).

So I wanted to make it portable and faster, and I tried a little bit with Tcl… but I didn’t like it at all. I don’t know why but it just looks awkward and ugly to me…

After that I looked at python and I think that this script, although simple, shows how you can immediately be productive with python. python is really worth learning. I just read the online tutorial and I started coding. Lots of fun, in other words!

Here is the code. You can also download it

  1 #!/usr/bin/python
  2 
  3 """
  4  MODULE:       hikereport
  5 
  6  AUTHOR(S):    Stefano Negri <ste.negri AT gmail.com>
  7 
  8  PURPOSE:      let the user draw a path, then report path's length,
  9                cumulative uphill climb, downhill climb, highest and lowest 
 10                 altitude, maximum and average uphill and downhill slope.
 11 
 12  COPYRIGHT:    (c) 2007 Stefano Negri
 13 
 14                This program is free software under the GNU General Public
 15                License (>=v2). Read the file COPYING that comes with GRASS
 16                for details.
 17 
 18  updates:
 19 
 20 24 December 2007: first version (bash)
 21 07 April 2008: first Python version
 22 08 April 2008: added "slope" option
 23 """
 24 
 25 import sys
 26 import os
 27 import subprocess
 28 import math
 29 
 30 class ParsingError(Exception):
 31     """
 32     Exception for errors while parsing r.profile's output
 33     """
 34     def __init__(self, value):
 35         self.value = value
 36     def __str__(self):
 37         return repr(self.value)
 38 
 39 def gmessage(message):
 40     """
 41     simple remapping of g.message
 42     """
 43     try:
 44         subprocess.call(["g.message", message])
 45     except OSError:
 46         sys.stderr.write(message + "\n" )
 47 
 48 def mapExists(mapname):
 49     """
 50     check if the map exists in current mapset (if so it returns True)
 51     """
 52     try:
 53         if subprocess.call(["g.findfile", "element=cell", "file=%s" % mapname]) != 0:
 54             gmessage('Raster map %s not found in mapset search path' % mapname)
 55             return False
 56         else:
 57             return True
 58     except OSError:
 59         return False
 60 
 61 def getTempFile():
 62     """
 63     ask grass a temporary file name and return it
 64     """
 65     return os.popen('g.tempfile pid=%d' % os.getpid()).readlines()[0]
 66 
 67 def getResolution(gmap):
 68     """
 69     get resolution from r.info for the given map
 70     """
 71 
 72     res = os.getenv("GIS_OPT_resolution" )
 73     if res == '':
 74         line = os.popen('r.info -s map=%s' % gmap).readlines()[0]
 75         res = line.split("=" )[1][:-1] #[:-1] removes trailing newline
 76     return res
 77 
 78 def getLen(line):
 79     """
 80     parses a line of r.profile's output and extracts the transect length
 81 
 82     >>> print getLen("a b 100 856" )
 83     100.0
 84 
 85     >>> print getLen("100 856" )
 86     100.0
 87     """
 88     fields = line.split(" " )
 89     if len(fields) == 2:
 90         return float(fields[0])
 91     elif len(fields) == 4:
 92         return float(fields[2])
 93     else:
 94         raise ParsingError("%s has an invalid number of fields" % line)
 95 
 96 def getZ(line):
 97     """
 98     parses a line of r.profile's output and extracts the height coordinate
 99 
100     >>> print getZ("a b 100 856" )
101     856.0
102 
103     >>> print getZ("100 856" )
104     856.0
105     """
106     fields = line.split(" " )
107     if len(fields) == 2:
108         return float(fields[1])
109     elif len(fields) == 4:
110         return float(fields[3])
111     else:
112         raise ParsingError("%s has an invalid number of fields" % line)
113 
114 def average(array):
115     """
116     given an array of numeric values, it returns its average
117     
118     >>> print average([20, 30, 70])
119     40
120     """
121     return sum(array)/len(array)
122 
123 def cumulativeUpHill(array):
124     """
125     given an altitude profile, expressed as an array of numeric values,
126     it returns the cumulative uphill
127 
128     >>> print cumulativeUpHill([10, 20, 30.5])
129     20.5
130 
131     >>> print cumulativeUpHill([10, 9, 8.1, 8.2, 7, 0, 4])
132     4.1
133 
134     >>> print cumulativeUpHill([2])
135     0
136     """
137     if len(array) < 2:
138         return 0
139 
140     i, j = 0, 1
141     total = 0
142     while j < len(array):
143         if array[j] > array[i]:
144             total += array[j] - array[i]
145         i, j = i+1, j+1
146 
147     return total
148 
149 def cumulativeDownHill(array):
150     """
151     given an altitude profile, expressed as an array of numeric values,
152     it returns the cumulative downhill
153 
154     >>> print  cumulativeDownHill([10, 20, 30.5])
155     0
156 
157     >>> print cumulativeDownHill([10, 9, 8.1, 8.2, 7, 0, 4])
158     10.1
159 
160     >>> print cumulativeDownHill([2])
161     0
162     """
163     if len(array) < 2:
164         return 0
165 
166     i, j = 0, 1
167     total = 0
168     while j < len(array):
169         if array[j] < array[i]:
170             total += array[i] - array[j]
171         i, j = i+1, j+1
172 
173     return total
174 
175 def extractArrays(outfile):
176     """
177     Extract an array of z coordinates and an array of lengths
178     from r.profile's output (contained in the file 'outfile')
179     """
180     lengthArr = []
181     heightArr = []
182     profile = open(outfile, 'r')
183     for curline in profile:
184         lengthArr.append(getLen(curline))
185         heightArr.append(getZ(curline))
186 
187     assert len(lengthArr) == len(heightArr), \
188             "An error has occurred while processing r.profile's output"
189 
190     #uncomment for debug
191 #    i = 0
192 #    for x, y in zip(lengthArr, heightArr):
193 #        print("%d: %f    %f" % (i, x, y))
194 #        i = i+1
195 
196     return heightArr, lengthArr
197 
198 def computeSlope(h, l, s):
199     """
200     return slope on single segment, 
201     expressed in % or degs according to parameter s
202 
203     >>> print computeSlope(10, 10, 'percent')
204     100
205 
206     >>> print computeSlope(10, 10, 'degrees')
207     45.0
208     """
209     if s == "percent":
210         return h / l * 100
211     if s == "degrees":
212         return math.degrees(math.atan(h / l))
213 
214 def buildSlopeArrays(ZArr, lenArr, slope):
215     """
216     Build an array of segments with upslope and downslope,
217     given ZArr, an array of heights, and lenArr, an array of lengths.
218     Resulting array contain slopes expressed in % or degrees, 
219     according to the value of parameter "slope"
220 
221     len(ZArr) == len(lenArr) while this is not necessarily true for
222     upSlopeArr and downSlopeArr
223     """
224 
225     assert (slope == "percent" or slope == "degrees" ), \
226             "slope preference incorrectly set"
227     upSlopes = []
228     downSlopes = []
229     i, j = 0, 1
230     #uncomment for debug:
231 #    print("%4s, %4s: (%5s - %11s), (%5s - %11s); slope percentage" % \
232 #        ("i", "j", "len", "z coord", "len", "z coord" ))
233 #    print("-" * 80)
234     while j < len(ZArr):
235         assert lenArr[j] != lenArr[i], "Corrupted profile"
236         segLen = lenArr[j] - lenArr[i]
237         if ZArr[j] > ZArr[i]:
238             heightDiff = ZArr[j] - ZArr[i]
239             upSlopes.append(computeSlope(heightDiff, segLen, slope))
240             #uncomment for debug:
241 #            print("%4d, %4d: (%5d - %f), (%5d - %f); up=%f" % \
242 #                (i, j, lenArr[i], ZArr[i], lenArr[j], ZArr[j], upSlopes[-1]))
243         else:
244             heightDiff = ZArr[i] - ZArr[j]
245             downSlopes.append(computeSlope(heightDiff, segLen, slope))
246             #uncomment for debug:
247 #            print("%4d, %4d: (%5d - %f), (%5d - %f); down=%f" % \
248 #                (i, j, lenArr[i], ZArr[i], lenArr[j], ZArr[j], downSlopes[-1]))
249         i, j = i+1, j+1
250     return upSlopes, downSlopes
251 
252 def main():
253     DEM = os.getenv("GIS_OPT_DEM" )
254     if not mapExists(DEM):
255         gmessage('Raster map %s not found in mapset search path' % DEM)
256         sys.exit(1)
257 
258     resolution = os.getenv("GIS_OPT_resolution" )
259     if resolution == '':
260         resolution = getResolution(DEM)
261         gmessage('using DEM map resolution: %s' % resolution)
262 
263     slope = os.getenv("GIS_OPT_slope" )
264     if not (slope == "percent" or slope == "degrees" ):
265         slope = "degrees"
266         gmessage('slope preference not set (available: "percent", "degrees" ): \
267                 using degrees')
268     if slope == "degrees":
269         slopeSymbol = "deg"
270     elif slope == "percent":
271         slopeSymbol = "%"
272 
273     outfile = getTempFile()
274     if outfile == '':
275         gmessage('Unable to create temporary file')
276         sys.exit(1)
277 
278     try:
279         subprocess.call(["r.profile", "-gi", "res=%s" % resolution, \
280                 "input=%s" % DEM, "output=%s" % outfile])
281         lenArr = [] #an array where all the segment lengths will be stored
282         ZArr = []   #an array where all the z coordinates will be stored
283         ZArr, lenArr = extractArrays(outfile)
284 
285         upSlopeArr = []   #up slope array
286         downSlopeArr = [] #down slope array
287         upSlopeArr, downSlopeArr = buildSlopeArrays(ZArr, lenArr, slope)
288 
289         minHeight = min(ZArr)
290         maxHeight = max(ZArr)
291         pathLength = lenArr[-1]
292         cumulativeUphill = cumulativeUpHill(ZArr)
293         cumulativeDownhill = cumulativeDownHill(ZArr)
294         maxUpSlope = max(upSlopeArr)
295         maxDownSlope = max(downSlopeArr)
296         avgUpSlope = average(upSlopeArr)
297         avgDownSlope = average(downSlopeArr)
298 
299         ##################################################
300         # output
301         ##################################################
302         gmessage("Path length:               %d" % pathLength)
303         gmessage("Cumulative uphill climb:   %d" % cumulativeUphill)
304         gmessage("Cumulative downhill climb: %d" % cumulativeDownhill)
305         gmessage("Highest altitude:          %d" % maxHeight)
306         gmessage("Lowest altitude:           %d" % minHeight)
307         gmessage("Maximum uphill slope:      %.1f%s" % \
308                 (maxUpSlope, slopeSymbol))
309         gmessage("Maximum downhill slope:    %.1f%s" % \
310                 (maxDownSlope, slopeSymbol))
311         gmessage("Average uphill slope:      %.1f%s" % \
312                 (avgUpSlope, slopeSymbol))
313         gmessage("Average downhill slope:    %.1f%s" % \
314                 (avgDownSlope, slopeSymbol))
315         gmessage("Profile:                   %s" % \
316                 outfile)
317 
318     except OSError:
319         gmessage("A fatal error has occurred while processing r.profile" )
320         sys.exit(1)
321     except ParsingError, pe:
322         gmessage(pe.value)
323         sys.exit(1)
324 
325 
326 if __name__ == "__main__":
327 
328     if os.getenv("GISBASE" ) == None:
329         gmessage("You must be in GRASS GIS to run this program" )
330         #if invoked from outside GRASS... at least run the tests  
331         gmessage("...running tests for this module:" )
332         import doctest
333         failures, tests = doctest.testmod(report=1)
334         gmessage("    %d failures in %d tests" % (failures, tests))
335         sys.exit(1)
336 
337     if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
338         os.execvp("g.parser", [sys.argv[0]] + sys.argv)
339     else:
340         main()
341 
342 
343 #%Module
344 #%  description: computes lenght and cumulative up and down climb for a given path
345 #%End
346 #%option
347 #%  key: DEM
348 #%  type: string
349 #%  gisprompt: old,cell,raster
350 #%  description: DEM map
351 #%  required : yes
352 #%end
353 #%option
354 #%  key: resolution
355 #%  type: integer
356 #%  description: resolution (default: DEM's resolution)
357 #%  required : no
358 #%end
359 #%option
360 #%  key: slope
361 #%  type: string
362 #%  description: express slope in percent or degrees? (default: degrees)
363 #%  required : no
364 #%end

I think (hope) it’s quite self explanatory. Something worth noting is the python way to do unit testing: you just write calls and expected results, see for example line 157, 158. Writing test code for such small modules might look useless but I can assure it is not. In my opinion unit tests are always a great help, and in python it’s really easy to write simple ones.

Usage

Anyway, how do you use it? First of all you must have a display open on a map of your choice. The script doesn’t open a display for you, it just fails if you don’t have one open. This is actually because the script relies on

r.profile

and that’s

r.profile

’s behaviour.

You can call the script giving three options: DEM, resolution and slope. Only DEM is mandatory.

DEM

the digital elevation model that contains the elevation information. Of course this is not necessarily the same map you have on display (for example I always use the CTR with 1:10000 scale on display, because it is so detailed and it’s easy to draw my path on it)

resolution

the resolution to use along the path; defaults to DEM’s resolution

slope

this has to be “percent” or “degrees”, one of these two strings. You can use it to specify (guess…) if you want slopes expressed in percent or degrees. It defaults to degrees

Once you issue the command you can start drawing your path on the screen, by left-clicking. When you right click your input session ends and the program outputs the results.