r3.in.xyz.py 9.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321
  1. #!/usr/bin/env python3
  2. ############################################################################
  3. #
  4. # MODULE: r3.in.xyz
  5. # AUTHOR: M. Hamish Bowman, Dunedin, New Zealand
  6. # PURPOSE: Run r.in.xyz in a loop for various z-levels and construct
  7. # a 3D raster. Unlike r.in.xyz, reading from stdin and z-scaling
  8. # won't work.
  9. #
  10. # COPYRIGHT: (c) 2011-2012 Hamish Bowman, and the GRASS Development Team
  11. # Port of r3.in.xyz(.sh) for GRASS 6.4.
  12. # This program is free software under the GNU General Public
  13. # License (>=v2). Read the file COPYING that comes with GRASS
  14. # for details.
  15. #
  16. # This program is distributed in the hope that it will be useful,
  17. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  19. # GNU General Public License for more details.
  20. #
  21. #############################################################################
  22. # %Module
  23. # % description: Create a 3D raster map from an assemblage of many coordinates using univariate statistics
  24. # % keyword: raster3d
  25. # % keyword: import
  26. # % keyword: voxel
  27. # % keyword: LIDAR
  28. # % keyword: statistics
  29. # % keyword: conversion
  30. # % keyword: aggregation
  31. # % keyword: binning
  32. # %End
  33. # %Flag
  34. # % key: s
  35. # % description: Scan data file for extent then exit
  36. # %End
  37. # %Flag
  38. # % key: g
  39. # % description: In scan mode, print using shell script style
  40. # %End
  41. # %Flag
  42. # % key: i
  43. # % description: Ignore broken lines
  44. # %End
  45. # %Option G_OPT_F_INPUT
  46. # % required: yes
  47. # % description: ASCII file containing input data
  48. # %End
  49. # %Option
  50. # % key: output
  51. # % type: string
  52. # % required: yes
  53. # % multiple: no
  54. # % key_desc: name
  55. # % description: Name for output raster map
  56. # % gisprompt: new,grid3,3d-raster
  57. # %End
  58. # %Option
  59. # % key: method
  60. # % type: string
  61. # % required: no
  62. # % multiple: no
  63. # % options: n,min,max,range,sum,mean,stddev,variance,coeff_var,median,percentile,skewness,trimmean
  64. # % description: Statistic to use for raster values
  65. # % answer: mean
  66. # % guisection: Statistic
  67. # %End
  68. # %Option
  69. # % key: type
  70. # % type: string
  71. # % required: no
  72. # % multiple: no
  73. # % options: float,double
  74. # % description: Storage type for resultant raster map
  75. # % answer: float
  76. # %End
  77. # %Option G_OPT_F_SEP
  78. # % guisection: Input
  79. # %End
  80. # %Option
  81. # % key: x
  82. # % type: integer
  83. # % required: no
  84. # % multiple: no
  85. # % description: Column number of x coordinates in input file (first column is 1)
  86. # % answer: 1
  87. # % guisection: Input
  88. # %End
  89. # %Option
  90. # % key: y
  91. # % type: integer
  92. # % required: no
  93. # % multiple: no
  94. # % description: Column number of y coordinates in input file
  95. # % answer: 2
  96. # % guisection: Input
  97. # %End
  98. # %Option
  99. # % key: z
  100. # % type: integer
  101. # % required: no
  102. # % multiple: no
  103. # % description: Column number of z coordinates in input file
  104. # % answer: 3
  105. # % guisection: Input
  106. # %End
  107. # %Option
  108. # % key: value_column
  109. # % type: integer
  110. # % required: no
  111. # % multiple: no
  112. # % label: Column number of data values in input file
  113. # % description: If not given or set to 0, the data points' z-values are used
  114. # % answer: 0
  115. # % guisection: Input
  116. # %End
  117. # %Option
  118. # % key: vrange
  119. # % type: double
  120. # % required: no
  121. # % key_desc: min,max
  122. # % description: Filter range for value column data (min,max)
  123. # %End
  124. # %option
  125. # % key: vscale
  126. # % type: double
  127. # % required: no
  128. # % multiple: no
  129. # % description: Scaling factor to apply to value column data
  130. # % answer: 1.0
  131. # %end
  132. # %Option
  133. # % key: percent
  134. # % type: integer
  135. # % required: no
  136. # % multiple: no
  137. # % options: 1-100
  138. # % description: Percent of map to keep in memory
  139. # % answer: 100
  140. # %End
  141. # %Option
  142. # % key: pth
  143. # % type: integer
  144. # % required: no
  145. # % multiple: no
  146. # % options: 1-100
  147. # % description: pth percentile of the values
  148. # % guisection: Statistic
  149. # %End
  150. # %Option
  151. # % key: trim
  152. # % type: double
  153. # % required: no
  154. # % multiple: no
  155. # % options: 0-50
  156. # % description: Discard <trim> percent of the smallest and <trim> percent of the largest observations
  157. # % guisection: Statistic
  158. # %End
  159. # %Option
  160. # % key: workers
  161. # % type: integer
  162. # % required: no
  163. # % multiple: no
  164. # % options: 1-256
  165. # % answer: 1
  166. # % description: Number of parallel processes to launch
  167. # %End
  168. import sys
  169. import os
  170. import atexit
  171. from grass.script import core as grass
  172. from grass.exceptions import CalledModuleError
  173. def cleanup():
  174. grass.run_command('g.remove', flags='f',
  175. type="rast", pattern='tmp.r3xyz.%d.*' % os.getpid(),
  176. quiet=True)
  177. def main():
  178. infile = options['input']
  179. output = options['output']
  180. method = options['method']
  181. dtype = options['type']
  182. fs = options['separator']
  183. x = options['x']
  184. y = options['y']
  185. z = options['z']
  186. value_column = options['value_column']
  187. vrange = options['vrange']
  188. vscale = options['vscale']
  189. percent = options['percent']
  190. pth = options['pth']
  191. trim = options['trim']
  192. workers = int(options['workers'])
  193. scan_only = flags['s']
  194. shell_style = flags['g']
  195. ignore_broken = flags['i']
  196. if workers == 1 and "WORKERS" in os.environ:
  197. workers = int(os.environ["WORKERS"])
  198. if not os.path.exists(infile):
  199. grass.fatal(_("Unable to read input file <%s>") % infile)
  200. addl_opts = {}
  201. if pth:
  202. addl_opts['pth'] = '%s' % pth
  203. if trim:
  204. addl_opts['trim'] = '%s' % trim
  205. if value_column:
  206. addl_opts['value_column'] = '%s' % value_column
  207. if vrange:
  208. addl_opts['vrange'] = '%s' % vrange
  209. if vscale:
  210. addl_opts['vscale'] = '%s' % vscale
  211. if ignore_broken:
  212. addl_opts['flags'] = 'i'
  213. if scan_only or shell_style:
  214. if shell_style:
  215. doShell = 'g'
  216. else:
  217. doShell = ''
  218. grass.run_command('r.in.xyz', flags='s' + doShell, input=infile,
  219. output='dummy', sep=fs, x=x, y=y, z=z,
  220. **addl_opts)
  221. sys.exit()
  222. if dtype == 'float':
  223. data_type = 'FCELL'
  224. else:
  225. data_type = 'DCELL'
  226. region = grass.region(region3d=True)
  227. if region['nsres'] != region['nsres3'] or region['ewres'] != region['ewres3']:
  228. grass.run_command('g.region', flags='3p')
  229. grass.fatal(_("The 2D and 3D region settings are different. Can not continue."))
  230. grass.verbose(_("Region bottom=%.15g top=%.15g vertical_cell_res=%.15g (%d depths)")
  231. % (region['b'], region['t'], region['tbres'], region['depths']))
  232. grass.verbose(_("Creating slices ..."))
  233. # to avoid a point which falls exactly on a top bound from being
  234. # considered twice, we shrink the
  235. # For the top slice we keep it though, as someone scanning the bounds
  236. # may have set the bounds exactly to the data extent (a bad idea, but
  237. # it happens..)
  238. eps = 1.0e-15
  239. # if there are thousands of depths hopefully this array doesn't get too
  240. # large and so we don't have to worry much about storing/looping through
  241. # all the finished process infos.
  242. proc = {}
  243. pout = {}
  244. depths = list(range(1, 1 + region['depths']))
  245. for i in depths:
  246. tmp_layer_name = 'tmp.r3xyz.%d.%s' % (os.getpid(), '%05d' % i)
  247. zrange_min = region['b'] + (region['tbres'] * (i - 1))
  248. if i < region['depths']:
  249. zrange_max = region['b'] + (region['tbres'] * i) - eps
  250. else:
  251. zrange_max = region['b'] + (region['tbres'] * i)
  252. # spawn depth layer import job in the background
  253. #grass.debug("slice %d, <%s> %% %d" % (band, image[band], band % workers))
  254. grass.message(_("Processing horizontal slice %d of %d [%.15g,%.15g) ...")
  255. % (i, region['depths'], zrange_min, zrange_max))
  256. proc[i] = grass.start_command('r.in.xyz', input=infile, output=tmp_layer_name,
  257. sep=fs, method=method, x=x, y=y, z=z,
  258. percent=percent, type=data_type,
  259. zrange='%.15g,%.15g' % (zrange_min, zrange_max),
  260. **addl_opts)
  261. grass.debug("i=%d, %%=%d (workers=%d)" % (i, i % workers, workers))
  262. # print sys.getsizeof(proc) # sizeof(proc array) [not so big]
  263. if i % workers == 0:
  264. # wait for the ones launched so far to finish
  265. for p_i in depths[:i]:
  266. pout[p_i] = proc[p_i].communicate()[0]
  267. if proc[p_i].wait() != 0:
  268. grass.fatal(_("Trouble importing data. Aborting."))
  269. # wait for jSobs to finish, collect any stray output
  270. for i in depths:
  271. pout[i] = proc[i].communicate()[0]
  272. if proc[i].wait() != 0:
  273. grass.fatal(_("Trouble importing data. Aborting."))
  274. del proc
  275. grass.verbose(_("Assembling 3D cube ..."))
  276. # input order: lower most strata first
  277. slices = grass.read_command('g.list', type='raster', sep=',',
  278. pattern='tmp.r3xyz.%d.*' % os.getpid()).rstrip(os.linesep)
  279. grass.debug(slices)
  280. try:
  281. grass.run_command('r.to.rast3', input=slices, output=output)
  282. except CalledModuleError:
  283. grass.message(_("Done. 3D raster map <%s> created.") % output)
  284. if __name__ == "__main__":
  285. options, flags = grass.parser()
  286. atexit.register(cleanup)
  287. main()