[GRASS-SVN] r52143 - in grass/trunk/scripts: . r3.in.xyz
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 19 17:42:25 PDT 2012
Author: hamish
Date: 2012-06-19 17:42:25 -0700 (Tue, 19 Jun 2012)
New Revision: 52143
Added:
grass/trunk/scripts/r3.in.xyz/
grass/trunk/scripts/r3.in.xyz/Makefile
grass/trunk/scripts/r3.in.xyz/r3.in.xyz.html
grass/trunk/scripts/r3.in.xyz/r3.in.xyz.py
Log:
initial port of r3.in.xyz from grass6 addons; raw and completely untested
Added: grass/trunk/scripts/r3.in.xyz/Makefile
===================================================================
--- grass/trunk/scripts/r3.in.xyz/Makefile (rev 0)
+++ grass/trunk/scripts/r3.in.xyz/Makefile 2012-06-20 00:42:25 UTC (rev 52143)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r3.in.xyz
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass/trunk/scripts/r3.in.xyz/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass/trunk/scripts/r3.in.xyz/r3.in.xyz.html
===================================================================
--- grass/trunk/scripts/r3.in.xyz/r3.in.xyz.html (rev 0)
+++ grass/trunk/scripts/r3.in.xyz/r3.in.xyz.html 2012-06-20 00:42:25 UTC (rev 52143)
@@ -0,0 +1,75 @@
+<h2>DESCRIPTION</h2>
+
+<em>r3.in.xyz</em> imports sparse XYZ data from an ASCII file into
+a 3D raster map (voxels). It does this by running the <em>r.in.xyz</em>
+module multiple times for different z-ranges and then assembling the
+slices with <em>r.to.rast3</em>.
+<p>
+See the <a href="r.in.xyz.html">r.in.xyz</a> help page for general
+parameter usage and tips.
+<p>
+The map is created using the rows, columns, and depths set by
+current region settings. Be sure to check and adjust these with
+the <em>g.region</em> module before performing the import.
+
+
+<h2>NOTES</h2>
+
+The 2D and 3D horizontal region resolutions must match. See the
+EXAMPLES section below.
+<p>
+Unlike <em>r.in.xyz</em>, reading from stdin and z-scaling are not
+possible.
+<p>
+To enable parallel processing support, set the <b>workers=</b> option
+to match the number of CPUs or CPU-cores avaiable on your system.
+<p>
+Points falling exactly on a vertical bound will belong to the depth
+band below them, except for points exactly on the top bound, which will
+belong to the top-most slice.
+<p>
+The script is expected to be nearly as efficient as if it was fully
+written in C.
+
+
+<h2>EXAMPLE</h2>
+
+Using the Serpent Mound dataset. (see the
+ <a href="http://grass.osgeo.org/wiki/LIDAR">GRASS LiDAR wiki page</a>)
+
+<div class="code"><pre>
+ #scan dataset for extent:
+ r3.in.xyz -s in=Serpent_Mound_Model_LAS_Data.txt out=dummy \
+ x=1 y=2 z=3 fs=space
+
+ # set the 2D and 3D regions:
+ g.region n=4323641.57 s=4320942.61 w=289020.90 e=290106.02 res=1 -a
+ g.region b=166 t=216 res3=1 tbres=5 -3 -p
+
+ r3.in.xyz in=Serpent_Mound_Model_LAS_Data.txt out=serpent3D \
+ method=mean x=1 y=2 z=3 fs=space type=float
+</pre></div>
+
+
+<h2>BUGS</h2>
+
+<em>r.to.rast3</em> always creates a <tt>double</tt> output map
+regardless of input.
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a><br>
+<a href="r.in.xyz.html">r.in.xyz</a><br>
+<a href="r.to.rast3.html">r.to.rast3</a><br>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Hamish Bowman<br>
+<i>Dunedin, New Zealand</i>
+
+<p>
+<i>Last changed: $Date$</i>
Property changes on: grass/trunk/scripts/r3.in.xyz/r3.in.xyz.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/scripts/r3.in.xyz/r3.in.xyz.py
===================================================================
--- grass/trunk/scripts/r3.in.xyz/r3.in.xyz.py (rev 0)
+++ grass/trunk/scripts/r3.in.xyz/r3.in.xyz.py 2012-06-20 00:42:25 UTC (rev 52143)
@@ -0,0 +1,283 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: r3.in.xyz
+# AUTHOR: M. Hamish Bowman, Dunedin, New Zealand
+# PURPOSE: Run r.in.xyz in a loop for various z-levels and construct
+# a 3D raster. Unlike r.in.xyz, reading from stdin and z-scaling
+# won't work.
+#
+# COPYRIGHT: (c) 2011-2012 Hamish Bowman, and the GRASS Development Team
+# Port of r3.in.xyz(.sh) for GRASS 6.4.
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+#############################################################################
+
+#%Module
+#% description: Create a 3D raster map from an assemblage of many coordinates using univariate statistics
+#% keywords: raster3d, import, voxel, LiDAR
+#%End
+#%Flag
+#% key: s
+#% description: Scan data file for extent then exit
+#%End
+#%Flag
+#% key: g
+#% description: In scan mode, print using shell script style
+#%End
+#%Flag
+#% key: i
+#% description: Ignore broken lines
+#%End
+#%Option
+#% key: input
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: name
+#% description: ASCII file containing input data
+#% gisprompt: old_file,file,input
+#%End
+#%Option
+#% key: output
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: name
+#% description: Name for output raster map
+#% gisprompt: new,grid3,3d-raster
+#%End
+#%Option
+#% key: method
+#% type: string
+#% required: no
+#% multiple: no
+#% options: n,min,max,range,sum,mean,stddev,variance,coeff_var,median,percentile,skewness,trimmean
+#% description: Statistic to use for raster values
+#% answer: mean
+#% guisection: Statistic
+#%End
+#%Option
+#% key: type
+#% type: string
+#% required: no
+#% multiple: no
+#% options: float,double
+#% description: Storage type for resultant raster map
+#% answer: float
+#%End
+#%Option
+#% key: fs
+#% type: string
+#% required: no
+#% multiple: no
+#% key_desc: character
+#% description: Field separator
+#% answer: |
+#% guisection: Input
+#%End
+#%Option
+#% key: x
+#% type: integer
+#% required: no
+#% multiple: no
+#% description: Column number of x coordinates in input file (first column is 1)
+#% answer: 1
+#% guisection: Input
+#%End
+#%Option
+#% key: y
+#% type: integer
+#% required: no
+#% multiple: no
+#% description: Column number of y coordinates in input file
+#% answer: 2
+#% guisection: Input
+#%End
+#%Option
+#% key: z
+#% type: integer
+#% required: no
+#% multiple: no
+#% description: Column number of data values in input file
+#% answer: 3
+#% guisection: Input
+#%End
+#%Option
+#% key: percent
+#% type: integer
+#% required: no
+#% multiple: no
+#% options: 1-100
+#% description: Percent of map to keep in memory
+#% answer: 100
+#%End
+#%Option
+#% key: pth
+#% type: integer
+#% required: no
+#% multiple: no
+#% options: 1-100
+#% description: pth percentile of the values
+#% guisection: Statistic
+#%End
+#%Option
+#% key: trim
+#% type: double
+#% required: no
+#% multiple: no
+#% options: 0-50
+#% description: Discard <trim> percent of the smallest and <trim> percent of the largest observations
+#% guisection: Statistic
+#%End
+#%Option
+#% key: workers
+#% type: integer
+#% required: no
+#% multiple: no
+#% options: 1-256
+#% answer: 1
+#% description: Number of parallel processes to launch
+#%End
+
+
+import sys
+import os
+from grass.script import core as grass
+
+
+def cleanup():
+ # FIXME for pid and regex
+ grass.run_command('g.mremove', flags = 'f', rast = 'tmp.r3xyz.$$.*', quiet = True)
+
+
+def main():
+ infile = options['input']
+ output = options['output']
+ method = options['method']
+ dtype = options['type']
+ fs = options['fs']
+ x = options['x']
+ y = options['y']
+ z = options['z']
+ percent = options['percent']
+ pth = options['pth']
+ trim = options['trim']
+ workers = options['workers']
+ scan_only = flags['s']
+ shell_style = flags['g']
+ ignore_broken = flags['i']
+
+ if workers is 1 and "WORKERS" in os.environ:
+ workers = int(os.environ["WORKERS"])
+
+ if not os.path.exists(infile):
+ grass.fatal(_("Unable to read input file <%s>") % infile)
+
+ if scan_only or shell_style:
+ if shell_style:
+ doShell = 'g'
+ else:
+ doShell = ''
+ grass.run_command('r.in.xyz', flags = 's' + doShell, input = infile,
+ output = 'dummy', fs = fs, x = x, y = y, z = z)
+ sys.exit()
+
+ addl_opts = []
+ if pth:
+ addl_opts.append("pth = '%s'" % pth)
+ if trim:
+ addl_opts.append("trim = '%s'" % trim)
+ if ignore_broken:
+ addl_opts.append("flags = 'i'")
+
+ if dtype is 'float':
+ data_type = 'FCELL'
+ else:
+ data_type = 'DCELL'
+
+ region = grass.region(region3d = True)
+
+ if region['nsres'] != region['nsres3'] or region['ewres'] != region['ewres3']:
+ grass.run_command('g.region', flags = '3p')
+ grass.error(_("The 2D and 3D region settings are different. Can not continue."))
+
+ grass.verbose(_("Region bottom=%.15g top=%.15g vertical_cell_res=%.15g (depths %d)"),
+ % (region['b'], region['t'], region['tbres'], region['depths']))
+
+ grass.verbose(_("Creating slices ..."))
+
+ # to avoid a point which falls exactly on a top bound from being
+ # considered twice, we shrink the
+ # For the top slice we keep it though, as someone scanning the bounds
+ # may have set the bounds exactly to the data extent (a bad idea, but
+ # it happens..)
+ eps = 1.0e-15
+
+ # if there are thousands of depths hopefully this array doesn't get too large
+ # and so we don't have to worry much about storing/looping through the finished process infos.
+ proc = {}
+
+ for i in range(1, 1 + depths):
+ i_str = '%05d' % i
+ tmp_layer_name = 'tmp.r3xyz.%d.%s' % (os.getpid(), i_str)
+
+ zrange_min = region['b'] + (region['tbres'] * (i-1))
+
+ if i < region['depths']:
+ zrange_max = region['b'] + (region['tbres'] * i) - eps
+ else:
+ zrange_max = region['b'] + (region['tbres'] * i)
+
+ # spawn depth layerimport job in the background
+ #grass.debug("slice %d, <%s> %% %d" % (band, image[band], band % workers))
+ grass.message(_("Processing horizontal slice %d of %d [%.15g,%.15g) ..."),
+ % (i, depths, zrange_min, zrange_max))
+
+ proc[i] = grass.start_command('r.in.xyz', input = infile, output = tmp_layer_name,
+ fs = fs, method = method, x = x, y = y, z = z,
+ percent = percent, type = dtype,
+ zrange = '%.15g,%.15g' % (zrange_min, zrange_max),
+ addl_opts)
+
+ if i % workers is 0:
+ # wait for the ones launched so far to finish
+ for p_i in range(1, 1 + depths)[:i]:
+ if not proc[p_i].stdout.closed:
+ pout[p_i] = proc[p_i].communicate()[0]
+ if proc[p_i].wait() is not 0:
+ grass.error(_("Trouble importing data. Aborting."))
+
+
+ # wait for jobs to finish, collect any stray output
+ for p_i in range(1, 1 + depths)[:i]:
+ if not proc[p_i].stdout.closed:
+ pout[p_i] = proc[p_i].communicate()[0]
+ if proc[p_i].wait() is not 0:
+ grass.error(_("Trouble importing data. Aborting."))
+
+ del proc
+
+ grass.verbose(_("Assembling 3D cube ..."))
+
+ #input order: lower most strata first
+ slices = grass.read_command('g.mlist', type = 'rast', sep = ',',
+ pattern = 'tmp.r3xyz.%d.*' % os.getpid())
+
+ if grass.run_command('r.to.rast3', input = slices, output = output) is 0:
+ grass.message(_("Done. 3D raster map <%s> created.") % output)
+
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ main()
+
+
Property changes on: grass/trunk/scripts/r3.in.xyz/r3.in.xyz.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list