[GRASS-SVN] r49592 - in grass-addons/grass6/raster3d: . r3.in.xyz
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 7 01:52:45 EST 2011
Author: hamish
Date: 2011-12-06 22:52:44 -0800 (Tue, 06 Dec 2011)
New Revision: 49592
Added:
grass-addons/grass6/raster3d/r3.in.xyz/
grass-addons/grass6/raster3d/r3.in.xyz/Makefile
grass-addons/grass6/raster3d/r3.in.xyz/description.html
grass-addons/grass6/raster3d/r3.in.xyz/r3.in.xyz
Log:
like r.in.xyz, but imports into a cube. parallel processing is supported
Copied: grass-addons/grass6/raster3d/r3.in.xyz/Makefile (from rev 49589, grass-addons/grass6/raster/r.in.xyz.auto/Makefile)
===================================================================
--- grass-addons/grass6/raster3d/r3.in.xyz/Makefile (rev 0)
+++ grass-addons/grass6/raster3d/r3.in.xyz/Makefile 2011-12-07 06:52:44 UTC (rev 49592)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r3.in.xyz
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass6/raster3d/r3.in.xyz/description.html
===================================================================
--- grass-addons/grass6/raster3d/r3.in.xyz/description.html (rev 0)
+++ grass-addons/grass6/raster3d/r3.in.xyz/description.html 2011-12-07 06:52:44 UTC (rev 49592)
@@ -0,0 +1,71 @@
+<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>
+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-addons/grass6/raster3d/r3.in.xyz/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Copied: grass-addons/grass6/raster3d/r3.in.xyz/r3.in.xyz (from rev 49589, grass-addons/grass6/raster/r.in.xyz.auto/r.in.xyz.auto)
===================================================================
--- grass-addons/grass6/raster3d/r3.in.xyz/r3.in.xyz (rev 0)
+++ grass-addons/grass6/raster3d/r3.in.xyz/r3.in.xyz 2011-12-07 06:52:44 UTC (rev 49592)
@@ -0,0 +1,303 @@
+#!/bin/sh
+############################################################################
+#
+# 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 Hamish Bowman, and the GRASS Development Team
+# 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
+
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." 1>&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+
+if [ ! -e "$GIS_OPT_INPUT" ] ; then
+ echo "Input file not found." 1>&2
+ exit 1
+fi
+
+
+if [ "$GIS_FLAG_S" -eq 1 ] || [ "$GIS_FLAG_G" -eq 1 ] ; then
+ if [ "$GIS_FLAG_G" -eq 1 ] ; then
+ SCAN_FLAG="-g"
+ else
+ SCAN_FLAG=""
+ fi
+
+ r.in.xyz input="$GIS_OPT_INPUT" output=dummy fs="$GIS_OPT_FS" \
+ x="$GIS_OPT_X" y="$GIS_OPT_Y" z="$GIS_OPT_Z" -s $SCAN_FLAG
+
+ exit
+fi
+
+#### check if we have awk
+if [ ! -x "`which awk`" ] ; then
+ g.message -e "awk required, please install awk or gawk first"
+ exit 1
+fi
+
+# set environment so that awk works properly in all locales
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+
+cleanup()
+{
+ g.message -v "Cleaning up ..."
+ g.mremove -f --quiet \
+ rast=`g.mlist type=rast pattern="tmp.r3xyz.$$.*" sep=,`
+}
+trap "cleanup" 2 3 15
+
+
+# set up opts only given if there
+ADDL_OPTS=""
+if [ -n "$GIS_OPT_PTH" ] ; then
+ ADDL_OPTS="pth=\"$GIS_OPT_PTH\""
+fi
+if [ -n "$GIS_OPT_TRIM" ] ; then
+ ADDL_OPTS="$ADDL_OPTS trim=\"$GIS_OPT_TRIM\""
+fi
+if [ "$GIS_FLAG_I" -eq 1 ] ; then
+ ADDL_OPTS="$ADDL_OPTS -i"
+fi
+if [ "$GIS_OPT_TYPE" = "float" ] ; then
+ DATA_TYPE="FCELL"
+else
+ DATA_TYPE="DCELL"
+fi
+
+
+eval `g.region -3 -g`
+
+if [ "$nsres" != "$nsres3" ] || [ "$ewres" != "$ewres3" ] ; then
+ g.region -3 -p | grep res
+ g.message -e "The 2D and 3D region settings are different. Can not continue."
+ exit 1
+fi
+
+g.message -v message="Region bottom=$b top=$t \
+ vertical_cell_res=$tbres ($depths depths)"
+
+g.message -v message="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..)
+EPSILON=1.0e-15
+# init
+EXIT_CODE=0
+
+for i in `seq "$depths"` ; do
+
+ # $b + $tbres * ($i-1)
+ zrange_min=`echo "$b $tbres $i" | awk '{printf("%.15g", $1 + ($2*($3 - 1)))}'`
+
+ if [ $i -lt $depths ] ; then
+ # $b + $tbres * ($i) - $EPSILON
+ zrange_max=`echo "$b $tbres $i $EPSILON" | awk '{printf("%.15g", $1 + ($2 * $3) - $4)}'`
+ else
+ # $b + $tbres * ($i)
+ zrange_max=`echo "$b $tbres $i" | awk '{printf("%.15g", $1 + ($2 * $3))}'`
+ fi
+
+ # import it
+ CMD="r.in.xyz input=\"$GIS_OPT_INPUT\" output=\"tmp.r3xyz.$$.$i\" \
+ fs=\"$GIS_OPT_FS\" method=\"$GIS_OPT_METHOD\" \
+ x=\"$GIS_OPT_X\" y=\"$GIS_OPT_Y\" z=\"$GIS_OPT_Z\" \
+ percent=\"$GIS_OPT_PERCENT\" type=\"$DATA_TYPE\" $ADDL_OPTS \
+ zrange=\"$zrange_min,$zrange_max\""
+
+ g.message -d debug=1 message="$CMD"
+
+ # parallel launching: could use GNU Parallel (an advanced form of xargs), but
+ # it may not be available. so we use a more generic approach
+ # see http://www.gnu.org/software/parallel/
+ # and http://grass.osgeo.org/wiki/OpenMP#Alternatives
+
+ # poor man's multi-threading for a multi-core CPU
+ MODULUS=`echo "$i 1 $GIS_OPT_WORKERS" | awk '{print $1 % ($2 * $3)}'`
+
+ if [ "$MODULUS" = "0" -o "$i" -eq "$depths" ] ; then
+ # stall to let the background jobs finish
+ g.message \
+ "Processing horizontal slice $i of $depths [$zrange_min,$zrange_max) ..."
+ eval $CMD
+ EXIT_CODE=$?
+ sleep 2
+ else
+ g.message \
+ "Launching horizontal slice $i of $depths [$zrange_min,$zrange_max) in parallel ..."
+ eval $CMD &
+ fi
+
+ if [ $? -ne 0 -o "$EXIT_CODE" -ne 0 ] ; then
+ #killall r.in.xyz
+ cleanup
+ g.message -e "Trouble importing data. Aborting."
+ exit 1
+ fi
+done
+
+
+g.message -v message="Assembling 3D cube ..."
+
+r.to.rast3 input=`g.mlist type=rast pattern="tmp.r3xyz.$$.*" sep=,` \
+ output="$GIS_OPT_OUTPUT"
+
+
+
+if [ $? -eq 0 ] ; then
+ g.message "Done. 3D raster map <$GIS_OPT_OUTPUT> created."
+fi
+
+cleanup
More information about the grass-commit
mailing list