Python script for comparing "zonal statistics" implementations/algorithms

zonalStatBenchmark.py 10KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267
  1. #!/usr/bin/python
  2. # -*- coding: utf-8 -*-
  3. import subprocess
  4. import sys
  5. import os
  6. from collections import defaultdict
  7. import re
  8. class zonalBenchmark:
  9. '''
  10. The script will provide a straight forward way to compare
  11. different zonal statistic implementations.
  12. '''
  13. def __init__(self):
  14. '''
  15. Current state of implementation:
  16. 1. oft-stat
  17. 2. pktools
  18. 3. Python/GDAL (perrygeo - gist)
  19. 4. GRASS
  20. 5. R
  21. 6. SAGA
  22. ( ^ those are implemented)
  23. 7. QGIS
  24. 8. rasterstats
  25. 9. starspan
  26. 10. R - multicore
  27. 11. PostGIS
  28. 12. GeoTools
  29. ( ^ those not yet)
  30. '''
  31. def checkInstallStatus():
  32. '''
  33. Check if the selected provider(s) are available on the
  34. (linux) machine.
  35. '''
  36. def getArguments(self, args):
  37. '''
  38. Get the links to the data which has been provided by the CLI
  39. and hand it to the execution. Checks also if the data exists.
  40. '''
  41. # show help for no arguments
  42. if len(args[0]) < 2:
  43. zonalBenchmark().showUsage()
  44. sys.exit()
  45. try:
  46. if args[0][1] == '-h' or args[0][1] == '--help':
  47. zonalBenchmark().showUsage()
  48. sys.exit()
  49. except IndexError:
  50. zonalBenchmark().showUsage()
  51. sys.exit()
  52. # get arguments for tools, split them, convert to int and sort
  53. try:
  54. toolsArgs = args[0][1].split('-')
  55. tools = sorted(map(int, toolsArgs), key=int)
  56. except (TypeError, ValueError) as e:
  57. zonalBenchmark(args).showUsage()
  58. sys.exit()
  59. # dict of implemented tools
  60. availableTools = {1:'Open Foris Tools',
  61. 2: 'pktools',
  62. 3: 'perrygeo',
  63. 4: 'GRASS',
  64. 5: 'R',
  65. 6: 'SAGA'}
  66. selectedTools = {}
  67. # create dict for selected tools
  68. for item in tools:
  69. if item in availableTools:
  70. selectedTools.update({item: availableTools[item]})
  71. # get input file
  72. workingDir = os.path.dirname(os.path.realpath(__file__))
  73. try:
  74. inputArg = args[0][2]
  75. inputFile = os.path.join(workingDir, inputArg)
  76. if not os.path.isfile(inputFile):
  77. print 'ERROR: Input file does not exist\You can run \'{0} --help\' for usage of the program.\n'.format(str(args[0][0]).lstrip('./'))
  78. sys.exit()
  79. except (IndexError, ValueError) as e:
  80. zonalBenchmark().showUsage()
  81. # get mask file
  82. try:
  83. maskArg = args[0][3]
  84. maskFile = os.path.join(workingDir, maskArg)
  85. if not os.path.isfile(maskFile):
  86. print 'ERROR: Mask file does not exist\nYou can run \'{0} --help\' for usage of the program.\n'.format(str(args[0][0]).lstrip('./'))
  87. sys.exit()
  88. except (IndexError, ValueError) as e:
  89. zonalBenchmark().showUsage()
  90. # get number of repetitions
  91. try:
  92. repArg = args[0][4]
  93. except (IndexError, ValueError) as e:
  94. zonalBenchmark().showUsage()
  95. return selectedTools, inputFile, maskFile, repArg
  96. def setCommands(self, selectedTools, inputFile, maskFile):
  97. '''
  98. Constructs the actual command(s) for the selected providers.
  99. '''
  100. selectedCommands = {}
  101. if 1 in selectedTools:
  102. maskFileName = os.path.basename(maskFile)
  103. maskTifName = os.path.dirname(maskFile) + '/' + maskFileName[:-4] + ".tif"
  104. maskShpName = os.path.dirname(maskFile) + '/' + maskFileName[:-4] + ".shp"
  105. f = open('bash_oft.sh', 'w')
  106. f.write('''resX=$(pkinfo -i "{0}" -ns | cut -d ' ' -f 2);echo $resX; resY=$(pkinfo -i "{0}" -nl | cut -d ' ' -f 2);echo $resY;xmin=$(pkinfo -i "{0}" -ulx | cut -d '=' -f 2);echo $xmin;xmax=$(pkinfo -i "{0}" -lrx | cut -d '=' -f 2);echo $xmax;ymax=$(pkinfo -i "{0}" -uly | cut -d '=' -f 2);echo $ymax;ymin=$(pkinfo -i "{0}" -lry | cut -d '=' -f 2);echo $ymin;gdal_rasterize -ts $resX $resY -te $xmin $ymin $xmax $ymax -a_nodata 0 -ot Byte -burn 1 -l "{2}" "{1}" "{3}"'''.format(inputFile, maskShpName, maskFileName[:-4], maskTifName))
  107. selectedCommands.update({1:'bash_script'})
  108. if 2 in selectedTools:
  109. #~ maskFileShp = maskFile.replace('tif','shp')
  110. pkCMD = 'pkextract -f \'ESRI Shapefile\' -s {0} -i {1} -o /tmp/pk_out.shp -polygon --rule mean'.format(maskFile, inputFile)
  111. selectedCommands.update({2:pkCMD})
  112. if 6 in selectedTools:
  113. sagaCMD = 'saga_cmd statistics_grid 5 -ZONES {0} -STATLIST {1} -OUTTAB /tmp/saga_out.txt'.format(maskFile, inputFile)
  114. selectedCommands.update({6:sagaCMD})
  115. return selectedCommands
  116. def executeCommands(self, commands, repetitions):
  117. '''
  118. Executes the command for the selected provider.
  119. '''
  120. timing = {}
  121. secondsAvg = []
  122. cpuAvg = []
  123. #~ '''/usr/bin/time -f "seconds: %e \ncpu-load: %P" /bin/bash -c 'resX=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -ns | cut -d '\'' '\'' -f 2);echo $resX; resY=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -nl | cut -d '\'' '\'' -f 2);echo $resY;xmin=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -ulx | cut -d '=' -f 2);echo $xmin;xmax=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -lrx | cut -d '=' -f 2);echo $xmax;ymax=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -uly | cut -d '=' -f 2);echo $ymax;ymin=$(pkinfo -i test_data/wc2.0_10m_tavg_07.tif -lry | cut -d '=' -f 2);echo $ymin;gdal_rasterize -ts $resX $resY -te $xmin $ymin $xmax $ymax -a_nodata 0 -ot Byte -burn 1 -l mask test_data/mask.shp test_data/out.tif'
  124. #~ '''
  125. for key, value in commands.items():
  126. reps = 0
  127. for i in range(0,int(repetitions)):
  128. try:
  129. # check if the command should be run via .sh / bash script
  130. if commands[key] == 'bash_script':
  131. output = subprocess.check_output(
  132. '''/usr/bin/time -f "seconds: %e \ncpu-load: %P" /bin/bash -c ./bash_oft.sh | tail -2'''.format(value),
  133. shell=True,
  134. stderr=subprocess.STDOUT
  135. )
  136. else:
  137. print "nope?"
  138. except subprocess.CalledProcessError as grepexc:
  139. print "error code", grepexc.returncode, grepexc.output
  140. # grep seconds and cpu-load from the string
  141. seconds = float((re.split(":|\n", output[0:28].replace(" ", "")))[1])
  142. cpuload = int(((re.split(":|\n", output[0:28].replace(" ", "")))[3]).replace("%", ""))
  143. # set output to 0 to prevent any 'misappropriation' of older outputs
  144. output = "sec: 0 \n cpu:0"
  145. # increment reps and fill lists for timing & cpu load
  146. secondsAvg.append(seconds)
  147. cpuAvg.append(cpuload)
  148. reps += 1
  149. # update the dict with the average timing and cpuload, as well as a repetition counter
  150. timing.update({key:[round(sum(secondsAvg) / len(secondsAvg), 3), sum(cpuAvg) / len(cpuAvg), reps]})
  151. # clear the lists for the next values
  152. secondsAvg[:] = []
  153. cpuAvg[:] = []
  154. print timing
  155. sys.exit()
  156. return timing
  157. def showOutput(self, output, selectedTools):
  158. '''
  159. Construct a pretty table to show a comparison of the programs
  160. on the command line.
  161. '''
  162. # merge the names of the selected tools with the numerical representation
  163. # save as a new dict
  164. dd = defaultdict(list)
  165. for d in (selectedTools, output):
  166. for key, value in d.iteritems():
  167. dd[key].extend([value])
  168. print dd
  169. out = ''
  170. out += "All selected tools ran successfully.\n\n"
  171. # print with underscore
  172. out += "\033[4m{:<8} {:<17} {:<10} {:<15}\033[0m\n".format('ID','Tool','Time', 'CPU_load')
  173. for k, v in dd.iteritems():
  174. out += "{:<8} {:<17} {:<10} {:<10}\n".format(k, dd[k][0], dd[k][1][0], dd[k][1][1])
  175. print out
  176. #return out
  177. def showUsage(self):
  178. #print zonalBenchmark.__init__.__doc__
  179. print u'''Usage:
  180. python zonalStatBenchmark [tools] [input raster] [input mask / shape]
  181. [tools] expects a chain of numbers, seperated by a hyphen;
  182. e.g. 1-2-6 for selecting oft-stat, pktools and SAGA
  183. [input raster] the raster file for input
  184. [input mask / shape] the raster or shape mask; be aware: some
  185. tools allow only raster input, some only
  186. vector input - so save a vector file and
  187. its rasterization with the same name in
  188. the same dir
  189. The following tools are available:
  190. 1. oft-stat
  191. 2. pktools
  192. 3̶.̶ ̶P̶y̶t̶h̶o̶n̶/̶G̶D̶A̶L̶ ̶(̶p̶e̶r̶r̶y̶g̶e̶o̶ ̶-̶ ̶g̶i̶s̶t̶)̶
  193. 4̶.̶ ̶G̶R̶A̶S̶S̶
  194. 5̶.̶ ̶R̶
  195. 6. SAGA
  196. '''
  197. # userArgs[0] = selected providers
  198. # userArgs[1] = input TIF
  199. # userArgs[2] = masking shape
  200. # userArgs[3] = number of repetitions
  201. args = []
  202. args.append(sys.argv)
  203. userArgs = zonalBenchmark().getArguments(args)
  204. commands = zonalBenchmark().setCommands(userArgs[0], userArgs[1], userArgs[2])
  205. timeOutput = zonalBenchmark().executeCommands(commands, userArgs[3])
  206. zonalBenchmark().showOutput(timeOutput, userArgs[0])