2008/1/3, stefano negri <[EMAIL PROTECTED]>: > Hello, > I have written for my use a simple script that lets the user draw a path on a > map (like r.profile, on which the script is built) and then outputs the > path's length, cumulative uphill and downhill climb, highest and lowest > altitudes, maximum and average slope (up and down). > Unfortunately I was not able to make it work with ash, because (in my > understanding) ash does not support arrays. I am also thinking about > rewriting everything in C or Tcl. > Te script is attached to this email, in case it can be of any interest to > anybody. > Regards and happy new year, > Stefano.
Hello, there is now a python version of basically the same script (and also removed a bug that was there in the bash version). Both python and corrected bash version are attached. The python version, of course, does not depend on external commands (awk, bc...) and is incomparably faster. Ideas and comments are welcome. If you think it could be worth going into the addons, I would be glad to add it there, but I don't have write access to the repository. As I have read some posts in the past about scripting GRASS with bash, Tcl and python, I feel I can give my 2 cents on that (for others that are considering learning one of the three), because before deciding for python I have also tried to learn some Tcl. I don't want to comment on the language itself, just on how easy it can be to learn. On this respect I think python is for sure the best to choose: I found it very natural, while I had exactly the opposite feeling with Tcl. That's just *my* experience... it could also be a matter of what kind of background you have: mine is mainly in Java and Delphi (but I have also used in the past C, a bit of Lisp, and bash). Best regards, Stefano.
#!/usr/bin/python
"""
MODULE: hikereport
AUTHOR(S): Stefano Negri <ste.negri AT gmail.com>
PURPOSE: let the user draw a path, then report path's length,
cumulative uphill climb, downhill climb, highest and lowest
altitude, maximum and average uphill and downhill slope.
COPYRIGHT: (c) 2007 Stefano Negri
This program is free software under the GNU General Public
License (>=v2). Read the file COPYING that comes with GRASS
for details.
updates:
24 December 2007: first version (bash)
07 April 2008: first Python version
08 April 2008: added "slope" option
"""
import sys
import os
import subprocess
import math
class ParsingError(Exception):
"""
Exception for errors while parsing r.profile's output
"""
def __init__(self, value):
self.value = value
def __str__(self):
return repr(self.value)
def gmessage(message):
"""
simple remapping of g.message
"""
try:
subprocess.call(["g.message", message])
except OSError:
sys.stderr.write(message + "\n")
def mapExists(mapname):
"""
check if the map exists in current mapset (if so it returns True)
"""
try:
if subprocess.call(["g.findfile", "element=cell", "file=%s" % mapname]) != 0:
gmessage('Raster map %s not found in mapset search path' % mapname)
return False
else:
return True
except OSError:
return False
def getTempFile():
"""
ask grass a temporary file name and return it
"""
return os.popen('g.tempfile pid=%d' % os.getpid()).readlines()[0]
def getResolution(gmap):
"""
get resolution from r.info for the given map
"""
res = os.getenv("GIS_OPT_resolution")
if res == '':
line = os.popen('r.info -s map=%s' % gmap).readlines()[0]
res = line.split("=")[1][:-1] #[:-1] removes trailing newline
return res
def getLen(line):
"""
parses a line of r.profile's output and extracts the transect length
>>> print getLen("a b 100 856")
100.0
>>> print getLen("100 856")
100.0
"""
fields = line.split(" ")
if len(fields) == 2:
return float(fields[0])
elif len(fields) == 4:
return float(fields[2])
else:
raise ParsingError("%s has an invalid number of fields" % line)
def getZ(line):
"""
parses a line of r.profile's output and extracts the height coordinate
>>> print getZ("a b 100 856")
856.0
>>> print getZ("100 856")
856.0
"""
fields = line.split(" ")
if len(fields) == 2:
return float(fields[1])
elif len(fields) == 4:
return float(fields[3])
else:
raise ParsingError("%s has an invalid number of fields" % line)
def average(array):
"""
given an array of numeric values, it returns its average
>>> print average([20, 30, 70])
40
"""
return sum(array)/len(array)
def cumulativeUpHill(array):
"""
given an altitude profile, expressed as an array of numeric values,
it returns the cumulative uphill
>>> print cumulativeUpHill([10, 20, 30.5])
20.5
>>> print cumulativeUpHill([10, 9, 8.1, 8.2, 7, 0, 4])
4.1
>>> print cumulativeUpHill([2])
0
"""
if len(array) < 2:
return 0
i, j = 0, 1
total = 0
while j < len(array):
if array[j] > array[i]:
total += array[j] - array[i]
i, j = i+1, j+1
return total
def cumulativeDownHill(array):
"""
given an altitude profile, expressed as an array of numeric values,
it returns the cumulative downhill
>>> print cumulativeDownHill([10, 20, 30.5])
0
>>> print cumulativeDownHill([10, 9, 8.1, 8.2, 7, 0, 4])
10.1
>>> print cumulativeDownHill([2])
0
"""
if len(array) < 2:
return 0
i, j = 0, 1
total = 0
while j < len(array):
if array[j] < array[i]:
total += array[i] - array[j]
i, j = i+1, j+1
return total
def extractArrays(outfile):
"""
Extract an array of z coordinates and an array of lengths
from r.profile's output (contained in the file 'outfile')
"""
lengthArr = []
heightArr = []
profile = open(outfile, 'r')
for curline in profile:
lengthArr.append(getLen(curline))
heightArr.append(getZ(curline))
assert len(lengthArr) == len(heightArr), \
"An error has occurred while processing r.profile's output"
#uncomment for debug
# i = 0
# for x, y in zip(lengthArr, heightArr):
# print("%d: %f %f" % (i, x, y))
# i = i+1
return heightArr, lengthArr
def computeSlope(h, l, s):
"""
return slope on single segment,
expressed in % or degs according to parameter s
>>> print computeSlope(10, 10, 'percent')
100
>>> print computeSlope(10, 10, 'degrees')
45.0
"""
if s == "percent":
return h / l * 100
if s == "degrees":
return math.degrees(math.atan(h / l))
def buildSlopeArrays(ZArr, lenArr, slope):
"""
Build an array of segments with upslope and downslope,
given ZArr, an array of heights, and lenArr, an array of lengths.
Resulting array contain slopes expressed in % or degrees,
according to the value of parameter "slope"
len(ZArr) == len(lenArr) while this is not necessarily true for
upSlopeArr and downSlopeArr
"""
assert (slope == "percent" or slope == "degrees"), \
"slope preference incorrectly set"
upSlopes = []
downSlopes = []
i, j = 0, 1
#uncomment for debug:
# print("%4s, %4s: (%5s - %11s), (%5s - %11s); slope percentage" % \
# ("i", "j", "len", "z coord", "len", "z coord"))
# print("-" * 80)
while j < len(ZArr):
assert lenArr[j] != lenArr[i], "Corrupted profile"
segLen = lenArr[j] - lenArr[i]
if ZArr[j] > ZArr[i]:
heightDiff = ZArr[j] - ZArr[i]
upSlopes.append(computeSlope(heightDiff, segLen, slope))
#uncomment for debug:
# print("%4d, %4d: (%5d - %f), (%5d - %f); up=%f" % \
# (i, j, lenArr[i], ZArr[i], lenArr[j], ZArr[j], upSlopes[-1]))
else:
heightDiff = ZArr[i] - ZArr[j]
downSlopes.append(computeSlope(heightDiff, segLen, slope))
#uncomment for debug:
# print("%4d, %4d: (%5d - %f), (%5d - %f); down=%f" % \
# (i, j, lenArr[i], ZArr[i], lenArr[j], ZArr[j], downSlopes[-1]))
i, j = i+1, j+1
return upSlopes, downSlopes
def main():
DEM = os.getenv("GIS_OPT_DEM")
if not mapExists(DEM):
gmessage('Raster map %s not found in mapset search path' % DEM)
sys.exit(1)
resolution = os.getenv("GIS_OPT_resolution")
if resolution == '':
resolution = getResolution(DEM)
gmessage('using DEM map resolution: %s' % resolution)
slope = os.getenv("GIS_OPT_slope")
if not (slope == "percent" or slope == "degrees"):
slope = "degrees"
gmessage('slope preference not set (available: "percent", "degrees"): \
using degrees')
if slope == "degrees":
slopeSymbol = "deg"
elif slope == "percent":
slopeSymbol = "%"
outfile = getTempFile()
if outfile == '':
gmessage('Unable to create temporary file')
sys.exit(1)
try:
subprocess.call(["r.profile", "-gi", "res=%s" % resolution, \
"input=%s" % DEM, "output=%s" % outfile])
lenArr = [] #an array where all the segment lengths will be stored
ZArr = [] #an array where all the z coordinates will be stored
ZArr, lenArr = extractArrays(outfile)
upSlopeArr = [] #up slope array
downSlopeArr = [] #down slope array
upSlopeArr, downSlopeArr = buildSlopeArrays(ZArr, lenArr, slope)
minHeight = min(ZArr)
maxHeight = max(ZArr)
pathLength = lenArr[-1]
cumulativeUphill = cumulativeUpHill(ZArr)
cumulativeDownhill = cumulativeDownHill(ZArr)
maxUpSlope = max(upSlopeArr)
maxDownSlope = max(downSlopeArr)
avgUpSlope = average(upSlopeArr)
avgDownSlope = average(downSlopeArr)
##################################################
# output
##################################################
gmessage("Path length: %d" % pathLength)
gmessage("Cumulative uphill climb: %d" % cumulativeUphill)
gmessage("Cumulative downhill climb: %d" % cumulativeDownhill)
gmessage("Highest altitude: %d" % maxHeight)
gmessage("Lowest altitude: %d" % minHeight)
gmessage("Maximum uphill slope: %.1f%s" % \
(maxUpSlope, slopeSymbol))
gmessage("Maximum downhill slope: %.1f%s" % \
(maxDownSlope, slopeSymbol))
gmessage("Average uphill slope: %.1f%s" % \
(avgUpSlope, slopeSymbol))
gmessage("Average downhill slope: %.1f%s" % \
(avgDownSlope, slopeSymbol))
gmessage("Profile: %s" % \
outfile)
except OSError:
gmessage("A fatal error has occurred while processing r.profile")
sys.exit(1)
except ParsingError, pe:
gmessage(pe.value)
sys.exit(1)
if __name__ == "__main__":
if os.getenv("GISBASE") == None:
gmessage("You must be in GRASS GIS to run this program")
#if invoked from outside GRASS... at least run the tests :-)
gmessage("...running tests for this module:")
import doctest
failures, tests = doctest.testmod(report=1)
gmessage(" %d failures in %d tests" % (failures, tests))
sys.exit(1)
if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
os.execvp("g.parser", [sys.argv[0]] + sys.argv)
else:
main()
#%Module
#% description: computes lenght and cumulative up and down climb for a given path
#%End
#%option
#% key: DEM
#% type: string
#% gisprompt: old,cell,raster
#% description: DEM map
#% required : yes
#%end
#%option
#% key: resolution
#% type: integer
#% description: resolution (default: DEM's resolution)
#% required : no
#%end
#%option
#% key: slope
#% type: string
#% description: express slope in percent or degrees? (default: degrees)
#% required : no
#%end
r.hikereport
Description: Binary data
_______________________________________________ grass-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-dev
