make r3.in.xyz opertational, add some notes about future improvements

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@52148 15284696-431f-4ddb-bdfa-cd5b030d7da7

scripts/Makefile

@@ -40,6 +40,7 @@ SUBDIRS = \
 	r.shaded.relief \
 	r.tileset \
 	r.unpack \
+	r3.in.xyz \
 	v.build.all \
 	v.centroids \
 	v.convert.all \

scripts/r3.in.xyz/r3.in.xyz.html

@@ -51,6 +51,30 @@ Using the Serpent Mound dataset. (see the
 </pre></div>
 
 
+<h2>TODO</h2>
+
+Currently this module stores only the z-value into the voxel,
+which is perhaps a little monochromatic. The 'n' statistic
+and simple presence,absence may still be of interest though.
+<p>
+Possible ideas for adding full <tt>x,y,z,value</tt> support
+to the module include piping the input file through a z-range
+gate filter in Python while simultaneously setting the
+<em>r.in.xyz</em> module's <tt>z=</tt> option to be the
+column of the value to be stored. This will be slower than a
+dedicated <tt>r.in.xyzv</tt> module would be, but less memory
+intensive and less complicated to code (fewer bugs). The question
+is how much slower would it be, and if that is acceptable. It
+seems like a reasonable approach for experimenting with a
+working prototype though.
+<p>
+Another idea (which is likely the way I'll go) is to add a
+new extra-column filter option to <em>r.in.xyz</em>. The
+data can then be filtered by z-range while simultaneously
+storing the data of another (possibly gate-ranged itself)
+column.
+
+
 <h2>BUGS</h2>
 
 <em>r.to.rast3</em> always creates a <tt>double</tt> output map

scripts/r3.in.xyz/r3.in.xyz.py

@@ -54,7 +54,7 @@
 #% description: Name for output raster map
 #% gisprompt: new,grid3,3d-raster
 #%End
-#%Option      
+#%Option
 #% key: method
 #% type: string
 #% required: no
@@ -63,54 +63,72 @@
 #% description: Statistic to use for raster values
 #% answer: mean
 #% guisection: Statistic
-#%End         
-#%Option      
-#% key: type  
+#%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    
+#%End
+#%Option
+#% key: fs
 #% type: string
 #% required: no
 #% multiple: no
 #% key_desc: character
 #% description: Field separator
-#% answer: |  
+#% answer: |
 #% guisection: Input
-#%End         
-#%Option      
-#% key: x     
+#%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  
+#% answer: 1
 #% guisection: Input
-#%End         
-#%Option      
-#% key: y     
+#%End
+#%Option
+#% key: y
 #% type: integer
 #% required: no
 #% multiple: no
 #% description: Column number of y coordinates in input file
-#% answer: 2  
+#% answer: 2
 #% guisection: Input
-#%End         
-#%Option      
-#% key: z     
+#%End
+#%Option
+#% key: z
+#% type: integer
+#% required: no
+#% multiple: no
+#% description: Column number of z coordinates in input file
+#% answer: 3
+#% guisection: Input
+#%End
+#%Option
+#% key: value
 #% type: integer
 #% required: no
 #% multiple: no
-#% description: Column number of data values in input file
-#% answer: 3  
+#% label: Column number of data values in input file
+#% description: If not given or set to 0, the data points' z-values are used
+#% answer: 0
 #% guisection: Input
-#%End         
-#%Option      
+#%End
+#%Option
+#% key: vrange
+#% type: double
+#% required: no
+#% key_desc: min,max
+#% description: Filter range for value data (min,max)
+#% guisection: Input
+#%End
+#%Option
 #% key: percent
 #% type: integer
 #% required: no
@@ -118,26 +136,26 @@
 #% options: 1-100
 #% description: Percent of map to keep in memory
 #% answer: 100
-#%End         
-#%Option      
-#% key: pth   
+#%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  
+#%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      
+#%End
+#%Option
 #% key: workers
 #% type: integer
 #% required: no
@@ -145,16 +163,18 @@
 #% options: 1-256
 #% answer: 1
 #% description: Number of parallel processes to launch
-#%End         
+#%End
 
 
 import sys
 import os
+import atexit
 from grass.script import core as grass
 
 
 def cleanup():
-    grass.run_command('g.mremove', flags = 'f', rast = 'tmp.r3xyz.%d.*' % os.getpid(),
+    grass.run_command('g.mremove', flags = 'f',
+		      rast = 'tmp.r3xyz.%d.*' % os.getpid(),
 		      quiet = True)
 
 
@@ -165,12 +185,14 @@ def main():
     dtype = options['type']
     fs = options['fs']
     x = options['x']
-    y = options['y']   
+    y = options['y']
     z = options['z']
+    value = int(options['value'])
+    vrange = options['vrange']
     percent = options['percent']
     pth = options['pth']
     trim = options['trim']
-    workers = options['workers']
+    workers = int(options['workers'])
     scan_only = flags['s']
     shell_style = flags['g']
     ignore_broken = flags['i']
@@ -181,6 +203,9 @@ def main():
     if not os.path.exists(infile):
 	grass.fatal(_("Unable to read input file <%s>") % infile)
 
+    if value is not 0:
+	grass.fatal(_("This functionality is not yet operational."))
+
     if scan_only or shell_style:
         if shell_style:
 	    doShell = 'g'
@@ -190,13 +215,13 @@ def main():
 			  output = 'dummy', fs = fs, x = x, y = y, z = z)
 	sys.exit()
 
-    addl_opts = []
+    addl_opts = {}
     if pth:
-        addl_opts.append("pth = '%s'" % pth)
+        addl_opts['pth'] = '%s' % pth
     if trim:
-        addl_opts.append("trim = '%s'" % trim)
+        addl_opts['trim'] = '%s' % trim
     if ignore_broken:
-        addl_opts.append("flags = 'i'")
+        addl_opts['flags'] = 'i'
 
     if dtype is 'float':
        data_type = 'FCELL'
@@ -209,7 +234,7 @@ def main():
         grass.run_command('g.region', flags = '3p')
         grass.fatal(_("The 2D and 3D region settings are different. Can not continue."))
 
-    grass.verbose(_("Region bottom=%.15g  top=%.15g  vertical_cell_res=%.15g  (depths %d)")
+    grass.verbose(_("Region bottom=%.15g  top=%.15g  vertical_cell_res=%.15g  (%d depths)")
        % (region['b'], region['t'], region['tbres'], region['depths']))
 
     grass.verbose(_("Creating slices ..."))
@@ -221,11 +246,15 @@ def main():
     # 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.
+    # 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
+    # all the finished process infos.
     proc = {}
+    pout = {}
+
+    depths = range(1, 1 + region['depths'])
 
-    for i in range(1, 1 + depths):
+    for i in depths:
 	tmp_layer_name = 'tmp.r3xyz.%d.%s' % (os.getpid(), '%05d' % i)
 
         zrange_min = region['b'] + (region['tbres'] * (i-1))
@@ -237,28 +266,33 @@ def main():
 
         # 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))
+        grass.message(_("Processing horizontal slice %d of %d [%.15g,%.15g) ...")
+		        % (i, region['depths'], zrange_min, zrange_max))
 
+	#if not value:
 	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 = data_type,
 				      zrange = '%.15g,%.15g' % (zrange_min, zrange_max),
-				      addl_opts)
+				      **addl_opts)
+	#else:
+	#    use new alternate gate-range filter option in r.in.xyz
+
+
+	grass.debug("i=%d, %%=%d  (workers=%d)" % (i, i % workers, workers))
+	#print sys.getsizeof(proc)  # sizeof(proc array)  [not so big]
 
 	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]
+	    for p_i in depths[:i]:
+		pout[p_i] = proc[p_i].communicate()[0]
 		if proc[p_i].wait() is not 0:
 		    grass.fatal(_("Trouble importing data. Aborting."))
 
 
-    # wait for jobs to finish, collect any stray output
-    for i in range(1, 1 + depths):
-	if not proc[i].stdout.closed:
-	    pout[i] = proc[i].communicate()[0]
+    # wait for jSobs to finish, collect any stray output
+    for i in depths:
+	pout[i] = proc[i].communicate()[0]
 	if proc[i].wait() is not 0:
 	    grass.fatal(_("Trouble importing data. Aborting."))
 
@@ -267,8 +301,9 @@ def main():
     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())
+    slices = grass.read_command('g.mlist', type = 'rast', fs = ',',
+				pattern = 'tmp.r3xyz.%d.*' % os.getpid()).rstrip('\n')
+    grass.debug(slices)
 
     if grass.run_command('r.to.rast3', input = slices, output = output) is 0:
 	grass.message(_("Done. 3D raster map <%s> created.") % output)