[GRASS-SVN] r30227 - in grass-addons: . raster raster/r.boxcount
raster/r.boxcount.sh vector vector/v.strahler
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Feb 18 05:12:36 EST 2008
Author: florian
Date: 2008-02-18 05:12:35 -0500 (Mon, 18 Feb 2008)
New Revision: 30227
Added:
grass-addons/raster/r.boxcount.sh/
grass-addons/raster/r.boxcount.sh/Makefile
grass-addons/raster/r.boxcount.sh/description.html
grass-addons/raster/r.boxcount.sh/r.boxcount.sh
grass-addons/raster/r.boxcount/
grass-addons/raster/r.boxcount/Makefile
grass-addons/raster/r.boxcount/bits.c
grass-addons/raster/r.boxcount/bits.h
grass-addons/raster/r.boxcount/description.html
grass-addons/raster/r.boxcount/file.c
grass-addons/raster/r.boxcount/file.h
grass-addons/raster/r.boxcount/main.c
grass-addons/raster/r.boxcount/record.h
grass-addons/raster/r.boxcount/sort.c
grass-addons/raster/r.boxcount/sort.h
grass-addons/vector/v.strahler/
grass-addons/vector/v.strahler/Makefile
grass-addons/vector/v.strahler/description.html
grass-addons/vector/v.strahler/forest2tree.c
grass-addons/vector/v.strahler/helper.c
grass-addons/vector/v.strahler/images/
grass-addons/vector/v.strahler/main.c
grass-addons/vector/v.strahler/r.broscoe.sh
grass-addons/vector/v.strahler/r.strahler.sh
grass-addons/vector/v.strahler/strahler.c
grass-addons/vector/v.strahler/strahler.h
grass-addons/vector/v.strahler/write.c
Removed:
grass-addons/r.boxcount.sh/
grass-addons/r.boxcount/
grass-addons/raster/r.boxcount.sh/Makefile
grass-addons/raster/r.boxcount.sh/description.html
grass-addons/raster/r.boxcount.sh/r.boxcount.sh
grass-addons/raster/r.boxcount/Makefile
grass-addons/raster/r.boxcount/bits.c
grass-addons/raster/r.boxcount/bits.h
grass-addons/raster/r.boxcount/description.html
grass-addons/raster/r.boxcount/file.c
grass-addons/raster/r.boxcount/file.h
grass-addons/raster/r.boxcount/main.c
grass-addons/raster/r.boxcount/record.h
grass-addons/raster/r.boxcount/sort.c
grass-addons/raster/r.boxcount/sort.h
grass-addons/v.strahler/
grass-addons/vector/v.strahler/Makefile
grass-addons/vector/v.strahler/description.html
grass-addons/vector/v.strahler/forest2tree.c
grass-addons/vector/v.strahler/helper.c
grass-addons/vector/v.strahler/main.c
grass-addons/vector/v.strahler/r.strahler.sh
grass-addons/vector/v.strahler/strahler.c
grass-addons/vector/v.strahler/strahler.h
grass-addons/vector/v.strahler/write.c
Log:
Moved v.strahler and r.boxcount* to vector and raster subdirectories
Copied: grass-addons/raster/r.boxcount (from rev 30225, grass-addons/r.boxcount)
Deleted: grass-addons/raster/r.boxcount/Makefile
===================================================================
--- grass-addons/r.boxcount/Makefile 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,12 +0,0 @@
-# fix this relative to include/
-# or use absolute path to the GRASS source code
-MODULE_TOPDIR = ../..
-
-PGM = r.boxcount
-
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
-
-include $(MODULE_TOPDIR)/include/Make/Module.make
-
-default: cmd
Copied: grass-addons/raster/r.boxcount/Makefile (from rev 30226, grass-addons/r.boxcount/Makefile)
===================================================================
--- grass-addons/raster/r.boxcount/Makefile (rev 0)
+++ grass-addons/raster/r.boxcount/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,12 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.boxcount
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Deleted: grass-addons/raster/r.boxcount/bits.c
===================================================================
--- grass-addons/r.boxcount/bits.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/bits.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,161 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- * PURPOSE:
- * These functions manipuluate individual bits and might
- * need modification on a Big-Endian machine (e.g SPARC, 68xxx).
- * This code was developed on an 80x86.
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-
-
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <limits.h>
-#include "bits.h"
-
-
-/**************************************************************************/
-
-unsigned Bits (unsigned int a, int k, int j)
-{
- /* Return the j bits which occur k bits from
- the right (i.e. the 8th bit occurs 7 bits
- from the right) */
-
- return (a >> k) & ~(~0 << j);
-}
-
-
-/**************************************************************************/
-
-int Bits_available (char type)
-{
- unsigned long int actual, expected;
- int bits;
-
- if (type == 'l')
- {
- actual = ULONG_MAX;
- expected = 4294967295UL;
- bits = 31;
- }
-
- if (type == 'i')
- {
- actual = UINT_MAX;
- expected = 4294967295UL;
- bits = 31;
- }
-
- if (type == 's')
- {
- actual = USHRT_MAX;
- expected = 65535;
- bits = 15;
- }
-
- if (actual != expected)
- {
- fprintf (stderr, "\n *** Word length of integer not as expected *** \n");
- fflush (stderr);
- bits = 0;
- }
-
- return bits;
-}
-
-
-/**************************************************************************/
-
-unsigned int Create_chunked_bit_string (unsigned int a,
- unsigned int b, int k)
-{
- /* This function takes 2 bit strings of length k and
- concatenates them. For example if a=1101 (13) and
- b=1001 (9) then the result will be 11011001 (217)
-
- Left shift a by k positions and then add b */
-
- return ((a << k) + b);
-}
-
-
-/**************************************************************************/
-
-unsigned int Create_cyclic_bit_string (unsigned int a,
- unsigned int b, int k)
-{
- /* This function combines 2 bit strings of length k
- into one string of length 2*k such that the individual
- bits of a and b are ordered:
- k bit of a -> 2k bit of bitstr
- k bit of b -> 2k-1 bit of bitstr
- k-1 bit of a -> 2k-2 bit of bitstr
- k-1 bit of b -> 2k-3 bit of bitstr
- k-2 bit of a -> 2k-4 bit of bitstr
- ...
- k-(k-1) bit of b -> 2k-(2k-1) bit of bitstr */
-
- unsigned int bitstr;
- int bit, ls;
-
- bitstr = 0;
- ls = 2 * k;
- for (bit = k - 1; bit >= 0; bit --)
- {
- ls --;
- bitstr += Bits (a, bit, 1) << ls;
- ls --;
- bitstr += Bits (b, bit, 1) << ls;
- }
- return bitstr;
-}
-
-
-/**************************************************************************/
-
-int Bits_in_pwr_two_can_sqr_and_store_in (char type)
-{
- int bits;
-
- bits = Bits_available (type);
-
-
- /* Word which is an even no. of bits can store
- 2^x * 2^x where x = bits / 2. Word which is
- an odd no. of bits can store 2^x * 2^x where
- x = (bits + 1) / 2 */
-
- if (fmod (bits, 2) != 1)
- return (bits / 2);
- else
- return ((int) ((bits + 1) / 2));
-}
Copied: grass-addons/raster/r.boxcount/bits.c (from rev 30226, grass-addons/r.boxcount/bits.c)
===================================================================
--- grass-addons/raster/r.boxcount/bits.c (rev 0)
+++ grass-addons/raster/r.boxcount/bits.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,161 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ * PURPOSE:
+ * These functions manipuluate individual bits and might
+ * need modification on a Big-Endian machine (e.g SPARC, 68xxx).
+ * This code was developed on an 80x86.
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <limits.h>
+#include "bits.h"
+
+
+/**************************************************************************/
+
+unsigned Bits (unsigned int a, int k, int j)
+{
+ /* Return the j bits which occur k bits from
+ the right (i.e. the 8th bit occurs 7 bits
+ from the right) */
+
+ return (a >> k) & ~(~0 << j);
+}
+
+
+/**************************************************************************/
+
+int Bits_available (char type)
+{
+ unsigned long int actual, expected;
+ int bits;
+
+ if (type == 'l')
+ {
+ actual = ULONG_MAX;
+ expected = 4294967295UL;
+ bits = 31;
+ }
+
+ if (type == 'i')
+ {
+ actual = UINT_MAX;
+ expected = 4294967295UL;
+ bits = 31;
+ }
+
+ if (type == 's')
+ {
+ actual = USHRT_MAX;
+ expected = 65535;
+ bits = 15;
+ }
+
+ if (actual != expected)
+ {
+ fprintf (stderr, "\n *** Word length of integer not as expected *** \n");
+ fflush (stderr);
+ bits = 0;
+ }
+
+ return bits;
+}
+
+
+/**************************************************************************/
+
+unsigned int Create_chunked_bit_string (unsigned int a,
+ unsigned int b, int k)
+{
+ /* This function takes 2 bit strings of length k and
+ concatenates them. For example if a=1101 (13) and
+ b=1001 (9) then the result will be 11011001 (217)
+
+ Left shift a by k positions and then add b */
+
+ return ((a << k) + b);
+}
+
+
+/**************************************************************************/
+
+unsigned int Create_cyclic_bit_string (unsigned int a,
+ unsigned int b, int k)
+{
+ /* This function combines 2 bit strings of length k
+ into one string of length 2*k such that the individual
+ bits of a and b are ordered:
+ k bit of a -> 2k bit of bitstr
+ k bit of b -> 2k-1 bit of bitstr
+ k-1 bit of a -> 2k-2 bit of bitstr
+ k-1 bit of b -> 2k-3 bit of bitstr
+ k-2 bit of a -> 2k-4 bit of bitstr
+ ...
+ k-(k-1) bit of b -> 2k-(2k-1) bit of bitstr */
+
+ unsigned int bitstr;
+ int bit, ls;
+
+ bitstr = 0;
+ ls = 2 * k;
+ for (bit = k - 1; bit >= 0; bit --)
+ {
+ ls --;
+ bitstr += Bits (a, bit, 1) << ls;
+ ls --;
+ bitstr += Bits (b, bit, 1) << ls;
+ }
+ return bitstr;
+}
+
+
+/**************************************************************************/
+
+int Bits_in_pwr_two_can_sqr_and_store_in (char type)
+{
+ int bits;
+
+ bits = Bits_available (type);
+
+
+ /* Word which is an even no. of bits can store
+ 2^x * 2^x where x = bits / 2. Word which is
+ an odd no. of bits can store 2^x * 2^x where
+ x = (bits + 1) / 2 */
+
+ if (fmod (bits, 2) != 1)
+ return (bits / 2);
+ else
+ return ((int) ((bits + 1) / 2));
+}
Deleted: grass-addons/raster/r.boxcount/bits.h
===================================================================
--- grass-addons/r.boxcount/bits.h 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/bits.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,40 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-#ifndef BITSH
-#define BITSH
-
-unsigned Bits (unsigned int, int, int);
-int Bits_available (char);
-unsigned int Create_chunked_bit_string (unsigned int,
- unsigned int, int);
-unsigned int Create_cyclic_bit_string (unsigned int,
- unsigned int, int);
-int Bits_in_pwr_two_can_sqr_and_store_in (char);
-
-#endif
Copied: grass-addons/raster/r.boxcount/bits.h (from rev 30226, grass-addons/r.boxcount/bits.h)
===================================================================
--- grass-addons/raster/r.boxcount/bits.h (rev 0)
+++ grass-addons/raster/r.boxcount/bits.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,40 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#ifndef BITSH
+#define BITSH
+
+unsigned Bits (unsigned int, int, int);
+int Bits_available (char);
+unsigned int Create_chunked_bit_string (unsigned int,
+ unsigned int, int);
+unsigned int Create_cyclic_bit_string (unsigned int,
+ unsigned int, int);
+int Bits_in_pwr_two_can_sqr_and_store_in (char);
+
+#endif
Deleted: grass-addons/raster/r.boxcount/description.html
===================================================================
--- grass-addons/r.boxcount/description.html 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,200 +0,0 @@
-<H2>r.boxcount</H2>
-
-
-NAME
- r.boxcount - Calculate fractal dimension by box counting.
- (GRASS Raster Program)
-
-SYNOPSIS
- r.boxcount
- r.boxcount help
- r.boxcount [-gnot] input=name [output=name] k=value
- [resolution=value] [saturation=value]
-
-DESCRIPTION
- r.boxcount calculates the box-counting (fractal) dimension
- of a binary raster map. The program counts the minimum
- number of boxes, NbE, of size 1,...,1/2^k that are required
- to cover all non-zero cells. It then uses this information
- to calculate the fractal dimension for each pair of box
- sizes. The results may be saved in a text file, in which
- case they may also be displayed in graph form by invoking
- the program Gnuplot (which, athough not part of GRASS, is
- widely available on UNIX-like systems).
-
- r.boxcount also calculates the box-counting dimension by
- linear regression over a range of box sizes determined by
- the parameters resolution and saturation. The largest box
- size that will be used is 1/resolution (i.e. this is the
- coarsest resolution). The smallest box size is set to that
- for which the results satisfy the condition NbE < saturation
- * the number of non-zero cells in the region (i.e. before
- saturation occurs).
-
-OPTIONS
- The program can be run either non-interactively or
- interactively. To run r.boxcount non-interactively, the
- user must specify at least an input file name and the value
- of k on the command line; any remaining parameters whose
- values are unspecified on the command line will be set to
- their default values (see below).
-
- To run the r.boxcount interactively the user should simply
- type r.boxcount on the command line, in which case the
- program will prompt for parameter values using the standard
- GRASS interface described in the manual entry for parser.
-
-
- Flags:
-
- -g Create a file containing the commands for
- Gnuplot to draw a graph of the results.
- This file will be created in the current
- working directory with the name
- output.gnu. This option requires that an
- output filename has been provided.
-
-
-
-GRASS 4.2.1 Baylor University 1
-
-
-
-
-
-
-r.boxcount <contrib> GRASS Reference Manual <contrib> r.boxcount
-
-
-
- -n Invoke Gnuplot to draw a graph of the
- results now (i.e. automatically). This
- option requires the -g flag.
-
- -o Silently overwrite the files output.dat
- and output.gnu if they already exist in
- the current working directory.
-
- -t Provide only a terse description of the
- results (i.e. provide no information about
- progress and only display the fractal
- dimension calculated by linear
- regression). This option applies to the
- on-screen output only; it does not affect
- that written to output.dat.
-
- Parameters:
-
- input=name The name of a raster map layer containing
- binary values 0 or 1.
-
- output=name Basename for the text file in which the
- program output will be stored. The file
- name.dat is created in the user's current
- working directory according to whether the
- flag -o has been set.
-
- k=value The program calculates the box counting
- dimension of the input map using box sizes
- from 1,...,1/2^k. Note that the current
- region will be covered by 1 box of size 1,
- 4 of size 1/2, 16 of size 1/4, etc. The
- program automatically skips box sizes that
- are smaller than the size of one map cell.
-
- resolution=name Set the largest box size to 1/resolution
- instead of 1. Values of resolution must
- be from the sequence 1, 2, 4, 8, 16, 32,
- 64, 128, 256,...
- Default: 4
-
- saturation=value Saturation is deemed to have occured when
- NbE divided by the number of non-zero
- cells in the map is >= saturation.
- Default: 0.2
-
-NON-INTERACTIVE PROGRAM USE
- The flag -n should not be used if the program is to be
- invoked from a shell script since Gnuplot will wait on the
- user pressing the RETURN key before exiting.
-
-
-
-
-
-GRASS 4.2.1 Baylor University 2
-
-
-
-
-
-
-r.boxcount <contrib> GRASS Reference Manual <contrib> r.boxcount
-
-
-
-NOTES
- r.boxcount may give erroneous results if the current region
- is not square, i.e. the number of rows and columns are not
- equal.
-
- r.boxcount works in the current geographic region with the
- current mask.
-
- The last row of the table of results does not contain a
- value for the fractal dimension, D, because this can only be
- calculated for pairs of box sizes. In the output file this
- missing value is replaced by the dummy value 99.999. In
- both cases the recorded values of D always relate to the box
- sizes recorded in the given row and the following row.
-
- It is a good idea to examine the table and/or graph of
- results before accepting the fractal dimension calculated by
- linear regression based on the default values of resolution
- and saturation. The user may find it helpful to study the
- results of applying r.boxcount to a theoretical fractal such
- as the Sierpinski Gasket. r.rifs may be used to produce
- raster maps containing this and many other well known
- fractals.
-
-SEE ALSO
- r.rifs, parser
-
- For an introduction to the box-counting method of
- calculating fractal dimension see: Peitgen, Jurgens and
- Saupe, 1992, Chaos and Fractals: New Frontiers of Science,
- Springer-Verlag: New York. Chapt. 4.
-
-ACKNOWLEDGEMENTS
- This program was written during the author's tenure of a
- Leverhulme Special Research Fellowhip.
-
-AUTHOR
- Mark Lake, Institute of Archaeology, University College
- London.
-
-NOTICE
- This program is part of the contrib section of the GRASS
- distribution. As such, it is externally contributed code
- that has not been examined or tested by the Office of GRASS
- Integration.
-
-
-
-
-
-
-
-
-
-
-GRASS 4.2.1 Baylor University 3
-
-
-
-
-<H2>AUTHOR</H2>
-
-GRASS Development Team
-
-<p><i>Last changed: $Date: 2005/04/19 16:12:34 $</i>
Copied: grass-addons/raster/r.boxcount/description.html (from rev 30226, grass-addons/r.boxcount/description.html)
===================================================================
--- grass-addons/raster/r.boxcount/description.html (rev 0)
+++ grass-addons/raster/r.boxcount/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,200 @@
+<H2>r.boxcount</H2>
+
+
+NAME
+ r.boxcount - Calculate fractal dimension by box counting.
+ (GRASS Raster Program)
+
+SYNOPSIS
+ r.boxcount
+ r.boxcount help
+ r.boxcount [-gnot] input=name [output=name] k=value
+ [resolution=value] [saturation=value]
+
+DESCRIPTION
+ r.boxcount calculates the box-counting (fractal) dimension
+ of a binary raster map. The program counts the minimum
+ number of boxes, NbE, of size 1,...,1/2^k that are required
+ to cover all non-zero cells. It then uses this information
+ to calculate the fractal dimension for each pair of box
+ sizes. The results may be saved in a text file, in which
+ case they may also be displayed in graph form by invoking
+ the program Gnuplot (which, athough not part of GRASS, is
+ widely available on UNIX-like systems).
+
+ r.boxcount also calculates the box-counting dimension by
+ linear regression over a range of box sizes determined by
+ the parameters resolution and saturation. The largest box
+ size that will be used is 1/resolution (i.e. this is the
+ coarsest resolution). The smallest box size is set to that
+ for which the results satisfy the condition NbE < saturation
+ * the number of non-zero cells in the region (i.e. before
+ saturation occurs).
+
+OPTIONS
+ The program can be run either non-interactively or
+ interactively. To run r.boxcount non-interactively, the
+ user must specify at least an input file name and the value
+ of k on the command line; any remaining parameters whose
+ values are unspecified on the command line will be set to
+ their default values (see below).
+
+ To run the r.boxcount interactively the user should simply
+ type r.boxcount on the command line, in which case the
+ program will prompt for parameter values using the standard
+ GRASS interface described in the manual entry for parser.
+
+
+ Flags:
+
+ -g Create a file containing the commands for
+ Gnuplot to draw a graph of the results.
+ This file will be created in the current
+ working directory with the name
+ output.gnu. This option requires that an
+ output filename has been provided.
+
+
+
+GRASS 4.2.1 Baylor University 1
+
+
+
+
+
+
+r.boxcount <contrib> GRASS Reference Manual <contrib> r.boxcount
+
+
+
+ -n Invoke Gnuplot to draw a graph of the
+ results now (i.e. automatically). This
+ option requires the -g flag.
+
+ -o Silently overwrite the files output.dat
+ and output.gnu if they already exist in
+ the current working directory.
+
+ -t Provide only a terse description of the
+ results (i.e. provide no information about
+ progress and only display the fractal
+ dimension calculated by linear
+ regression). This option applies to the
+ on-screen output only; it does not affect
+ that written to output.dat.
+
+ Parameters:
+
+ input=name The name of a raster map layer containing
+ binary values 0 or 1.
+
+ output=name Basename for the text file in which the
+ program output will be stored. The file
+ name.dat is created in the user's current
+ working directory according to whether the
+ flag -o has been set.
+
+ k=value The program calculates the box counting
+ dimension of the input map using box sizes
+ from 1,...,1/2^k. Note that the current
+ region will be covered by 1 box of size 1,
+ 4 of size 1/2, 16 of size 1/4, etc. The
+ program automatically skips box sizes that
+ are smaller than the size of one map cell.
+
+ resolution=name Set the largest box size to 1/resolution
+ instead of 1. Values of resolution must
+ be from the sequence 1, 2, 4, 8, 16, 32,
+ 64, 128, 256,...
+ Default: 4
+
+ saturation=value Saturation is deemed to have occured when
+ NbE divided by the number of non-zero
+ cells in the map is >= saturation.
+ Default: 0.2
+
+NON-INTERACTIVE PROGRAM USE
+ The flag -n should not be used if the program is to be
+ invoked from a shell script since Gnuplot will wait on the
+ user pressing the RETURN key before exiting.
+
+
+
+
+
+GRASS 4.2.1 Baylor University 2
+
+
+
+
+
+
+r.boxcount <contrib> GRASS Reference Manual <contrib> r.boxcount
+
+
+
+NOTES
+ r.boxcount may give erroneous results if the current region
+ is not square, i.e. the number of rows and columns are not
+ equal.
+
+ r.boxcount works in the current geographic region with the
+ current mask.
+
+ The last row of the table of results does not contain a
+ value for the fractal dimension, D, because this can only be
+ calculated for pairs of box sizes. In the output file this
+ missing value is replaced by the dummy value 99.999. In
+ both cases the recorded values of D always relate to the box
+ sizes recorded in the given row and the following row.
+
+ It is a good idea to examine the table and/or graph of
+ results before accepting the fractal dimension calculated by
+ linear regression based on the default values of resolution
+ and saturation. The user may find it helpful to study the
+ results of applying r.boxcount to a theoretical fractal such
+ as the Sierpinski Gasket. r.rifs may be used to produce
+ raster maps containing this and many other well known
+ fractals.
+
+SEE ALSO
+ r.rifs, parser
+
+ For an introduction to the box-counting method of
+ calculating fractal dimension see: Peitgen, Jurgens and
+ Saupe, 1992, Chaos and Fractals: New Frontiers of Science,
+ Springer-Verlag: New York. Chapt. 4.
+
+ACKNOWLEDGEMENTS
+ This program was written during the author's tenure of a
+ Leverhulme Special Research Fellowhip.
+
+AUTHOR
+ Mark Lake, Institute of Archaeology, University College
+ London.
+
+NOTICE
+ This program is part of the contrib section of the GRASS
+ distribution. As such, it is externally contributed code
+ that has not been examined or tested by the Office of GRASS
+ Integration.
+
+
+
+
+
+
+
+
+
+
+GRASS 4.2.1 Baylor University 3
+
+
+
+
+<H2>AUTHOR</H2>
+
+GRASS Development Team
+
+<p><i>Last changed: $Date: 2005/04/19 16:12:34 $</i>
Deleted: grass-addons/raster/r.boxcount/file.c
===================================================================
--- grass-addons/r.boxcount/file.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/file.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,271 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-
-#include <stdlib.h>
-#include <string.h>
-#include "file.h"
-#include "math.h"
-
-
-void Strip_path (char*, char*);
-
-/**************************************************************************/
-
-FILE* Create_file (char *name, char *suffix, char *message, int overwrite)
-{
- FILE* stream;
-
- strcat (name, suffix);
- stream = fopen (name, "r");
- if ((stream != NULL) && (!overwrite))
- {
- sprintf (message, "File %s exits ", name);
- fclose (stream);
- return NULL;
- }
- else
- {
- stream = fopen (name, "w");
- if (stream == NULL)
- sprintf (message, "Can't create file %s ", name);
- }
-
- return (stream);
-}
-
-
-/**************************************************************************/
-
-void Functional_print (FILE *stream, tRecord *record,
- int k, float *coef)
-{
- int m;
-
- fprintf (stream, "\n 1/Box_size Occupied Log_1/Box_Size Log_Occupied D");
- for (m = 0; m <= k; m ++)
- {
- if (record [m].size != 0)
- {
- fprintf (stream, "\n%12lu ", (unsigned long int) (pow (2, k - m)));
- fprintf (stream, "%12lu ", record [m].occupied);
- fprintf (stream, " %6.3f ", record [m].log_reciprocal_size);
- fprintf (stream, " %6.3f", record [m].log_occupied);
-
- /* Fractal dimension only calculated for pairs
- of points... */
-
- if (m < k)
- fprintf (stream, " %6.3f", record [m].d);
- else
-
- /* ...so skip the last one and print a dummy value */
-
- fprintf (stream, " %6.3f", 99.999);
- }
- }
- fflush (stream);
-}
-
-
-/**************************************************************************/
-
-void Pretty_print (FILE *stream, tRecord *record,
- int k, int r, int points_used,
- float *coef)
-{
- int m;
-
- fprintf (stream, "\n\nResults:");
- fprintf (stream, "\n 1/Box size Occupied Log 1/Box Size Log Occupied D");
- for (m = 0; m <= k; m ++)
- {
- fprintf (stream, "\n%12lu ", (unsigned long int) (pow (2, k - m)));
- if (record [m].size == 0)
- fprintf (stream, " -------------- too few map cells ----------------");
- else
- {
- fprintf (stream, "%12lu ", record [m].occupied);
- fprintf (stream, " %6.3f ", record [m].log_reciprocal_size);
- fprintf (stream, " %6.3f", record [m].log_occupied);
-
- /* Fractal dimension only calculated for pairs
- of points... */
-
- if (m < k)
- fprintf (stream, " %6.3f", record [m].d);
- else
-
- /* ...so skip the last one */
-
- fprintf (stream, " ------");
- }
- }
- if (points_used >= 2)
- {
- fprintf (stream, "\n\nD by regression on 1/box sizes %lu to %lu (inc.) is: %6.3f\n",
- (unsigned long int) pow (2, r),
- (unsigned long int) pow (2, r + points_used - 1),
- coef [1]);
- }
- fflush (stream);
-}
-
-
-/**************************************************************************/
-
-void Print_gnuplot_commands (FILE *stream,
- char *out_dat_name,
- tRecord *record,
- int k,
- int r, int points_used,
- float *coef)
-{
- int m;
- float x, y, max;
- char message [196];
- char filename [128];
-
- /* Strip path from out_dat_name (for tcltkgrass) */
-
- Strip_path (out_dat_name, filename);
-
- sprintf (message, "%s (D by regression on log10 1/s = %5.3f to %5.3f (inc.) is: %6.3f)",
- filename,
- record [k - r].log_reciprocal_size,
- record [k - r - points_used + 1].log_reciprocal_size,
- coef [1]);
-
-
- /* Find greatest x val and add a bit to ensure that plot area
- will be large enough to accomodate last label */
-
- max = 0.0;
- for (m = 0; m <= (k - 1); m ++)
- {
- if (record [m].log_reciprocal_size > max)
- max = record [m].log_reciprocal_size;
- }
-
- /* Ensure that x axis ends on an integer multiple of 1 or 0.5 */
-
- x = (float) ceil ((double) max + 0.2);
- if ((x - max) > 0.5)
- x -= 0.5;
-
- /* Setup the x and y axes */
-
- fprintf (stream, "set xlabel 'log10 1/s'");
- fprintf (stream, "\nset ylabel 'log10 NbE'");
- fprintf (stream, "\nset xrange [0:%f]", x);
- fprintf (stream, "\nset yrange [0:*]");
-
-
- /* Plot fractal dimensions of pairs of points as labels
- halfway between them. Skip biggest box since fractal
- dimension is only calculated for pairs of points */
-
- for (m = 0; m <= (k - 1); m ++)
- {
- if (record [m].size != 0)
- {
- x = record [m + 1].log_reciprocal_size
- + ((record [m].log_reciprocal_size
- - record [m + 1].log_reciprocal_size) / 2);
- y = record [m + 1].log_occupied
- + ((record [m].log_occupied
- - record [m + 1].log_occupied) / 2);
-
- /* Subtract a tiny amount from y to ensure
- that the label clears the plotted line */
-
- fprintf (stream, "\nset label '%6.3f' at %f,%f left",
- record [m].d, x, y - 0.1);
- }
- }
-
-
- /* Issue instructions to plot graph */
-
- fprintf (stream, "\nplot '%s' using ($3):($4) title '%s' with linespoints",
- out_dat_name, message);
- fprintf (stream, "\npause -1 'Hit RETURN to delete graph'");
-}
-
-
-/**************************************************************************/
-
-void Terse_print (FILE *stream,
- int r, int points_used,
- float *coef)
-{
- if (points_used >= 2)
- {
- fprintf (stream, "D by regression on 1/box sizes %lu to %lu (inc.) is: %6.3f\n",
- (unsigned long int) pow (2, r),
- (unsigned long int) pow (2, r + points_used - 1),
- coef [1]);
- }
- else
- {
- /* we need 12 tokens in output for r.boxcount.sh ... */
- fprintf (stream, "D cannot be calculated - using NO DATA value . . -1\n");
- }
- fflush (stream);
-}
-
-/**************************************************************************/
-
-void Strip_path (char *fqfn, char *fn)
-{
- int pos, start, end, i;
-
- end = strlen (fqfn);
- pos = end;
- while ((fqfn [pos] != '/') && (pos >= 0))
- pos --;
- if (pos == -1)
- {
- strcpy (fn, fqfn);
- }
- else
- {
- start = pos + 1;
- i = 0;
- for (pos = start; pos <= end; pos ++)
- {
- fn [i] = fqfn [pos];
- i ++;
- }
- fn [i] = '\0';
- }
-}
-
-
-
-
Copied: grass-addons/raster/r.boxcount/file.c (from rev 30226, grass-addons/r.boxcount/file.c)
===================================================================
--- grass-addons/raster/r.boxcount/file.c (rev 0)
+++ grass-addons/raster/r.boxcount/file.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,271 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include "file.h"
+#include "math.h"
+
+
+void Strip_path (char*, char*);
+
+/**************************************************************************/
+
+FILE* Create_file (char *name, char *suffix, char *message, int overwrite)
+{
+ FILE* stream;
+
+ strcat (name, suffix);
+ stream = fopen (name, "r");
+ if ((stream != NULL) && (!overwrite))
+ {
+ sprintf (message, "File %s exits ", name);
+ fclose (stream);
+ return NULL;
+ }
+ else
+ {
+ stream = fopen (name, "w");
+ if (stream == NULL)
+ sprintf (message, "Can't create file %s ", name);
+ }
+
+ return (stream);
+}
+
+
+/**************************************************************************/
+
+void Functional_print (FILE *stream, tRecord *record,
+ int k, float *coef)
+{
+ int m;
+
+ fprintf (stream, "\n 1/Box_size Occupied Log_1/Box_Size Log_Occupied D");
+ for (m = 0; m <= k; m ++)
+ {
+ if (record [m].size != 0)
+ {
+ fprintf (stream, "\n%12lu ", (unsigned long int) (pow (2, k - m)));
+ fprintf (stream, "%12lu ", record [m].occupied);
+ fprintf (stream, " %6.3f ", record [m].log_reciprocal_size);
+ fprintf (stream, " %6.3f", record [m].log_occupied);
+
+ /* Fractal dimension only calculated for pairs
+ of points... */
+
+ if (m < k)
+ fprintf (stream, " %6.3f", record [m].d);
+ else
+
+ /* ...so skip the last one and print a dummy value */
+
+ fprintf (stream, " %6.3f", 99.999);
+ }
+ }
+ fflush (stream);
+}
+
+
+/**************************************************************************/
+
+void Pretty_print (FILE *stream, tRecord *record,
+ int k, int r, int points_used,
+ float *coef)
+{
+ int m;
+
+ fprintf (stream, "\n\nResults:");
+ fprintf (stream, "\n 1/Box size Occupied Log 1/Box Size Log Occupied D");
+ for (m = 0; m <= k; m ++)
+ {
+ fprintf (stream, "\n%12lu ", (unsigned long int) (pow (2, k - m)));
+ if (record [m].size == 0)
+ fprintf (stream, " -------------- too few map cells ----------------");
+ else
+ {
+ fprintf (stream, "%12lu ", record [m].occupied);
+ fprintf (stream, " %6.3f ", record [m].log_reciprocal_size);
+ fprintf (stream, " %6.3f", record [m].log_occupied);
+
+ /* Fractal dimension only calculated for pairs
+ of points... */
+
+ if (m < k)
+ fprintf (stream, " %6.3f", record [m].d);
+ else
+
+ /* ...so skip the last one */
+
+ fprintf (stream, " ------");
+ }
+ }
+ if (points_used >= 2)
+ {
+ fprintf (stream, "\n\nD by regression on 1/box sizes %lu to %lu (inc.) is: %6.3f\n",
+ (unsigned long int) pow (2, r),
+ (unsigned long int) pow (2, r + points_used - 1),
+ coef [1]);
+ }
+ fflush (stream);
+}
+
+
+/**************************************************************************/
+
+void Print_gnuplot_commands (FILE *stream,
+ char *out_dat_name,
+ tRecord *record,
+ int k,
+ int r, int points_used,
+ float *coef)
+{
+ int m;
+ float x, y, max;
+ char message [196];
+ char filename [128];
+
+ /* Strip path from out_dat_name (for tcltkgrass) */
+
+ Strip_path (out_dat_name, filename);
+
+ sprintf (message, "%s (D by regression on log10 1/s = %5.3f to %5.3f (inc.) is: %6.3f)",
+ filename,
+ record [k - r].log_reciprocal_size,
+ record [k - r - points_used + 1].log_reciprocal_size,
+ coef [1]);
+
+
+ /* Find greatest x val and add a bit to ensure that plot area
+ will be large enough to accomodate last label */
+
+ max = 0.0;
+ for (m = 0; m <= (k - 1); m ++)
+ {
+ if (record [m].log_reciprocal_size > max)
+ max = record [m].log_reciprocal_size;
+ }
+
+ /* Ensure that x axis ends on an integer multiple of 1 or 0.5 */
+
+ x = (float) ceil ((double) max + 0.2);
+ if ((x - max) > 0.5)
+ x -= 0.5;
+
+ /* Setup the x and y axes */
+
+ fprintf (stream, "set xlabel 'log10 1/s'");
+ fprintf (stream, "\nset ylabel 'log10 NbE'");
+ fprintf (stream, "\nset xrange [0:%f]", x);
+ fprintf (stream, "\nset yrange [0:*]");
+
+
+ /* Plot fractal dimensions of pairs of points as labels
+ halfway between them. Skip biggest box since fractal
+ dimension is only calculated for pairs of points */
+
+ for (m = 0; m <= (k - 1); m ++)
+ {
+ if (record [m].size != 0)
+ {
+ x = record [m + 1].log_reciprocal_size
+ + ((record [m].log_reciprocal_size
+ - record [m + 1].log_reciprocal_size) / 2);
+ y = record [m + 1].log_occupied
+ + ((record [m].log_occupied
+ - record [m + 1].log_occupied) / 2);
+
+ /* Subtract a tiny amount from y to ensure
+ that the label clears the plotted line */
+
+ fprintf (stream, "\nset label '%6.3f' at %f,%f left",
+ record [m].d, x, y - 0.1);
+ }
+ }
+
+
+ /* Issue instructions to plot graph */
+
+ fprintf (stream, "\nplot '%s' using ($3):($4) title '%s' with linespoints",
+ out_dat_name, message);
+ fprintf (stream, "\npause -1 'Hit RETURN to delete graph'");
+}
+
+
+/**************************************************************************/
+
+void Terse_print (FILE *stream,
+ int r, int points_used,
+ float *coef)
+{
+ if (points_used >= 2)
+ {
+ fprintf (stream, "D by regression on 1/box sizes %lu to %lu (inc.) is: %6.3f\n",
+ (unsigned long int) pow (2, r),
+ (unsigned long int) pow (2, r + points_used - 1),
+ coef [1]);
+ }
+ else
+ {
+ /* we need 12 tokens in output for r.boxcount.sh ... */
+ fprintf (stream, "D cannot be calculated - using NO DATA value . . -1\n");
+ }
+ fflush (stream);
+}
+
+/**************************************************************************/
+
+void Strip_path (char *fqfn, char *fn)
+{
+ int pos, start, end, i;
+
+ end = strlen (fqfn);
+ pos = end;
+ while ((fqfn [pos] != '/') && (pos >= 0))
+ pos --;
+ if (pos == -1)
+ {
+ strcpy (fn, fqfn);
+ }
+ else
+ {
+ start = pos + 1;
+ i = 0;
+ for (pos = start; pos <= end; pos ++)
+ {
+ fn [i] = fqfn [pos];
+ i ++;
+ }
+ fn [i] = '\0';
+ }
+}
+
+
+
+
Deleted: grass-addons/raster/r.boxcount/file.h
===================================================================
--- grass-addons/r.boxcount/file.h 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/file.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,41 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-#ifndef FILEH
-#define FILEH
-
-#include <stdio.h>
-#include "record.h"
-
-FILE* Create_file (char*, char*, char*, int);
-void Pretty_print (FILE*, tRecord*, int, int, int, float*);
-void Functional_print (FILE*, tRecord*, int, float*);
-void Print_gnuplot_commands (FILE*, char*, tRecord*, int, int, int, float*);
-void Terse_print (FILE*, int, int, float*);
-
-#endif
Copied: grass-addons/raster/r.boxcount/file.h (from rev 30226, grass-addons/r.boxcount/file.h)
===================================================================
--- grass-addons/raster/r.boxcount/file.h (rev 0)
+++ grass-addons/raster/r.boxcount/file.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,41 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#ifndef FILEH
+#define FILEH
+
+#include <stdio.h>
+#include "record.h"
+
+FILE* Create_file (char*, char*, char*, int);
+void Pretty_print (FILE*, tRecord*, int, int, int, float*);
+void Functional_print (FILE*, tRecord*, int, float*);
+void Print_gnuplot_commands (FILE*, char*, tRecord*, int, int, int, float*);
+void Terse_print (FILE*, int, int, float*);
+
+#endif
Deleted: grass-addons/raster/r.boxcount/main.c
===================================================================
--- grass-addons/r.boxcount/main.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/main.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,772 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- * PURPOSE:
- *
- * 1) Calculates the capacity, or box-counting, fractal dimension
- * of a binary raster map. The box-counting dimension is
- * calculated for each pair of box sizes and by linear
- * regression over a range of box sizes.
- *
- * OPTIONS:
- *
- * 1) Functional options determine which box sizes are used
- * in the linear regression.
- *
- * 2) Output options allow detailed results to be reported to
- * stdout, saved to a file and graphed using gnuplot.
- *
- * METHOD:
- *
- * 1) For an introduction to the box-counting method of
- * calculating fractal dimension see:
- * Peitgen, Jurgens and Saupe, 1992, 'Chaos and
- * Fractals: New Frontiers of Science', Springer-Verlag:
- * New York. Chapt. 4.
- *
- * 2) This program is based on the algorithm described in:
- * Liebovitch and Toth, 1989, 'A Fast Algorithm to
- * Determine Fractal Dimension by Box Counting'. In
- * Physics Letters A, vol. 141, pp.386--390.
- *
- * Note, however, that it uses an improved method
- * the same or similar to that suggested by Daniel Kaplan
- * (see note on p.389).
- *
- * 3) For more detailed information see the file 'Method.ps'.
- *
- * FILES:
- *
- * 1) main.c - most of box counting algorithm
- * and calculates the fractal dimension.
- *
- * 2) sort.c - a radix exchange sort used during box counting.
- *
- * 3) bits.c - functions that manipulate bits of a word.
- *
- * 4) file.c - functions that output the results.
- *
- * PORTABILITY:
- *
- * 1) The functions in 'bits.c' manipuluate individual bits
- * and might need modification on a Big-Endian machine
- * (e.g SPARC, 68xxx). This code was developed on an 80586.
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-
-#define MAIN
-
-#define kCells_in_block 10000 /* No of cells per block of memory */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include <limits.h>
-#include <unistd.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "sort.h"
-#include "file.h"
-#include "record.h"
-
-/*
- * global function declaration
- */
-extern CELL f_c(CELL);
-extern FCELL f_f(FCELL);
-extern DCELL f_d(DCELL);
-
-/*
- * function definitions
- */
-
-CELL c_calc(CELL x)
-{
- /* we do nothing exciting here */
- return x;
-}
-
-FCELL f_calc(FCELL x)
-{
- /* we do nothing exciting here */
- return x;
-}
-
-DCELL d_calc(DCELL x)
-{
- /* we do nothing exciting here */
- return x;
-}
-
-/*
- * main function
- * it copies raster input raster file, calling the appropriate function for each
- * data type (CELL, DCELL, FCELL)
- */
-int main(int argc, char *argv[])
-{
-
-/* from r.example */
-
- struct Cell_head cellhd; /* it stores region information,
- and header information of rasters */
- char *name; /* input raster name */
- char *result; /* output raster name */
- char *mapset; /* mapset name */
- void *inrast; /* input buffer */
- unsigned char *outrast; /* output buffer */
- int infd, outfd; /* file descriptor */
- int verbose;
- RASTER_MAP_TYPE data_type; /* type of the map (CELL/DCELL/...) */
- struct History history; /* holds meta-data (title, comments,..) */
-
- struct GModule *module; /* GRASS module for parsing arguments */
-
- struct Flag *flag1; /* flags */
-
-/* from r.boxcount42 */
-
- FILE *out_dat_str = NULL;
- FILE *out_gnuplot_str = NULL;
- char out_dat_name [128], out_gnuplot_name [128];
- int in_map_fd;
- char *current_mapset;
- int projection;
- int array_south, array_north;
- int array_east, array_west;
- int nrows, ncols, left;
- CELL *row_of_cells;
- tRecord *record;
- unsigned long int no_boxes, occupied_boxes;
- unsigned long int cell, no_cells;
- int cells_in_this_block, no_blocks;
- long int block;
- int r, k, max_k, step, m, i, j, n, max;
- int points_used;
- int col, row;
- unsigned int x, y, M, *keys;
- float row_factor, col_factor;
- float **f, a[2][3], coef[2], t;
- float s;
- int tmp1;
- float tmp2;
- char message [192];
- struct Option *input, *output, *k_opt, *r_opt, *s_opt;
- struct Flag *gnuplot, *now, *overwrite, *terse;
-
-
- /***********************************************************************/
- /* Preliminary things */
- /***********************************************************************/
-
- /* Initialize the GIS calls */
-
- G_gisinit(argv[0]);
- G_sleep_on_error (0);
-
- /* initialize module */
- /*
- module = G_define_module();
- module->keywords = _("keyword1, keyword2, keyword3");
- module->description = _("My first raster module");
-*/
- /* Define the different options as defined in gis.h */
- input = G_define_standard_option(G_OPT_R_INPUT);
-
- output = G_define_option ();
- output->key = "output";
- output->type = TYPE_STRING;
- output->required = NO;
- output->description = "Output text file for data points";
-
- k_opt = G_define_option ();
- k_opt->key = "k";
- k_opt->type = TYPE_INTEGER;
- k_opt->required = YES;
- k_opt->description = "Max 1/box size is 2^k ";
- s_opt = G_define_option ();
- s_opt->key = "saturation";
- s_opt->type = TYPE_DOUBLE;
- s_opt->description =
- "Occupied boxes as fraction of data points at saturation";
- s_opt->answer = "0.2";
- r_opt = G_define_option ();
- r_opt->key = "resolution";
- r_opt->type = TYPE_INTEGER;
- r_opt->description =
- "Smallest 1/box size to use in regression (= 1, 2, 4,...)";
- r_opt->answer = "4";
-
- /* Define the different flags */
-
- gnuplot = G_define_flag ();
- gnuplot->key = 'g';
- gnuplot->description = "Create file with instructions for Gnuplot ";
- now = G_define_flag ();
- now->key = 'n';
- now->description = "Gnuplot graph now (i.e. on exit from r.boxcount) ";
- overwrite = G_define_flag ();
- overwrite->key = 'o';
- overwrite->description = "Overwrite output file if it exists ";
- terse = G_define_flag ();
- terse->key = 't';
- terse->description =
- "Terse output only (i.e. only report D calculated by regression) ";
-
- /* Parse the command line arguments, or get parameters if none */
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
-
- /* Check that database is UTM or x-y referenced */
-/* fkindl: why is that?? - I kicked it out and it works with M28, too... */
-
-/*
- projection = G_projection ();
- if (projection != PROJECTION_UTM && projection != PROJECTION_XY)
- G_fatal_error ("Can't use database\n\t\t- not UTM or X-Y");
- */
-
- /* Open the raster map for reading */
-
- /* returns NULL if the map was not found in any mapset,
- * mapset name otherwise */
-
- name=input->answer;
- mapset = G_find_cell2(name, "");
- if (mapset == NULL)
- G_fatal_error(_("cell file [%s] not found"), name);
-
-
- /* determine the inputmap type (CELL/FCELL/DCELL) */
- data_type = G_raster_map_type(name, mapset);
-
- /* G_open_cell_old - returns file destriptor (>0) */
- if ((in_map_fd = G_open_cell_old(name, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), name);
-
- /* controlling, if we can open input raster */
- if (G_get_cellhd(name, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), name);
- G_debug(3, "number of rows %d", cellhd.rows);
-
-
-
-
- /* If required, open the output files for writing */
-
- if (output->answer != NULL)
- {
- strcpy (out_dat_name, output->answer);
- out_dat_str = Create_file (out_dat_name, ".dat",
- message, overwrite->answer);
- if (out_dat_str == NULL)
- G_fatal_error (message);
-
-
- if (gnuplot->answer)
- {
- strcpy (out_gnuplot_name, output->answer);
- out_gnuplot_str = Create_file (out_gnuplot_name, ".gnu",
- message, overwrite->answer);
- if (out_gnuplot_str == NULL)
- G_fatal_error (message);
- }
- }
-
-
-
-
- /* Calculate array limits */
-
- array_south = G_window_rows () - 1;
- array_north = 0;
- array_west = 0;
- array_east = G_window_cols () - 1;
- nrows = G_window_rows ();
- ncols = G_window_cols ();
-
-
- /* Allocate memory for one row of map */
-
- row_of_cells = G_allocate_cell_buf ();
-
-
- /* Copy the numeric parameters from the parser */
-
- sscanf (k_opt->answer, "%d", &k);
- sscanf (s_opt->answer, "%f", &s);
- sscanf (r_opt->answer, "%d", &r);
-
-
-
- /* Convert r to the value of m which produces it.
- I.e. r (as given) = 2^r (as calculated here).
- This value of r determines the lowest value of
- m used when choosing which box sizes to include
- in the linear regression (e.g. m = r, r+1, r+2,..., k) */
-
- tmp2 = frexp ((double) r, &r);
- r -= 1; /* Because frexp returns 0.5 <= exp < 1 */
- if (r < 0)
- G_fatal_error ("Problem with r ");
- else
- if (r > k)
- G_fatal_error ("r must be <= 2^k "); /* That is, r as given */
-
-
- /* Check that k is not too big. There are two
- limits:
- 1) The no. of bits available in each element
- of the array `keys'
- 2) It must be possible to store the max. no. of
- boxes (2^k * 2^k) in `cell' which, since it is an
- array subscript, must necessarily be an integer
- When an unsigned int and an unsigned long int are
- the same size then 2) is the ultimate limit, but
- we check both 1) and 2) just in case */
-
- max_k = Bits_in_pwr_two_can_sqr_and_store_in ('l') ; /* Check 2) */
- if ((Bits_available ('i') / 2) < max_k)
- max_k = Bits_available ('i') / 2; /* Check 1) */
-
- if (k > max_k)
- {
- sprintf (message, "Sorry, k is too big - %d is max. for this machine ",
- max_k);
- G_fatal_error (message);
- }
-
-
- /* Allocate memory for structure used to record results */
-
- record = (tRecord*) malloc ((k + 1) * sizeof (tRecord));
- if (record == NULL)
- G_fatal_error ("Unable to allocate memory for record structure ");
-
-
- /* Calculate reduction factors for normalisation.
- N.B., the result of applying these factors is truncated to the
- integer component. This results in normalisation
- between 0, ..., 2^k -1. The reason for subtracting
- a very small value (0.01) here is to ensure that the
- coordinates are spread across all boxes. If the factors
- were calculated using 2^k - 1 then the box with coordinates
- 2^k -1, 2^k -1 would only contain the point with coordinates
- ncols, nrows. */
-
- row_factor = (nrows - 1) / (pow (2, k) - 0.01);
- col_factor = (ncols - 1) / (pow (2, k) - 0.01);
-
-
- /***********************************************************************/
- /* Normalise coordinates of every map cell which has a non-zero
- value */
- /***********************************************************************/
-
- /* Allocate memory for array of keys to store the normalised
- coordinates of every map cell which has a non-zero value.
- Do this by allocating memory in blocks adequate for kCells_in_block
- keys. The rationale here is to avoid reading from disk twice:
- once to count non-zero cells and then again to normalise them
- after having allocated memory according to the number counted.
- A linked-list might seem the obvious solution, but it suffers from
- two problems: it ultimately requires more memory, and it
- complicates the sort process. */
-
- block = kCells_in_block * sizeof (unsigned int);
- keys = (unsigned int*) malloc (block);
- if (keys == NULL)
- G_fatal_error ("Not enough memory");
-
-
- /* Trawl through disk file, reallocating memory after
- every kCells_in_block non-zero cells */
-
- if (!terse->answer)
- fprintf (stderr, "\nNormalising coordinates ");
- no_blocks = 1;
- cells_in_this_block = 0;
- cell = 0;
- for (row = 0; row < nrows; row ++)
- {
- G_get_map_row (in_map_fd, row_of_cells, row);
- for (col = 0; col < ncols; col ++)
- {
- if (row_of_cells [col] != 0)
- {
- x = (unsigned int) (col / col_factor);
- y = (unsigned int) (row / row_factor);
-
- /* Combine coordinates into one string
- of 2k bits such that highest bit
- of x is followed by the highest bit of
- y, etc. (see bits.c for a detailed
- explanation) */
-
- keys [cell] = Create_cyclic_bit_string (x, y, k);
- cell ++;
-
- /* Check that allocated memory has not been exhausted */
-
- cells_in_this_block ++;
- if (cells_in_this_block == kCells_in_block)
- {
- no_blocks ++;
- cells_in_this_block = 0;
- keys = (unsigned int*) realloc ((void*)
- keys, no_blocks * block);
- if (keys == NULL)
- {
- sprintf (message, "Not enough memory for cells %d+ ",
- cell);
- G_fatal_error (message);
- }
- }
- }
- }
- }
- no_cells = cell;
- if (!terse->answer)
- fprintf (stderr, " - done "); fflush (stderr);
-
-
- /***********************************************************************/
- /* Sort all keys */
- /***********************************************************************/
-
- if (!terse->answer)
- fprintf (stderr, "\nSorting coordinates ");
- Radix_exchange_sort_array (keys,
- 0, no_cells -1, (2 * k) - 1);
-
-#ifdef DEBUG
- for (cell = 0; cell < no_cells; cell ++)
- fprintf (stdout, "\ncell = %12ld, key = %u ",
- cell, keys [cell]); fflush (stdout);
-#endif
-
- if (!terse->answer)
- fprintf (stderr, " - done "); fflush (stderr);
-
-
- /***********************************************************************/
- /* Count the minimum number of boxes of size 1/2^m,
- where m=k, k-1, k-2,...,0, required to cover every
- map cell which has a non-zero value. Note that the
- minimum resolution, r, and saturation, s, are ignored
- here - they only constrain the box sizes used during
- regression */
- /***********************************************************************/
-
- if (!terse->answer)
- fprintf (stderr, "\nCounting occupied boxes ");
- for (m = 0; m <= k; m ++)
- {
-
-#ifdef DEBUG
- no_boxes = (unsigned long int) (pow (2, k - m) * pow (2, k - m));
- fprintf (stderr, "\n%12lu boxes ", no_boxes);
-#endif
-
- /* If the number of boxes in a row or column is,
- for this m, greater than the number of map cells... */
-
- if ((pow (2, k) / pow (2, m) > nrows) ||
- (pow (2, k) / pow (2, m) > ncols))
- {
- /* Record that count was skipped for this m */
-
- record [m].size = 0;
- }
- else
- {
- /* Record 1/box size for this m */
-
- record [m].size = pow (2, k - m);
-
- /* Calculate mask, M, for each coordinate x and y... */
-
- M = 0;
- for (i = k; i > m; i --)
- M += pow (2, i - 1);
-
- /* ...and combine using same cyclic ordering
- that was used to create the keys */
-
- M = Create_cyclic_bit_string (M, M, k);
-
-#ifdef DEBUG
- fprintf (stdout, "\n\nm = %d, M = %u ", m, M); fflush (stdout);
-#endif
-
- /* Mask all keys in 'keys' */
-
- for (cell = 0; cell < no_cells; cell ++)
- keys [cell] = keys [cell] & M;
-
- /* Count number of changes in the array 'keys'.
- Remember that the cyclic ordering of bits
- in the keys and mask ensures that progressive
- masking, e.g 11111111, 11111100, 11110000,...
- does not unsort the array */
-
- occupied_boxes = 1;
- for (cell = 1; cell < no_cells; cell ++)
- {
- if (keys [cell] != keys [cell - 1])
- occupied_boxes ++;
- }
- record [m].occupied = occupied_boxes;
- }
- }
- if (!terse->answer)
- fprintf (stderr, " - done" );
-
- /* Free memory */
-
- free (keys);
-
-
- /***********************************************************************/
- /* Calculate log values for all m */
- /***********************************************************************/
-
- for (m = 0; m <= k; m ++)
- {
- record [m].log_occupied = log10 (record [m].occupied);
- record [m].log_reciprocal_size = log10 (record [m].size);
- }
-
-
- /***********************************************************************/
- /* Calculate fractal dimension for each pair of data points */
- /***********************************************************************/
-
- if (!terse->answer)
- fprintf (stderr, "\nCalculating dimensions ");
- for (m = 0; m < k; m++)
- {
- record [m].d = (record [m].log_occupied - record [m + 1].log_occupied) /
- (record [m].log_reciprocal_size - record [m + 1].log_reciprocal_size);
- }
-
-
- /***********************************************************************/
- /* Calculate the fractal dimension by least squares regression
- on box sizes 1/2^m for m = r, r+1, r+2,...,m_max. m_max
- is determined by the saturation fraction, s, such that
- m_max = k if s = 1.0, but m_max < k if s < 1.0. */
- /***********************************************************************/
-
- /* Set up arrays for vectors used to calculate matrix for
- Gaussian elimination */
-
- f = (float**) malloc (3 * sizeof (float*));
- if (f == NULL)
- G_fatal_error ("Not enough memory ");
- else
- for (i = 0; i < 3; i++)
- {
- f[i] = (float*) malloc ((k - r + 1) * sizeof (float));
- if (f == NULL)
- G_fatal_error ("Not enough memory ");
- }
-
- /* Compute vectors appropriate for function of the form
- y = a*x + b. Ignore data for m < r (resolution too
- poor) and for occupied > no_cells * s (minimum number of boxes
- required to cover non-zero map cells is too close to number
- of non-zero map cells - i.e. saturation reached) */
-
- points_used = 0;
- m = k - r;;
- while (( m >= 0) && (record [m].occupied < (no_cells * s))
- && record [m].size != 0)
- {
- f[0][m] = 1.0; /* b */
- f[1][m] = record [m].log_reciprocal_size; /* a*x vals */
- f[2][m] = record [m].log_occupied; /* y vals */
-
-#ifdef DEBUG
- fprintf (stderr, "\nm = %d f[0] = %6.3f, f[1] = %6.3f, f[2] = %6.3f", m, f[0][m], f[1][m], f[2][m]); fflush (stderr);
-#endif
-
- m --;
- points_used ++;
- }
-
- if (points_used >= 2)
- {
- /* Compute matrix for Gaussian elimination */
-
- for (i = 0; i <= 1; i++)
- {
- for (j = 0; j <= 2; j++)
- {
- t = 0.0;
- for (n = (k - r); n > m; n --)
- t += f[i][n] * f[j][n];
- a [i][j] = t;
-
-#ifdef DEBUG
- fprintf (stderr, "\na[%d][%d] = %f", i, j, a[i][j]); fflush (stderr);
-#endif
-
- }
- }
-
- /* Forward-elimination phase */
-
- for (i = 0; i <= 1; i++)
- {
- max = i;
- for (j = i + 1; j <= 1; j++)
- {
- if (abs (a [j][i]) > abs (a [max][i]))
- max = j;
- }
- for (n = i; n <= 2; n++)
- {
- t = a [i][n];
- a [i][n] = a [max][n];
- a [max][n] = t;
- }
- for (j = i + 1; j <= 1; j++)
- for (n = 2; n >= i; n--)
- a [j][n] -= a [i][n] * a [j][i] / a [i][i];
- }
-
- /* Backward-substitution phase.
- The coefficients a is returned in coef [1] and
- the coefficient b in coef [0] */
-
- for (j = 1; j >= 0; j --)
- {
- t = 0.0;
- for (n = j + 1; n <= 1; n++)
- t += a [j][n] * coef [n];
- coef [j] = (a [j][2] - t) / a [j][j];
- }
-
- if (!terse->answer)
- fprintf (stderr, " - done ");
- }
-
- else
- {
- if (!terse->answer)
- {
- fprintf (stderr, " - two few data points");
- fprintf (stderr,
- "\n for least squares regression");
- fprintf (stderr,
- "\n (try increasing s) ");
- }
- }
-
-#ifdef DEBUG
- fprintf (stderr, "\n\n eqn is %6.3f x + %6.3f ", coef [1], coef [0]);
- fflush (stderr);
-#endif
-
-
-
- /***********************************************************************/
- /* Report results */
- /***********************************************************************/
-
- /* If requested, output results to files */
-
- if (out_dat_str != NULL)
- {
- if (!terse->answer)
- fprintf (stderr, "\nSaving results to file '%s' ", out_dat_name);
- Functional_print (out_dat_str, record, k, coef);
- if (out_gnuplot_str != NULL)
- {
- if (!terse->answer)
- fprintf (stderr, "\nSaving commands to file '%s' ",
- out_gnuplot_name);
- Print_gnuplot_commands (out_gnuplot_str, out_dat_name,
- record, k, r, points_used, coef);
- }
- }
-
-
- /* Output results to screen in pretty or terse form */
-
- if (!terse->answer)
- Pretty_print (stdout, record, k, r, points_used, coef);
- else
- Terse_print (stdout, r, points_used, coef);
-
-
- /* Tidy up and, if requested, produce gnuplot graph */
-
- G_close_cell (in_map_fd);
-
- if (out_dat_str != NULL)
- {
- fclose (out_dat_str);
- if (out_gnuplot_str != NULL)
- {
- fclose (out_gnuplot_str);
-
- /* Spawn gnuplot in place of this process.
- Note that we don't get here unless the command file
- was created */
-
- if (now->answer)
- {
-
- /* This is the neat way of spawning Gnuplot
- because it allows control from the GRASS shell.
- But we don't use it becuase it isn't compatible
- with tcltkgrass... */
-
- /* fprintf (stderr, "\n");
- execlp ("gnuplot", "gnuplot", out_gnuplot_name, 0);*/
-
- /* ...so we do this instead (which is messy because it
- adds another xterm). */
-
- char command_line [256];
- sprintf (command_line,
- "xterm -e gnuplot %s &", out_gnuplot_name);
- system (command_line);
- }
- }
- }
-
-
- /* We only get here if not doing gnuplot */
-
- if (!terse->answer)
- fprintf (stderr, "\nJob finished\n");
-
- return (0);
-}
Copied: grass-addons/raster/r.boxcount/main.c (from rev 30226, grass-addons/r.boxcount/main.c)
===================================================================
--- grass-addons/raster/r.boxcount/main.c (rev 0)
+++ grass-addons/raster/r.boxcount/main.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,772 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ * PURPOSE:
+ *
+ * 1) Calculates the capacity, or box-counting, fractal dimension
+ * of a binary raster map. The box-counting dimension is
+ * calculated for each pair of box sizes and by linear
+ * regression over a range of box sizes.
+ *
+ * OPTIONS:
+ *
+ * 1) Functional options determine which box sizes are used
+ * in the linear regression.
+ *
+ * 2) Output options allow detailed results to be reported to
+ * stdout, saved to a file and graphed using gnuplot.
+ *
+ * METHOD:
+ *
+ * 1) For an introduction to the box-counting method of
+ * calculating fractal dimension see:
+ * Peitgen, Jurgens and Saupe, 1992, 'Chaos and
+ * Fractals: New Frontiers of Science', Springer-Verlag:
+ * New York. Chapt. 4.
+ *
+ * 2) This program is based on the algorithm described in:
+ * Liebovitch and Toth, 1989, 'A Fast Algorithm to
+ * Determine Fractal Dimension by Box Counting'. In
+ * Physics Letters A, vol. 141, pp.386--390.
+ *
+ * Note, however, that it uses an improved method
+ * the same or similar to that suggested by Daniel Kaplan
+ * (see note on p.389).
+ *
+ * 3) For more detailed information see the file 'Method.ps'.
+ *
+ * FILES:
+ *
+ * 1) main.c - most of box counting algorithm
+ * and calculates the fractal dimension.
+ *
+ * 2) sort.c - a radix exchange sort used during box counting.
+ *
+ * 3) bits.c - functions that manipulate bits of a word.
+ *
+ * 4) file.c - functions that output the results.
+ *
+ * PORTABILITY:
+ *
+ * 1) The functions in 'bits.c' manipuluate individual bits
+ * and might need modification on a Big-Endian machine
+ * (e.g SPARC, 68xxx). This code was developed on an 80586.
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#define MAIN
+
+#define kCells_in_block 10000 /* No of cells per block of memory */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <limits.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "sort.h"
+#include "file.h"
+#include "record.h"
+
+/*
+ * global function declaration
+ */
+extern CELL f_c(CELL);
+extern FCELL f_f(FCELL);
+extern DCELL f_d(DCELL);
+
+/*
+ * function definitions
+ */
+
+CELL c_calc(CELL x)
+{
+ /* we do nothing exciting here */
+ return x;
+}
+
+FCELL f_calc(FCELL x)
+{
+ /* we do nothing exciting here */
+ return x;
+}
+
+DCELL d_calc(DCELL x)
+{
+ /* we do nothing exciting here */
+ return x;
+}
+
+/*
+ * main function
+ * it copies raster input raster file, calling the appropriate function for each
+ * data type (CELL, DCELL, FCELL)
+ */
+int main(int argc, char *argv[])
+{
+
+/* from r.example */
+
+ struct Cell_head cellhd; /* it stores region information,
+ and header information of rasters */
+ char *name; /* input raster name */
+ char *result; /* output raster name */
+ char *mapset; /* mapset name */
+ void *inrast; /* input buffer */
+ unsigned char *outrast; /* output buffer */
+ int infd, outfd; /* file descriptor */
+ int verbose;
+ RASTER_MAP_TYPE data_type; /* type of the map (CELL/DCELL/...) */
+ struct History history; /* holds meta-data (title, comments,..) */
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+
+ struct Flag *flag1; /* flags */
+
+/* from r.boxcount42 */
+
+ FILE *out_dat_str = NULL;
+ FILE *out_gnuplot_str = NULL;
+ char out_dat_name [128], out_gnuplot_name [128];
+ int in_map_fd;
+ char *current_mapset;
+ int projection;
+ int array_south, array_north;
+ int array_east, array_west;
+ int nrows, ncols, left;
+ CELL *row_of_cells;
+ tRecord *record;
+ unsigned long int no_boxes, occupied_boxes;
+ unsigned long int cell, no_cells;
+ int cells_in_this_block, no_blocks;
+ long int block;
+ int r, k, max_k, step, m, i, j, n, max;
+ int points_used;
+ int col, row;
+ unsigned int x, y, M, *keys;
+ float row_factor, col_factor;
+ float **f, a[2][3], coef[2], t;
+ float s;
+ int tmp1;
+ float tmp2;
+ char message [192];
+ struct Option *input, *output, *k_opt, *r_opt, *s_opt;
+ struct Flag *gnuplot, *now, *overwrite, *terse;
+
+
+ /***********************************************************************/
+ /* Preliminary things */
+ /***********************************************************************/
+
+ /* Initialize the GIS calls */
+
+ G_gisinit(argv[0]);
+ G_sleep_on_error (0);
+
+ /* initialize module */
+ /*
+ module = G_define_module();
+ module->keywords = _("keyword1, keyword2, keyword3");
+ module->description = _("My first raster module");
+*/
+ /* Define the different options as defined in gis.h */
+ input = G_define_standard_option(G_OPT_R_INPUT);
+
+ output = G_define_option ();
+ output->key = "output";
+ output->type = TYPE_STRING;
+ output->required = NO;
+ output->description = "Output text file for data points";
+
+ k_opt = G_define_option ();
+ k_opt->key = "k";
+ k_opt->type = TYPE_INTEGER;
+ k_opt->required = YES;
+ k_opt->description = "Max 1/box size is 2^k ";
+ s_opt = G_define_option ();
+ s_opt->key = "saturation";
+ s_opt->type = TYPE_DOUBLE;
+ s_opt->description =
+ "Occupied boxes as fraction of data points at saturation";
+ s_opt->answer = "0.2";
+ r_opt = G_define_option ();
+ r_opt->key = "resolution";
+ r_opt->type = TYPE_INTEGER;
+ r_opt->description =
+ "Smallest 1/box size to use in regression (= 1, 2, 4,...)";
+ r_opt->answer = "4";
+
+ /* Define the different flags */
+
+ gnuplot = G_define_flag ();
+ gnuplot->key = 'g';
+ gnuplot->description = "Create file with instructions for Gnuplot ";
+ now = G_define_flag ();
+ now->key = 'n';
+ now->description = "Gnuplot graph now (i.e. on exit from r.boxcount) ";
+ overwrite = G_define_flag ();
+ overwrite->key = 'o';
+ overwrite->description = "Overwrite output file if it exists ";
+ terse = G_define_flag ();
+ terse->key = 't';
+ terse->description =
+ "Terse output only (i.e. only report D calculated by regression) ";
+
+ /* Parse the command line arguments, or get parameters if none */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* Check that database is UTM or x-y referenced */
+/* fkindl: why is that?? - I kicked it out and it works with M28, too... */
+
+/*
+ projection = G_projection ();
+ if (projection != PROJECTION_UTM && projection != PROJECTION_XY)
+ G_fatal_error ("Can't use database\n\t\t- not UTM or X-Y");
+ */
+
+ /* Open the raster map for reading */
+
+ /* returns NULL if the map was not found in any mapset,
+ * mapset name otherwise */
+
+ name=input->answer;
+ mapset = G_find_cell2(name, "");
+ if (mapset == NULL)
+ G_fatal_error(_("cell file [%s] not found"), name);
+
+
+ /* determine the inputmap type (CELL/FCELL/DCELL) */
+ data_type = G_raster_map_type(name, mapset);
+
+ /* G_open_cell_old - returns file destriptor (>0) */
+ if ((in_map_fd = G_open_cell_old(name, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), name);
+
+ /* controlling, if we can open input raster */
+ if (G_get_cellhd(name, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), name);
+ G_debug(3, "number of rows %d", cellhd.rows);
+
+
+
+
+ /* If required, open the output files for writing */
+
+ if (output->answer != NULL)
+ {
+ strcpy (out_dat_name, output->answer);
+ out_dat_str = Create_file (out_dat_name, ".dat",
+ message, overwrite->answer);
+ if (out_dat_str == NULL)
+ G_fatal_error (message);
+
+
+ if (gnuplot->answer)
+ {
+ strcpy (out_gnuplot_name, output->answer);
+ out_gnuplot_str = Create_file (out_gnuplot_name, ".gnu",
+ message, overwrite->answer);
+ if (out_gnuplot_str == NULL)
+ G_fatal_error (message);
+ }
+ }
+
+
+
+
+ /* Calculate array limits */
+
+ array_south = G_window_rows () - 1;
+ array_north = 0;
+ array_west = 0;
+ array_east = G_window_cols () - 1;
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+
+
+ /* Allocate memory for one row of map */
+
+ row_of_cells = G_allocate_cell_buf ();
+
+
+ /* Copy the numeric parameters from the parser */
+
+ sscanf (k_opt->answer, "%d", &k);
+ sscanf (s_opt->answer, "%f", &s);
+ sscanf (r_opt->answer, "%d", &r);
+
+
+
+ /* Convert r to the value of m which produces it.
+ I.e. r (as given) = 2^r (as calculated here).
+ This value of r determines the lowest value of
+ m used when choosing which box sizes to include
+ in the linear regression (e.g. m = r, r+1, r+2,..., k) */
+
+ tmp2 = frexp ((double) r, &r);
+ r -= 1; /* Because frexp returns 0.5 <= exp < 1 */
+ if (r < 0)
+ G_fatal_error ("Problem with r ");
+ else
+ if (r > k)
+ G_fatal_error ("r must be <= 2^k "); /* That is, r as given */
+
+
+ /* Check that k is not too big. There are two
+ limits:
+ 1) The no. of bits available in each element
+ of the array `keys'
+ 2) It must be possible to store the max. no. of
+ boxes (2^k * 2^k) in `cell' which, since it is an
+ array subscript, must necessarily be an integer
+ When an unsigned int and an unsigned long int are
+ the same size then 2) is the ultimate limit, but
+ we check both 1) and 2) just in case */
+
+ max_k = Bits_in_pwr_two_can_sqr_and_store_in ('l') ; /* Check 2) */
+ if ((Bits_available ('i') / 2) < max_k)
+ max_k = Bits_available ('i') / 2; /* Check 1) */
+
+ if (k > max_k)
+ {
+ sprintf (message, "Sorry, k is too big - %d is max. for this machine ",
+ max_k);
+ G_fatal_error (message);
+ }
+
+
+ /* Allocate memory for structure used to record results */
+
+ record = (tRecord*) malloc ((k + 1) * sizeof (tRecord));
+ if (record == NULL)
+ G_fatal_error ("Unable to allocate memory for record structure ");
+
+
+ /* Calculate reduction factors for normalisation.
+ N.B., the result of applying these factors is truncated to the
+ integer component. This results in normalisation
+ between 0, ..., 2^k -1. The reason for subtracting
+ a very small value (0.01) here is to ensure that the
+ coordinates are spread across all boxes. If the factors
+ were calculated using 2^k - 1 then the box with coordinates
+ 2^k -1, 2^k -1 would only contain the point with coordinates
+ ncols, nrows. */
+
+ row_factor = (nrows - 1) / (pow (2, k) - 0.01);
+ col_factor = (ncols - 1) / (pow (2, k) - 0.01);
+
+
+ /***********************************************************************/
+ /* Normalise coordinates of every map cell which has a non-zero
+ value */
+ /***********************************************************************/
+
+ /* Allocate memory for array of keys to store the normalised
+ coordinates of every map cell which has a non-zero value.
+ Do this by allocating memory in blocks adequate for kCells_in_block
+ keys. The rationale here is to avoid reading from disk twice:
+ once to count non-zero cells and then again to normalise them
+ after having allocated memory according to the number counted.
+ A linked-list might seem the obvious solution, but it suffers from
+ two problems: it ultimately requires more memory, and it
+ complicates the sort process. */
+
+ block = kCells_in_block * sizeof (unsigned int);
+ keys = (unsigned int*) malloc (block);
+ if (keys == NULL)
+ G_fatal_error ("Not enough memory");
+
+
+ /* Trawl through disk file, reallocating memory after
+ every kCells_in_block non-zero cells */
+
+ if (!terse->answer)
+ fprintf (stderr, "\nNormalising coordinates ");
+ no_blocks = 1;
+ cells_in_this_block = 0;
+ cell = 0;
+ for (row = 0; row < nrows; row ++)
+ {
+ G_get_map_row (in_map_fd, row_of_cells, row);
+ for (col = 0; col < ncols; col ++)
+ {
+ if (row_of_cells [col] != 0)
+ {
+ x = (unsigned int) (col / col_factor);
+ y = (unsigned int) (row / row_factor);
+
+ /* Combine coordinates into one string
+ of 2k bits such that highest bit
+ of x is followed by the highest bit of
+ y, etc. (see bits.c for a detailed
+ explanation) */
+
+ keys [cell] = Create_cyclic_bit_string (x, y, k);
+ cell ++;
+
+ /* Check that allocated memory has not been exhausted */
+
+ cells_in_this_block ++;
+ if (cells_in_this_block == kCells_in_block)
+ {
+ no_blocks ++;
+ cells_in_this_block = 0;
+ keys = (unsigned int*) realloc ((void*)
+ keys, no_blocks * block);
+ if (keys == NULL)
+ {
+ sprintf (message, "Not enough memory for cells %d+ ",
+ cell);
+ G_fatal_error (message);
+ }
+ }
+ }
+ }
+ }
+ no_cells = cell;
+ if (!terse->answer)
+ fprintf (stderr, " - done "); fflush (stderr);
+
+
+ /***********************************************************************/
+ /* Sort all keys */
+ /***********************************************************************/
+
+ if (!terse->answer)
+ fprintf (stderr, "\nSorting coordinates ");
+ Radix_exchange_sort_array (keys,
+ 0, no_cells -1, (2 * k) - 1);
+
+#ifdef DEBUG
+ for (cell = 0; cell < no_cells; cell ++)
+ fprintf (stdout, "\ncell = %12ld, key = %u ",
+ cell, keys [cell]); fflush (stdout);
+#endif
+
+ if (!terse->answer)
+ fprintf (stderr, " - done "); fflush (stderr);
+
+
+ /***********************************************************************/
+ /* Count the minimum number of boxes of size 1/2^m,
+ where m=k, k-1, k-2,...,0, required to cover every
+ map cell which has a non-zero value. Note that the
+ minimum resolution, r, and saturation, s, are ignored
+ here - they only constrain the box sizes used during
+ regression */
+ /***********************************************************************/
+
+ if (!terse->answer)
+ fprintf (stderr, "\nCounting occupied boxes ");
+ for (m = 0; m <= k; m ++)
+ {
+
+#ifdef DEBUG
+ no_boxes = (unsigned long int) (pow (2, k - m) * pow (2, k - m));
+ fprintf (stderr, "\n%12lu boxes ", no_boxes);
+#endif
+
+ /* If the number of boxes in a row or column is,
+ for this m, greater than the number of map cells... */
+
+ if ((pow (2, k) / pow (2, m) > nrows) ||
+ (pow (2, k) / pow (2, m) > ncols))
+ {
+ /* Record that count was skipped for this m */
+
+ record [m].size = 0;
+ }
+ else
+ {
+ /* Record 1/box size for this m */
+
+ record [m].size = pow (2, k - m);
+
+ /* Calculate mask, M, for each coordinate x and y... */
+
+ M = 0;
+ for (i = k; i > m; i --)
+ M += pow (2, i - 1);
+
+ /* ...and combine using same cyclic ordering
+ that was used to create the keys */
+
+ M = Create_cyclic_bit_string (M, M, k);
+
+#ifdef DEBUG
+ fprintf (stdout, "\n\nm = %d, M = %u ", m, M); fflush (stdout);
+#endif
+
+ /* Mask all keys in 'keys' */
+
+ for (cell = 0; cell < no_cells; cell ++)
+ keys [cell] = keys [cell] & M;
+
+ /* Count number of changes in the array 'keys'.
+ Remember that the cyclic ordering of bits
+ in the keys and mask ensures that progressive
+ masking, e.g 11111111, 11111100, 11110000,...
+ does not unsort the array */
+
+ occupied_boxes = 1;
+ for (cell = 1; cell < no_cells; cell ++)
+ {
+ if (keys [cell] != keys [cell - 1])
+ occupied_boxes ++;
+ }
+ record [m].occupied = occupied_boxes;
+ }
+ }
+ if (!terse->answer)
+ fprintf (stderr, " - done" );
+
+ /* Free memory */
+
+ free (keys);
+
+
+ /***********************************************************************/
+ /* Calculate log values for all m */
+ /***********************************************************************/
+
+ for (m = 0; m <= k; m ++)
+ {
+ record [m].log_occupied = log10 (record [m].occupied);
+ record [m].log_reciprocal_size = log10 (record [m].size);
+ }
+
+
+ /***********************************************************************/
+ /* Calculate fractal dimension for each pair of data points */
+ /***********************************************************************/
+
+ if (!terse->answer)
+ fprintf (stderr, "\nCalculating dimensions ");
+ for (m = 0; m < k; m++)
+ {
+ record [m].d = (record [m].log_occupied - record [m + 1].log_occupied) /
+ (record [m].log_reciprocal_size - record [m + 1].log_reciprocal_size);
+ }
+
+
+ /***********************************************************************/
+ /* Calculate the fractal dimension by least squares regression
+ on box sizes 1/2^m for m = r, r+1, r+2,...,m_max. m_max
+ is determined by the saturation fraction, s, such that
+ m_max = k if s = 1.0, but m_max < k if s < 1.0. */
+ /***********************************************************************/
+
+ /* Set up arrays for vectors used to calculate matrix for
+ Gaussian elimination */
+
+ f = (float**) malloc (3 * sizeof (float*));
+ if (f == NULL)
+ G_fatal_error ("Not enough memory ");
+ else
+ for (i = 0; i < 3; i++)
+ {
+ f[i] = (float*) malloc ((k - r + 1) * sizeof (float));
+ if (f == NULL)
+ G_fatal_error ("Not enough memory ");
+ }
+
+ /* Compute vectors appropriate for function of the form
+ y = a*x + b. Ignore data for m < r (resolution too
+ poor) and for occupied > no_cells * s (minimum number of boxes
+ required to cover non-zero map cells is too close to number
+ of non-zero map cells - i.e. saturation reached) */
+
+ points_used = 0;
+ m = k - r;;
+ while (( m >= 0) && (record [m].occupied < (no_cells * s))
+ && record [m].size != 0)
+ {
+ f[0][m] = 1.0; /* b */
+ f[1][m] = record [m].log_reciprocal_size; /* a*x vals */
+ f[2][m] = record [m].log_occupied; /* y vals */
+
+#ifdef DEBUG
+ fprintf (stderr, "\nm = %d f[0] = %6.3f, f[1] = %6.3f, f[2] = %6.3f", m, f[0][m], f[1][m], f[2][m]); fflush (stderr);
+#endif
+
+ m --;
+ points_used ++;
+ }
+
+ if (points_used >= 2)
+ {
+ /* Compute matrix for Gaussian elimination */
+
+ for (i = 0; i <= 1; i++)
+ {
+ for (j = 0; j <= 2; j++)
+ {
+ t = 0.0;
+ for (n = (k - r); n > m; n --)
+ t += f[i][n] * f[j][n];
+ a [i][j] = t;
+
+#ifdef DEBUG
+ fprintf (stderr, "\na[%d][%d] = %f", i, j, a[i][j]); fflush (stderr);
+#endif
+
+ }
+ }
+
+ /* Forward-elimination phase */
+
+ for (i = 0; i <= 1; i++)
+ {
+ max = i;
+ for (j = i + 1; j <= 1; j++)
+ {
+ if (abs (a [j][i]) > abs (a [max][i]))
+ max = j;
+ }
+ for (n = i; n <= 2; n++)
+ {
+ t = a [i][n];
+ a [i][n] = a [max][n];
+ a [max][n] = t;
+ }
+ for (j = i + 1; j <= 1; j++)
+ for (n = 2; n >= i; n--)
+ a [j][n] -= a [i][n] * a [j][i] / a [i][i];
+ }
+
+ /* Backward-substitution phase.
+ The coefficients a is returned in coef [1] and
+ the coefficient b in coef [0] */
+
+ for (j = 1; j >= 0; j --)
+ {
+ t = 0.0;
+ for (n = j + 1; n <= 1; n++)
+ t += a [j][n] * coef [n];
+ coef [j] = (a [j][2] - t) / a [j][j];
+ }
+
+ if (!terse->answer)
+ fprintf (stderr, " - done ");
+ }
+
+ else
+ {
+ if (!terse->answer)
+ {
+ fprintf (stderr, " - two few data points");
+ fprintf (stderr,
+ "\n for least squares regression");
+ fprintf (stderr,
+ "\n (try increasing s) ");
+ }
+ }
+
+#ifdef DEBUG
+ fprintf (stderr, "\n\n eqn is %6.3f x + %6.3f ", coef [1], coef [0]);
+ fflush (stderr);
+#endif
+
+
+
+ /***********************************************************************/
+ /* Report results */
+ /***********************************************************************/
+
+ /* If requested, output results to files */
+
+ if (out_dat_str != NULL)
+ {
+ if (!terse->answer)
+ fprintf (stderr, "\nSaving results to file '%s' ", out_dat_name);
+ Functional_print (out_dat_str, record, k, coef);
+ if (out_gnuplot_str != NULL)
+ {
+ if (!terse->answer)
+ fprintf (stderr, "\nSaving commands to file '%s' ",
+ out_gnuplot_name);
+ Print_gnuplot_commands (out_gnuplot_str, out_dat_name,
+ record, k, r, points_used, coef);
+ }
+ }
+
+
+ /* Output results to screen in pretty or terse form */
+
+ if (!terse->answer)
+ Pretty_print (stdout, record, k, r, points_used, coef);
+ else
+ Terse_print (stdout, r, points_used, coef);
+
+
+ /* Tidy up and, if requested, produce gnuplot graph */
+
+ G_close_cell (in_map_fd);
+
+ if (out_dat_str != NULL)
+ {
+ fclose (out_dat_str);
+ if (out_gnuplot_str != NULL)
+ {
+ fclose (out_gnuplot_str);
+
+ /* Spawn gnuplot in place of this process.
+ Note that we don't get here unless the command file
+ was created */
+
+ if (now->answer)
+ {
+
+ /* This is the neat way of spawning Gnuplot
+ because it allows control from the GRASS shell.
+ But we don't use it becuase it isn't compatible
+ with tcltkgrass... */
+
+ /* fprintf (stderr, "\n");
+ execlp ("gnuplot", "gnuplot", out_gnuplot_name, 0);*/
+
+ /* ...so we do this instead (which is messy because it
+ adds another xterm). */
+
+ char command_line [256];
+ sprintf (command_line,
+ "xterm -e gnuplot %s &", out_gnuplot_name);
+ system (command_line);
+ }
+ }
+ }
+
+
+ /* We only get here if not doing gnuplot */
+
+ if (!terse->answer)
+ fprintf (stderr, "\nJob finished\n");
+
+ return (0);
+}
Deleted: grass-addons/raster/r.boxcount/record.h
===================================================================
--- grass-addons/r.boxcount/record.h 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/record.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,41 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-#ifndef RECORDH
-#define RECORDH
-
-typedef struct
-{
- unsigned long int occupied;
- float log_occupied;
- float size;
- float log_reciprocal_size;
- float d;
-} tRecord;
-
-#endif
Copied: grass-addons/raster/r.boxcount/record.h (from rev 30226, grass-addons/r.boxcount/record.h)
===================================================================
--- grass-addons/raster/r.boxcount/record.h (rev 0)
+++ grass-addons/raster/r.boxcount/record.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,41 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#ifndef RECORDH
+#define RECORDH
+
+typedef struct
+{
+ unsigned long int occupied;
+ float log_occupied;
+ float size;
+ float log_reciprocal_size;
+ float d;
+} tRecord;
+
+#endif
Deleted: grass-addons/raster/r.boxcount/sort.c
===================================================================
--- grass-addons/r.boxcount/sort.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/sort.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,75 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-
-#include <stdlib.h>
-#include <stdio.h>
-#include "bits.h"
-
-
-/**************************************************************************/
-
-void Radix_exchange_sort_array (unsigned int *keys,
- long int begin, long int end,
- int bits)
-{
- /* See R. Sedgewick, 1990, Algorithms in C, Addison-Wesley:
- Reading, MA, chapt. 10 for a discussion of this sort
- algorithm */
-
- long int upper, lower;
- unsigned int old_key;
-
-
- if (end > begin && bits >= 0)
- {
- lower = begin;
- upper = end;
-
- while (lower != upper)
- {
- while ((Bits (keys [lower], bits, 1) == 0)
- && (lower < upper))
- lower ++;
-
- while ((Bits (keys [upper], bits, 1) != 0)
- && (upper > lower))
- upper --;
-
- old_key = keys [lower];
- keys [lower] = keys [upper];
- keys [upper] = old_key;
- }
-
- if (Bits (keys [end], bits, 1) == 0)
- upper ++;
-
- Radix_exchange_sort_array (keys, begin, upper - 1, bits - 1);
- Radix_exchange_sort_array (keys, upper, end, bits - 1);
- }
-}
Copied: grass-addons/raster/r.boxcount/sort.c (from rev 30226, grass-addons/r.boxcount/sort.c)
===================================================================
--- grass-addons/raster/r.boxcount/sort.c (rev 0)
+++ grass-addons/raster/r.boxcount/sort.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,75 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include "bits.h"
+
+
+/**************************************************************************/
+
+void Radix_exchange_sort_array (unsigned int *keys,
+ long int begin, long int end,
+ int bits)
+{
+ /* See R. Sedgewick, 1990, Algorithms in C, Addison-Wesley:
+ Reading, MA, chapt. 10 for a discussion of this sort
+ algorithm */
+
+ long int upper, lower;
+ unsigned int old_key;
+
+
+ if (end > begin && bits >= 0)
+ {
+ lower = begin;
+ upper = end;
+
+ while (lower != upper)
+ {
+ while ((Bits (keys [lower], bits, 1) == 0)
+ && (lower < upper))
+ lower ++;
+
+ while ((Bits (keys [upper], bits, 1) != 0)
+ && (upper > lower))
+ upper --;
+
+ old_key = keys [lower];
+ keys [lower] = keys [upper];
+ keys [upper] = old_key;
+ }
+
+ if (Bits (keys [end], bits, 1) == 0)
+ upper ++;
+
+ Radix_exchange_sort_array (keys, begin, upper - 1, bits - 1);
+ Radix_exchange_sort_array (keys, upper, end, bits - 1);
+ }
+}
Deleted: grass-addons/raster/r.boxcount/sort.h
===================================================================
--- grass-addons/r.boxcount/sort.h 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount/sort.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,34 +0,0 @@
-/****************************************************************************
- *
- * MODULE: r.boxcount
- * AUTHOR(S):
- *
- * Original author:
- * Mark Lake 14/5/99
- *
- * University College London
- * Institute of Archaeology
- * 31-34 Gordon Square
- * London. WC1H 0PY
- * email: mark.lake at ucl.ac.uk
- *
- * Adaptations for grass63:
- * Florian Kindl, 2006-10-02
- * University of Innsbruck
- * Institute of Geography
- * email: florian.kindl at uibk.ac.at
- *
- *
- * COPYRIGHT: (C) 2008 by the authors
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-#ifndef SORTH
-#define SORTH
-
-void Radix_exchange_sort_array (unsigned int*, long int, long int, int);
-
-#endif
Copied: grass-addons/raster/r.boxcount/sort.h (from rev 30226, grass-addons/r.boxcount/sort.h)
===================================================================
--- grass-addons/raster/r.boxcount/sort.h (rev 0)
+++ grass-addons/raster/r.boxcount/sort.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,34 @@
+/****************************************************************************
+ *
+ * MODULE: r.boxcount
+ * AUTHOR(S):
+ *
+ * Original author:
+ * Mark Lake 14/5/99
+ *
+ * University College London
+ * Institute of Archaeology
+ * 31-34 Gordon Square
+ * London. WC1H 0PY
+ * email: mark.lake at ucl.ac.uk
+ *
+ * Adaptations for grass63:
+ * Florian Kindl, 2006-10-02
+ * University of Innsbruck
+ * Institute of Geography
+ * email: florian.kindl at uibk.ac.at
+ *
+ *
+ * COPYRIGHT: (C) 2008 by the authors
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#ifndef SORTH
+#define SORTH
+
+void Radix_exchange_sort_array (unsigned int*, long int, long int, int);
+
+#endif
Copied: grass-addons/raster/r.boxcount.sh (from rev 30225, grass-addons/r.boxcount.sh)
Deleted: grass-addons/raster/r.boxcount.sh/Makefile
===================================================================
--- grass-addons/r.boxcount.sh/Makefile 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount.sh/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,7 +0,0 @@
-MODULE_TOPDIR = ../..
-
-PGM=r.boxcount.sh
-
-include $(MODULE_TOPDIR)/include/Make/Script.make
-
-default: script
Copied: grass-addons/raster/r.boxcount.sh/Makefile (from rev 30226, grass-addons/r.boxcount.sh/Makefile)
===================================================================
--- grass-addons/raster/r.boxcount.sh/Makefile (rev 0)
+++ grass-addons/raster/r.boxcount.sh/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.boxcount.sh
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Deleted: grass-addons/raster/r.boxcount.sh/description.html
===================================================================
--- grass-addons/r.boxcount.sh/description.html 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount.sh/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,56 +0,0 @@
-<H2>DESCRIPTION</H2>
-
-
-<EM>r.boxcount</EM>Allows the user to study how the boxcounting fractal dimension varies across a raster map.
-<br>
-<pre>
-boxcount.sh is a shell script which allows the user to study how the
-boxcounting fractal dimension varies across a raster map. To to do
-this it repeatedly calls r.boxcount in a moving window. The moving
-window is a square region of a size specified by the user. The centre
-of this window is repeatedly moved by an increment also specified by
-the user. The region used for analysis is the current region. The
-region covered by output data is smaller by whatever distance is
-required to prevent an edge effect given the specified window
-size. Note that since the script changes the GRASS region other
-concurrent GRASS processes could be adversely affected. boxcount.sh
-does, however, restore the starting region once analysis is complete.
-
-boxcount.sh produces an ascii text file listing the fractal dimension
-for each location over which the moving window was centred. This
-output can be converted into a raster map as follows:
-1) Process the output using awk to remove blank lines and
- extraneous text, e.g.
- awk '/D/ { print $12}' outputfile > newfile
-2) Copy the header information and paste it at the
- start of newfile
-3) Import the ascii raster data using r.in.ascii, e.g.
- r.in.ascii input=newfile output=rastermap mult=1000
-
-INSTALLATION
-
-boxcount.sh is useless unless you have already installed r.boxcount.
-
-1) Copy the source file (this file) into the scripts directory:
-
-cp boxcount.sh $GISBASE/scripts
-
-2) Make the source file executable:
-
-chmod ugo+x $GISBASE/scripts/boxcount.sh
-
-3) Have fun!
- boxcount.sh mapname
-
-ACKNOWLEDGEMENTS
-
-This program was written during the author's tenure of a
-Leverhulme Special Research Fellowhip.
-</pre>
-
-
-<H2>AUTHOR</H2>
-
-Mark Lake, 1/9/99
-adaptations for GRASS GIS 6.3: Florian Kindl, 2006-10-04
-<p><i>Last changed: $Date: 2006/10/04 16:50:17 $</i>
Copied: grass-addons/raster/r.boxcount.sh/description.html (from rev 30226, grass-addons/r.boxcount.sh/description.html)
===================================================================
--- grass-addons/raster/r.boxcount.sh/description.html (rev 0)
+++ grass-addons/raster/r.boxcount.sh/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,56 @@
+<H2>DESCRIPTION</H2>
+
+
+<EM>r.boxcount</EM>Allows the user to study how the boxcounting fractal dimension varies across a raster map.
+<br>
+<pre>
+boxcount.sh is a shell script which allows the user to study how the
+boxcounting fractal dimension varies across a raster map. To to do
+this it repeatedly calls r.boxcount in a moving window. The moving
+window is a square region of a size specified by the user. The centre
+of this window is repeatedly moved by an increment also specified by
+the user. The region used for analysis is the current region. The
+region covered by output data is smaller by whatever distance is
+required to prevent an edge effect given the specified window
+size. Note that since the script changes the GRASS region other
+concurrent GRASS processes could be adversely affected. boxcount.sh
+does, however, restore the starting region once analysis is complete.
+
+boxcount.sh produces an ascii text file listing the fractal dimension
+for each location over which the moving window was centred. This
+output can be converted into a raster map as follows:
+1) Process the output using awk to remove blank lines and
+ extraneous text, e.g.
+ awk '/D/ { print $12}' outputfile > newfile
+2) Copy the header information and paste it at the
+ start of newfile
+3) Import the ascii raster data using r.in.ascii, e.g.
+ r.in.ascii input=newfile output=rastermap mult=1000
+
+INSTALLATION
+
+boxcount.sh is useless unless you have already installed r.boxcount.
+
+1) Copy the source file (this file) into the scripts directory:
+
+cp boxcount.sh $GISBASE/scripts
+
+2) Make the source file executable:
+
+chmod ugo+x $GISBASE/scripts/boxcount.sh
+
+3) Have fun!
+ boxcount.sh mapname
+
+ACKNOWLEDGEMENTS
+
+This program was written during the author's tenure of a
+Leverhulme Special Research Fellowhip.
+</pre>
+
+
+<H2>AUTHOR</H2>
+
+Mark Lake, 1/9/99
+adaptations for GRASS GIS 6.3: Florian Kindl, 2006-10-04
+<p><i>Last changed: $Date: 2006/10/04 16:50:17 $</i>
Deleted: grass-addons/raster/r.boxcount.sh/r.boxcount.sh
===================================================================
--- grass-addons/r.boxcount.sh/r.boxcount.sh 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/raster/r.boxcount.sh/r.boxcount.sh 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,261 +0,0 @@
-#!/bin/bash
-
-###############################################################################
-#
-# MODULE: boxcount.sh
-#
-# AUTHORS:
-#
-# Original author:
-# Mark Lake 1/9/99
-# University College London
-# Institute of Archaeology
-# 31-34 Gordon Square
-# London. WC1H 0PY
-# email: mark.lake at ucl.ac.uk
-#
-# Adaptations for grass63:
-# Florian Kindl, 2006-10-02
-# University of Innsbruck
-# Institute of Geography
-# email: florian.kindl at uibk.ac.at
-#
-# PURPOSE: Study how the boxcounting fractal dimension varies across a raster map.
-# COPYRIGHT: (C) 2008 by the authors.
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-###############################################################################
-
-
-#%Module
-#% description: Study how the boxcounting fractal dimension varies across a raster map.
-#% keywords: raster, fractality
-#%End
-#
-#%option
-#% key: input
-#% type: string
-#% gisprompt: in,cell,raster
-#% description: map to be analysed
-#% required : yes
-#%end
-#
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: name for fractal dimension map
-#% required : yes
-#%end
-#
-#% option
-#% key: boxsize
-#% type: integer
-#% description: size of moving window (16,32,64,128,256,512... are best)
-#% required : yes
-#% end
-#
-#% option
-#% key: gridsize
-#% type: integer
-#% description: distance between centre of each window
-#% required : yes
-#% end
-#
-#% option
-#% key: k
-#% type: integer
-#% description: Max 1/box size is 2^k
-#% required : yes
-#% end
-#
-#% option
-#% key: saturation
-#% type: double
-#% description: Occupied boxes as fraction of data points at saturation
-#% required : yes
-#% end
-#
-#% option
-#% key: resolution
-#% type: integer
-#% description: Smallest 1/box size to use in regression (= 1, 2, 4,...)
-#% required : yes
-#% end
-
-
-# Check whether running GRASS
-
-if [ -z $GISRC ]
- then
- echo "Sorry, you are not running GRASS " 1>&2
- exit 1
-fi
-
-# parsing arguments
-
-if [ "$1" != "@ARGS_PARSED@" ]
-then
- exec g.parser "$0" "$@"
-fi
-
-
-# check if we have awk
- if [ ! -x "`which awk`" ] ; then
- echo "ERROR: awk required, please install awk or gawk first" 1>&2
- exit 1
- fi
-
-# Read GRASS environment variables
-eval `g.gisenv`
-
-
-# Read current region into n, s, e, w
-eval `g.region -g`
-
-# is equivalent to this:
-#n=`g.region -p | awk ' /north:/ { print $2 }'`
-#s=`g.region -p | awk ' /south:/ { print $2 }'`
-#e=`g.region -p | awk ' /east:/ { print $2 }'`
-#w=`g.region -p | awk ' /west:/ { print $2 }'`
-
-
-# Save map to be analysed
-binarymap=$GIS_OPT_INPUT
-
-# Get output file name
-record=$GIS_OPT_OUTPUT
-
-# Get parameters for r.boxcount
-k=$GIS_OPT_K
-saturation=$GIS_OPT_SATURATION
-resolution=$GIS_OPT_RESOLUTION
-
-#echo "$GIS_OPT_SATURATION"
-#echo "$k $saturation $resolution"
-
-
-# create tempfiles
-tempfile="`g.tempfile $$`_${record}_fract"
-ascirast="`g.tempfile $$`_${record}_fract_rast"
-#tempfile="/tmp/fract"
-#ascirast="/tmp/fract_rast"
-
-# Get boxsize (i.e. region size for r.boxcount)
-boxsize=$GIS_OPT_BOXSIZE
-
-twiceboxsize=`expr $boxsize + $boxsize`
-halfboxsize=`expr $boxsize / 2`
-
-# Get gridsize (i.e. distance between centre of each region
-# used for boxcounting). Can be greater or less than boxsize.
-gridsize=$GIS_OPT_GRIDSIZE
-
-halfgridsize=`expr $gridsize / 2`
-
-# Calculate appropriate region to eliminate edge effect
-
-if [ $boxsize -ge $gridsize ]
-then
- startnorth=`expr $n - $halfboxsize`
- stopnorth=`expr $s + $halfboxsize + $gridsize`
- startwest=`expr $w + $halfboxsize`
- stopwest=`expr $e - $halfboxsize - $gridsize`
-else
- startnorth=`expr $n - $halfgridsize`
- stopnorth=`expr $s + $halfgridsize + $gridsize`
- startwest=`expr $w + $halfgridsize`
- stopwest=`expr $e - $halfgridsize - $gridsize`
-fi
-
-# Loop over grid
-
-newgridnorth=$startnorth
-row=0
-countnorth=0
-while [ $newgridnorth -ge $stopnorth ]
-do
- row=`expr $row + 1`
- newgridnorth=`expr $startnorth - $countnorth`
- countnorth=`expr $countnorth + $gridsize`
-
- col=0
- countwest=0
- newgridwest=$startwest
- while [ $newgridwest -le $stopwest ]
- do
-
- col=`expr $col + 1`
- echo -ne "\r" 1>&2
- echo -n "row = $row col = $col " 1>&2
- newgridwest=`expr $startwest + $countwest`
- countwest=`expr $countwest + $gridsize`
-
-# Calculate box for this point on the grid
-
- newnorth=`expr $newgridnorth + $halfboxsize`
- newsouth=`expr $newgridnorth - $halfboxsize`
- neweast=`expr $newgridwest + $halfboxsize`
- newwest=`expr $newgridwest - $halfboxsize`
-
-# eval `g.region n=$newnorth s=$newsouth e=$neweast w=$newwest`
- g.region n=$newnorth s=$newsouth e=$neweast w=$newwest
-
-# eval `r.boxcount input=$1 k=9 res=4 sat=1 -t >> /tmp/fract_tmp`
- #echo " "
- echo "row $row col $col" >> $tempfile
- r.boxcount -t input=$binarymap k=$k res=$resolution sat=$saturation >> $tempfile
-
- done
-done
-
-# Provide info. about appropriate header for use with r.in.ascii
-
-#Header for r.in.ascii"
-mapnorth=`expr $startnorth + $halfgridsize`
-echo "north: $mapnorth" > ${ascirast}
-mapsouth=`expr $newgridnorth - $halfgridsize`
-echo "south: $mapsouth" >>${ascirast}
-mapeast=`expr $newgridwest + $halfgridsize`
-echo "east: $mapeast" >>${ascirast}
-mapwest=`expr $startwest - $halfgridsize`
-echo "west: $mapwest" >>${ascirast}
-echo "rows: $row" >>${ascirast}
-echo "cols: $col" >>${ascirast}
-
-
-# cat the result into header file
-awk '/D/ { print $12}' $tempfile >> ${ascirast}
-
-# check if we got data:
-lines=0
-lines=`cat $tempfile|wc -l`
-echo $lines 1>&2
-
-if [ $lines -le 0 ]
-then
- echo "WARNING: Window too small - could not calculate D!" 1>&2
- echo "Please choose a bigger moving window. Stopped." 1>&2
- rm -f $tempfile ${ascirast}
- exit
-else
- # success:
- echo "Importing results..." 1>&2
- #import it into GRASS421:
- #r.in.ascii in=${ascirast} out=$record mult=100
-
- #import it into GRASS 6:
- r.in.ascii in=${ascirast} out=$record nv=-1
-
- #remove the tmp.files
- #rm -f $tempfile ${ascirast}
-
- echo "The resulting map showing the fractal dimension is called $record" 1>&2
- echo "Finished." 1>&2
-fi
-
-# Reset region
-g.region n=$n s=$s e=$e w=$w
-
-exit 0
Copied: grass-addons/raster/r.boxcount.sh/r.boxcount.sh (from rev 30226, grass-addons/r.boxcount.sh/r.boxcount.sh)
===================================================================
--- grass-addons/raster/r.boxcount.sh/r.boxcount.sh (rev 0)
+++ grass-addons/raster/r.boxcount.sh/r.boxcount.sh 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,261 @@
+#!/bin/bash
+
+###############################################################################
+#
+# MODULE: boxcount.sh
+#
+# AUTHORS:
+#
+# Original author:
+# Mark Lake 1/9/99
+# University College London
+# Institute of Archaeology
+# 31-34 Gordon Square
+# London. WC1H 0PY
+# email: mark.lake at ucl.ac.uk
+#
+# Adaptations for grass63:
+# Florian Kindl, 2006-10-02
+# University of Innsbruck
+# Institute of Geography
+# email: florian.kindl at uibk.ac.at
+#
+# PURPOSE: Study how the boxcounting fractal dimension varies across a raster map.
+# COPYRIGHT: (C) 2008 by the authors.
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+###############################################################################
+
+
+#%Module
+#% description: Study how the boxcounting fractal dimension varies across a raster map.
+#% keywords: raster, fractality
+#%End
+#
+#%option
+#% key: input
+#% type: string
+#% gisprompt: in,cell,raster
+#% description: map to be analysed
+#% required : yes
+#%end
+#
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: name for fractal dimension map
+#% required : yes
+#%end
+#
+#% option
+#% key: boxsize
+#% type: integer
+#% description: size of moving window (16,32,64,128,256,512... are best)
+#% required : yes
+#% end
+#
+#% option
+#% key: gridsize
+#% type: integer
+#% description: distance between centre of each window
+#% required : yes
+#% end
+#
+#% option
+#% key: k
+#% type: integer
+#% description: Max 1/box size is 2^k
+#% required : yes
+#% end
+#
+#% option
+#% key: saturation
+#% type: double
+#% description: Occupied boxes as fraction of data points at saturation
+#% required : yes
+#% end
+#
+#% option
+#% key: resolution
+#% type: integer
+#% description: Smallest 1/box size to use in regression (= 1, 2, 4,...)
+#% required : yes
+#% end
+
+
+# Check whether running GRASS
+
+if [ -z $GISRC ]
+ then
+ echo "Sorry, you are not running GRASS " 1>&2
+ exit 1
+fi
+
+# parsing arguments
+
+if [ "$1" != "@ARGS_PARSED@" ]
+then
+ exec g.parser "$0" "$@"
+fi
+
+
+# check if we have awk
+ if [ ! -x "`which awk`" ] ; then
+ echo "ERROR: awk required, please install awk or gawk first" 1>&2
+ exit 1
+ fi
+
+# Read GRASS environment variables
+eval `g.gisenv`
+
+
+# Read current region into n, s, e, w
+eval `g.region -g`
+
+# is equivalent to this:
+#n=`g.region -p | awk ' /north:/ { print $2 }'`
+#s=`g.region -p | awk ' /south:/ { print $2 }'`
+#e=`g.region -p | awk ' /east:/ { print $2 }'`
+#w=`g.region -p | awk ' /west:/ { print $2 }'`
+
+
+# Save map to be analysed
+binarymap=$GIS_OPT_INPUT
+
+# Get output file name
+record=$GIS_OPT_OUTPUT
+
+# Get parameters for r.boxcount
+k=$GIS_OPT_K
+saturation=$GIS_OPT_SATURATION
+resolution=$GIS_OPT_RESOLUTION
+
+#echo "$GIS_OPT_SATURATION"
+#echo "$k $saturation $resolution"
+
+
+# create tempfiles
+tempfile="`g.tempfile $$`_${record}_fract"
+ascirast="`g.tempfile $$`_${record}_fract_rast"
+#tempfile="/tmp/fract"
+#ascirast="/tmp/fract_rast"
+
+# Get boxsize (i.e. region size for r.boxcount)
+boxsize=$GIS_OPT_BOXSIZE
+
+twiceboxsize=`expr $boxsize + $boxsize`
+halfboxsize=`expr $boxsize / 2`
+
+# Get gridsize (i.e. distance between centre of each region
+# used for boxcounting). Can be greater or less than boxsize.
+gridsize=$GIS_OPT_GRIDSIZE
+
+halfgridsize=`expr $gridsize / 2`
+
+# Calculate appropriate region to eliminate edge effect
+
+if [ $boxsize -ge $gridsize ]
+then
+ startnorth=`expr $n - $halfboxsize`
+ stopnorth=`expr $s + $halfboxsize + $gridsize`
+ startwest=`expr $w + $halfboxsize`
+ stopwest=`expr $e - $halfboxsize - $gridsize`
+else
+ startnorth=`expr $n - $halfgridsize`
+ stopnorth=`expr $s + $halfgridsize + $gridsize`
+ startwest=`expr $w + $halfgridsize`
+ stopwest=`expr $e - $halfgridsize - $gridsize`
+fi
+
+# Loop over grid
+
+newgridnorth=$startnorth
+row=0
+countnorth=0
+while [ $newgridnorth -ge $stopnorth ]
+do
+ row=`expr $row + 1`
+ newgridnorth=`expr $startnorth - $countnorth`
+ countnorth=`expr $countnorth + $gridsize`
+
+ col=0
+ countwest=0
+ newgridwest=$startwest
+ while [ $newgridwest -le $stopwest ]
+ do
+
+ col=`expr $col + 1`
+ echo -ne "\r" 1>&2
+ echo -n "row = $row col = $col " 1>&2
+ newgridwest=`expr $startwest + $countwest`
+ countwest=`expr $countwest + $gridsize`
+
+# Calculate box for this point on the grid
+
+ newnorth=`expr $newgridnorth + $halfboxsize`
+ newsouth=`expr $newgridnorth - $halfboxsize`
+ neweast=`expr $newgridwest + $halfboxsize`
+ newwest=`expr $newgridwest - $halfboxsize`
+
+# eval `g.region n=$newnorth s=$newsouth e=$neweast w=$newwest`
+ g.region n=$newnorth s=$newsouth e=$neweast w=$newwest
+
+# eval `r.boxcount input=$1 k=9 res=4 sat=1 -t >> /tmp/fract_tmp`
+ #echo " "
+ echo "row $row col $col" >> $tempfile
+ r.boxcount -t input=$binarymap k=$k res=$resolution sat=$saturation >> $tempfile
+
+ done
+done
+
+# Provide info. about appropriate header for use with r.in.ascii
+
+#Header for r.in.ascii"
+mapnorth=`expr $startnorth + $halfgridsize`
+echo "north: $mapnorth" > ${ascirast}
+mapsouth=`expr $newgridnorth - $halfgridsize`
+echo "south: $mapsouth" >>${ascirast}
+mapeast=`expr $newgridwest + $halfgridsize`
+echo "east: $mapeast" >>${ascirast}
+mapwest=`expr $startwest - $halfgridsize`
+echo "west: $mapwest" >>${ascirast}
+echo "rows: $row" >>${ascirast}
+echo "cols: $col" >>${ascirast}
+
+
+# cat the result into header file
+awk '/D/ { print $12}' $tempfile >> ${ascirast}
+
+# check if we got data:
+lines=0
+lines=`cat $tempfile|wc -l`
+echo $lines 1>&2
+
+if [ $lines -le 0 ]
+then
+ echo "WARNING: Window too small - could not calculate D!" 1>&2
+ echo "Please choose a bigger moving window. Stopped." 1>&2
+ rm -f $tempfile ${ascirast}
+ exit
+else
+ # success:
+ echo "Importing results..." 1>&2
+ #import it into GRASS421:
+ #r.in.ascii in=${ascirast} out=$record mult=100
+
+ #import it into GRASS 6:
+ r.in.ascii in=${ascirast} out=$record nv=-1
+
+ #remove the tmp.files
+ #rm -f $tempfile ${ascirast}
+
+ echo "The resulting map showing the fractal dimension is called $record" 1>&2
+ echo "Finished." 1>&2
+fi
+
+# Reset region
+g.region n=$n s=$s e=$e w=$w
+
+exit 0
Copied: grass-addons/vector/v.strahler (from rev 30225, grass-addons/v.strahler)
Deleted: grass-addons/vector/v.strahler/Makefile
===================================================================
--- grass-addons/v.strahler/Makefile 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,17 +0,0 @@
-
-MODULE_TOPDIR = ../..
-
-PGM = v.strahler
-
-LIBES = $(DISPLAYLIB) $(RASTERLIB) $(VECTLIB) $(GRAPHLIB) $(DBMILIB) $(GISLIB)
-DEPENDENCIES = $(DISPLAYDEP) $(RASTERDEP) $(VECTDEP) $(GISDEP)
-EXTRA_INC = $(VECT_INC)
-EXTRA_CFLAGS = $(VECT_CFLAGS)
-
-include $(MODULE_TOPDIR)/include/Make/Module.make
-
-default: cmd
-
-
-
-
Copied: grass-addons/vector/v.strahler/Makefile (from rev 30226, grass-addons/v.strahler/Makefile)
===================================================================
--- grass-addons/vector/v.strahler/Makefile (rev 0)
+++ grass-addons/vector/v.strahler/Makefile 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,17 @@
+
+MODULE_TOPDIR = ../..
+
+PGM = v.strahler
+
+LIBES = $(DISPLAYLIB) $(RASTERLIB) $(VECTLIB) $(GRAPHLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(DISPLAYDEP) $(RASTERDEP) $(VECTDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
+
+
+
Deleted: grass-addons/vector/v.strahler/description.html
===================================================================
--- grass-addons/v.strahler/description.html 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,87 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>v.strahler</em> calculates the Strahler Order for all lines of a given dendritic network. The input vector map must be free of cycles. For the elaboration a new imported network or a network extracted from DEM by using <em><a HREF="r.watershed.html">r.watershed</a></em> can be used (in this case the topology has to be cleaned manually). More than one tree in the input data is allowed. No given flow direction is needed. To find the outlet of each tree, a DEM must be given.
-
-<h2>NOTES</h2>
-<h3>Problems</h3>
-This program is in beta status. It has the following shortcomings: <br />
-- The input data has to be topologically clean. Use <em>sloppy</em>=0.0 <br />
-- Source code comments are not doxydized. <br />
-
-
-<h3>How the algorithm works</h3>
-In a first step, <em>v.strahler</em> identifies all separate networks in the input dataset. That means, all connected lines are assigned a common <em>Basin ID</em>.<br />
-Consequently, the lowest leaf of each tree is identified as outlet. An auxiliary DEM is needed at this point. <br />
-Strahler ordering begins at each leaf of the tree with order N=1 (excluding the outlet). At a confluence, the order N(r) of the resulting stream is equal to the highest order N(max) of the joining streams or is raised by one if there are two or more joining streams of order N(max). <br />
-The algorithm returns an ASCII text file with columns: <em>Category</em> (from input map), <em>Line</em> (topology), <em>Basin</em>, <em>Order</em>; and the output map has the Strahler Order value instead of "category" for each line, and no connection with the database.
-
-<h2>EXAMPLE</h2>
-The input map (vector on DEM): <br />
-<br />
-<img src="images/example_input.jpg"> <br />
-<br />
-<br />
-An example of the sintax in GRASS shell: <br />
-<br />
-<div class="code"><pre>
-GRASS 6.3.0RC1 (ED50_Z33):~ > v.strahler input=esp_1 at mapset output=esp_1 dem=dem_20 at mapset
-txout=/home/mapset/esp_1 sloppy=0 layer=1 --overwrite 2>pippo
-</pre></div>
-<br />
-<br />
-An extract from the resultant text file: <br />
-<br />
-<div class="code"><pre>
-== Result of Strahler Order ==
- Category: Line: Basin: Order:
- 22 1 1 4
- 73 2 1 1
- 25 3 1 4
- 27 4 1 4
- 39 5 1 1
- 48 6 1 1
- 56 7 1 4
- 55 8 1 2
- 88 9 1 4
- 59 10 1 1
- 60 11 1 1
- 61 12 1 2
- 83 13 1 2
- 82 14 1 1
- 91 15 1 1
- 87 16 1 4
- 95 17 1 2
- 96 18 1 1
- 106 19 1 1
- 102 20 1 1
- 104 21 1 4
- 111 22 1 1
- 112 23 1 3
- 121 24 1 3
- 126 25 1 1
- 128 26 1 3
- 137 27 1 3
-</pre></div>
-<br />
-<br />
-The output map (vector on DEM) form a snapshot of GRASS display: <br />
-<br />
-<img src="images/example_output.jpg"> <br />
-<br />
-<br />
-The output map (vector only) from Qgis (different colors for different orders): <br />
-<br />
-<img src="images/Qgis.jpg"> <br />
-<br />
-<br />
-<h2>SEE ALSO</h2>
-<em><a href="r.watershed.html">r.watershed</a></em><BR>
-
-<h2>AUTHOR</h2>
-
-Florian Kindl, Univ. Innsbruck.<br />
-<br />
-Modified by: Ivan Marchesini and Annalisa Minelli, Univ. Perugia.<br />
-
-
-<p><i>Last changed: $Date: 2007/11/25 15:29:16 $</i>
Copied: grass-addons/vector/v.strahler/description.html (from rev 30226, grass-addons/v.strahler/description.html)
===================================================================
--- grass-addons/vector/v.strahler/description.html (rev 0)
+++ grass-addons/vector/v.strahler/description.html 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,87 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.strahler</em> calculates the Strahler Order for all lines of a given dendritic network. The input vector map must be free of cycles. For the elaboration a new imported network or a network extracted from DEM by using <em><a HREF="r.watershed.html">r.watershed</a></em> can be used (in this case the topology has to be cleaned manually). More than one tree in the input data is allowed. No given flow direction is needed. To find the outlet of each tree, a DEM must be given.
+
+<h2>NOTES</h2>
+<h3>Problems</h3>
+This program is in beta status. It has the following shortcomings: <br />
+- The input data has to be topologically clean. Use <em>sloppy</em>=0.0 <br />
+- Source code comments are not doxydized. <br />
+
+
+<h3>How the algorithm works</h3>
+In a first step, <em>v.strahler</em> identifies all separate networks in the input dataset. That means, all connected lines are assigned a common <em>Basin ID</em>.<br />
+Consequently, the lowest leaf of each tree is identified as outlet. An auxiliary DEM is needed at this point. <br />
+Strahler ordering begins at each leaf of the tree with order N=1 (excluding the outlet). At a confluence, the order N(r) of the resulting stream is equal to the highest order N(max) of the joining streams or is raised by one if there are two or more joining streams of order N(max). <br />
+The algorithm returns an ASCII text file with columns: <em>Category</em> (from input map), <em>Line</em> (topology), <em>Basin</em>, <em>Order</em>; and the output map has the Strahler Order value instead of "category" for each line, and no connection with the database.
+
+<h2>EXAMPLE</h2>
+The input map (vector on DEM): <br />
+<br />
+<img src="images/example_input.jpg"> <br />
+<br />
+<br />
+An example of the sintax in GRASS shell: <br />
+<br />
+<div class="code"><pre>
+GRASS 6.3.0RC1 (ED50_Z33):~ > v.strahler input=esp_1 at mapset output=esp_1 dem=dem_20 at mapset
+txout=/home/mapset/esp_1 sloppy=0 layer=1 --overwrite 2>pippo
+</pre></div>
+<br />
+<br />
+An extract from the resultant text file: <br />
+<br />
+<div class="code"><pre>
+== Result of Strahler Order ==
+ Category: Line: Basin: Order:
+ 22 1 1 4
+ 73 2 1 1
+ 25 3 1 4
+ 27 4 1 4
+ 39 5 1 1
+ 48 6 1 1
+ 56 7 1 4
+ 55 8 1 2
+ 88 9 1 4
+ 59 10 1 1
+ 60 11 1 1
+ 61 12 1 2
+ 83 13 1 2
+ 82 14 1 1
+ 91 15 1 1
+ 87 16 1 4
+ 95 17 1 2
+ 96 18 1 1
+ 106 19 1 1
+ 102 20 1 1
+ 104 21 1 4
+ 111 22 1 1
+ 112 23 1 3
+ 121 24 1 3
+ 126 25 1 1
+ 128 26 1 3
+ 137 27 1 3
+</pre></div>
+<br />
+<br />
+The output map (vector on DEM) form a snapshot of GRASS display: <br />
+<br />
+<img src="images/example_output.jpg"> <br />
+<br />
+<br />
+The output map (vector only) from Qgis (different colors for different orders): <br />
+<br />
+<img src="images/Qgis.jpg"> <br />
+<br />
+<br />
+<h2>SEE ALSO</h2>
+<em><a href="r.watershed.html">r.watershed</a></em><BR>
+
+<h2>AUTHOR</h2>
+
+Florian Kindl, Univ. Innsbruck.<br />
+<br />
+Modified by: Ivan Marchesini and Annalisa Minelli, Univ. Perugia.<br />
+
+
+<p><i>Last changed: $Date: 2007/11/25 15:29:16 $</i>
Deleted: grass-addons/vector/v.strahler/forest2tree.c
===================================================================
--- grass-addons/v.strahler/forest2tree.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/forest2tree.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,110 +0,0 @@
-#include "strahler.h"
-
-/* helper functions for ForestToTrees:
- push_line: push number of adjacent line on stack
- pop_line: return number of top line on stack
- we could do this with some Vect_list but we don't need
- to handle entire lines, only their identifiers */
-
-#define STACK_SIZE 65535
-int sp = 0; /* stack position */
-int stack[STACK_SIZE]; /* generous stack size */
-
-void push_line( int line )
-{
- if ( sp < STACK_SIZE ) {
- stack[sp++] = line;
- } else {
- G_fatal_error( "ForestToTrees: stack full" );
- }
-}
-
-int pop_line( void )
-{
- if ( sp > 0 ) {
- return stack[--sp];
- } else {
- G_debug( 3, "ForestToTrees: stack empty, tree finished" );
- return 0;
- }
-}
-
-
-/* or write our own lfind: nexttree() */
-/* start at offset, walk array up until bsnid == 0 and return that offset */
-
-int nexttree( DBBUF *dbbuf, int offset, int upper )
-{
- while ((dbbuf[offset].bsnid != 0) && (offset < upper) ) {
- ++offset;
- G_debug( 4, "bsnid[%d] is %d", offset, dbbuf[offset].bsnid);
- }
- return offset;
-}
-
-/* Find trees, assign basinID */
-int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf )
-{
- int tree, tree_finished, forest_done; /* ID of tree (basin), processing status */
- int l, n, d, degr, node; /* iterators */
- int nlines, offset; /* number of lines, position in dbbuf */
- int cline, aline; /* currently harvested line, adjacent line */
- int lnodes[2]; /* nodes of cline */
-
-
- G_debug( 1, "reached StrahForestToTrees" );
-
- nlines = Vect_get_num_lines( In );
-
- forest_done = 0;
- tree = 1; /* we start with tree no. 1 */
- offset = cline = 1; /* we start at line 1 */
-
- while ( forest_done == 0 ) {
- G_debug( 3, "\nProcessing Tree %d", tree );
- tree_finished = 0; /* initialize status indicator for this tree */
- while ( tree_finished == 0 ) {
- /* procedure to get adjacent lines and push them on the stack */
- Vect_get_line_nodes( In, cline, &lnodes[0], &lnodes[1] );
- G_debug( 4, "nodes for line %d: %d %d", cline, lnodes[0], lnodes[1] );
-
- for ( n=0; n<=1; n++) { /* loop through array[2] of fnode and tnode */
- node = lnodes[n];
- /*degr = Vect_get_node_n_lines( In, node );*/
- degr = StrahGetDegr( In, node );
- G_debug( 4, "degr %d for node %d", degr, node );
- for (d = 0; d < degr; d++) {
- /*aline = abs( Vect_get_node_line( In, node, d ) );*/
- aline = abs( StrahGetNodeLine( In, node, d ) );
- if ( dbbuf[aline].bsnid == 0 ) { /* push line on stack only if not done yet */
- G_debug( 4, "pushing line %d on sp %d", aline, sp);
- push_line( aline );
- G_debug( 4, "pushed, sp is now %d", sp);
- dbbuf[aline].bsnid = tree; /* line is done when pushed on stack */
- dbbuf[aline].line = aline;
- }
- }
- }
-
- /* pop a line from stack */
- cline = (int)pop_line();
- G_debug( 4, "popped line %d", cline );
- /* tree is finished when pop_line returned 0 */
- if ( cline == 0 ) {
- tree_finished = 1;
- }
- }
-
- /* get first cline of next tree */
- offset = cline = nexttree( dbbuf, ++offset, nlines );
- G_debug( 3, "cline/offset for tree %d is %d", tree, cline );
- if ( offset == nlines ) {
- forest_done = 1; /* what if last line still is 0 - it is a tree consisting of this line only - it shall keep bsnid 0 */
- } else {
- tree++; /* continue with next tree */
- }
- }
-
- return tree;
-}
-
Copied: grass-addons/vector/v.strahler/forest2tree.c (from rev 30226, grass-addons/v.strahler/forest2tree.c)
===================================================================
--- grass-addons/vector/v.strahler/forest2tree.c (rev 0)
+++ grass-addons/vector/v.strahler/forest2tree.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,110 @@
+#include "strahler.h"
+
+/* helper functions for ForestToTrees:
+ push_line: push number of adjacent line on stack
+ pop_line: return number of top line on stack
+ we could do this with some Vect_list but we don't need
+ to handle entire lines, only their identifiers */
+
+#define STACK_SIZE 65535
+int sp = 0; /* stack position */
+int stack[STACK_SIZE]; /* generous stack size */
+
+void push_line( int line )
+{
+ if ( sp < STACK_SIZE ) {
+ stack[sp++] = line;
+ } else {
+ G_fatal_error( "ForestToTrees: stack full" );
+ }
+}
+
+int pop_line( void )
+{
+ if ( sp > 0 ) {
+ return stack[--sp];
+ } else {
+ G_debug( 3, "ForestToTrees: stack empty, tree finished" );
+ return 0;
+ }
+}
+
+
+/* or write our own lfind: nexttree() */
+/* start at offset, walk array up until bsnid == 0 and return that offset */
+
+int nexttree( DBBUF *dbbuf, int offset, int upper )
+{
+ while ((dbbuf[offset].bsnid != 0) && (offset < upper) ) {
+ ++offset;
+ G_debug( 4, "bsnid[%d] is %d", offset, dbbuf[offset].bsnid);
+ }
+ return offset;
+}
+
+/* Find trees, assign basinID */
+int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf )
+{
+ int tree, tree_finished, forest_done; /* ID of tree (basin), processing status */
+ int l, n, d, degr, node; /* iterators */
+ int nlines, offset; /* number of lines, position in dbbuf */
+ int cline, aline; /* currently harvested line, adjacent line */
+ int lnodes[2]; /* nodes of cline */
+
+
+ G_debug( 1, "reached StrahForestToTrees" );
+
+ nlines = Vect_get_num_lines( In );
+
+ forest_done = 0;
+ tree = 1; /* we start with tree no. 1 */
+ offset = cline = 1; /* we start at line 1 */
+
+ while ( forest_done == 0 ) {
+ G_debug( 3, "\nProcessing Tree %d", tree );
+ tree_finished = 0; /* initialize status indicator for this tree */
+ while ( tree_finished == 0 ) {
+ /* procedure to get adjacent lines and push them on the stack */
+ Vect_get_line_nodes( In, cline, &lnodes[0], &lnodes[1] );
+ G_debug( 4, "nodes for line %d: %d %d", cline, lnodes[0], lnodes[1] );
+
+ for ( n=0; n<=1; n++) { /* loop through array[2] of fnode and tnode */
+ node = lnodes[n];
+ /*degr = Vect_get_node_n_lines( In, node );*/
+ degr = StrahGetDegr( In, node );
+ G_debug( 4, "degr %d for node %d", degr, node );
+ for (d = 0; d < degr; d++) {
+ /*aline = abs( Vect_get_node_line( In, node, d ) );*/
+ aline = abs( StrahGetNodeLine( In, node, d ) );
+ if ( dbbuf[aline].bsnid == 0 ) { /* push line on stack only if not done yet */
+ G_debug( 4, "pushing line %d on sp %d", aline, sp);
+ push_line( aline );
+ G_debug( 4, "pushed, sp is now %d", sp);
+ dbbuf[aline].bsnid = tree; /* line is done when pushed on stack */
+ dbbuf[aline].line = aline;
+ }
+ }
+ }
+
+ /* pop a line from stack */
+ cline = (int)pop_line();
+ G_debug( 4, "popped line %d", cline );
+ /* tree is finished when pop_line returned 0 */
+ if ( cline == 0 ) {
+ tree_finished = 1;
+ }
+ }
+
+ /* get first cline of next tree */
+ offset = cline = nexttree( dbbuf, ++offset, nlines );
+ G_debug( 3, "cline/offset for tree %d is %d", tree, cline );
+ if ( offset == nlines ) {
+ forest_done = 1; /* what if last line still is 0 - it is a tree consisting of this line only - it shall keep bsnid 0 */
+ } else {
+ tree++; /* continue with next tree */
+ }
+ }
+
+ return tree;
+}
+
Deleted: grass-addons/vector/v.strahler/helper.c
===================================================================
--- grass-addons/v.strahler/helper.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/helper.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,77 +0,0 @@
-#include "strahler.h"
-
-/* Replacement for Vect_get_node_n_lines in sloppy mode
-Gets adjacent lines by coordinate */
-
-int StrahGetDegr( struct Map_info *In, int node )
-{
- extern double sloppy;
- int degr, onode, j;
- double x, y;
- BOUND_BOX *Box;
- struct ilist *list;
-
- if ( sloppy != 0.0 ) {
- Vect_get_node_coor( In, node, &x, &y, NULL );
-
- Box->N = y + sloppy;
- Box->S = y - sloppy;
- Box->W = x + sloppy;
- Box->E = x - sloppy;
-
- list = Vect_new_list();
-
-/*
- printf( "N %f\n", Box->N );
- printf( "S %f\n", Box->S );
- printf( "W %f\n", Box->W );
- printf( "E %f\n", Box->E );
-*/
-
- Vect_select_nodes_by_box( In, Box, list );
- G_debug( 3, " %d nodes selected", list->n_values );
- /* is always 0 -why? */
- for ( j = 0; j < list->n_values; j++ ) {
- onode = abs( list->value[j] );
- G_debug( 3, "List %d: %d\n",j , onode );
- }
-/*
-Vect_list_append ( StArcs, line );
-Vect_get_line_nodes ( &Map, line, &node1, &node2);
-Vect_list_append ( StNodes, node1 );
-Vect_list_append ( StNodes, node2 );
-*/
-
- } else {
- degr = Vect_get_node_n_lines( In, node );
- }
-
- G_debug( 4, "node %d degr %d sloppy %f\n", node, degr, sloppy );
- return degr;
-}
-
-int StrahGetNodeLine( struct Map_info *In, int node, int d ) {
-
- extern double sloppy;
- int aline;
-
- if ( sloppy != 0.0 ) {
- G_debug(4, "StrahGetNodeLine at node %d", node);
- /*
- get all nodes within sloppy from node -> from table StrahGetDegr has written
- get all lines for each found node -> Vect_get_node_line(In, foundnode, df)
- for ( fn=0; fn<foundnodes.length; fn++ ) {
- degrf = Vect_get_node_n_lines( In, fn );
- for (df = 0; df < degrf; df++) {
- foundline = abs( Vect_get_node_line( In, fn, df ) );
- addtolist(foundline);
- }
- }
- return line d in list of foundlines;
- */
- } else {
- /* method if junction is only 1 node */
- aline = abs( Vect_get_node_line( In, node, d ) );
- }
- return aline;
-}
Copied: grass-addons/vector/v.strahler/helper.c (from rev 30226, grass-addons/v.strahler/helper.c)
===================================================================
--- grass-addons/vector/v.strahler/helper.c (rev 0)
+++ grass-addons/vector/v.strahler/helper.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,77 @@
+#include "strahler.h"
+
+/* Replacement for Vect_get_node_n_lines in sloppy mode
+Gets adjacent lines by coordinate */
+
+int StrahGetDegr( struct Map_info *In, int node )
+{
+ extern double sloppy;
+ int degr, onode, j;
+ double x, y;
+ BOUND_BOX *Box;
+ struct ilist *list;
+
+ if ( sloppy != 0.0 ) {
+ Vect_get_node_coor( In, node, &x, &y, NULL );
+
+ Box->N = y + sloppy;
+ Box->S = y - sloppy;
+ Box->W = x + sloppy;
+ Box->E = x - sloppy;
+
+ list = Vect_new_list();
+
+/*
+ printf( "N %f\n", Box->N );
+ printf( "S %f\n", Box->S );
+ printf( "W %f\n", Box->W );
+ printf( "E %f\n", Box->E );
+*/
+
+ Vect_select_nodes_by_box( In, Box, list );
+ G_debug( 3, " %d nodes selected", list->n_values );
+ /* is always 0 -why? */
+ for ( j = 0; j < list->n_values; j++ ) {
+ onode = abs( list->value[j] );
+ G_debug( 3, "List %d: %d\n",j , onode );
+ }
+/*
+Vect_list_append ( StArcs, line );
+Vect_get_line_nodes ( &Map, line, &node1, &node2);
+Vect_list_append ( StNodes, node1 );
+Vect_list_append ( StNodes, node2 );
+*/
+
+ } else {
+ degr = Vect_get_node_n_lines( In, node );
+ }
+
+ G_debug( 4, "node %d degr %d sloppy %f\n", node, degr, sloppy );
+ return degr;
+}
+
+int StrahGetNodeLine( struct Map_info *In, int node, int d ) {
+
+ extern double sloppy;
+ int aline;
+
+ if ( sloppy != 0.0 ) {
+ G_debug(4, "StrahGetNodeLine at node %d", node);
+ /*
+ get all nodes within sloppy from node -> from table StrahGetDegr has written
+ get all lines for each found node -> Vect_get_node_line(In, foundnode, df)
+ for ( fn=0; fn<foundnodes.length; fn++ ) {
+ degrf = Vect_get_node_n_lines( In, fn );
+ for (df = 0; df < degrf; df++) {
+ foundline = abs( Vect_get_node_line( In, fn, df ) );
+ addtolist(foundline);
+ }
+ }
+ return line d in list of foundlines;
+ */
+ } else {
+ /* method if junction is only 1 node */
+ aline = abs( Vect_get_node_line( In, node, d ) );
+ }
+ return aline;
+}
Copied: grass-addons/vector/v.strahler/images (from rev 30226, grass-addons/v.strahler/images)
Deleted: grass-addons/vector/v.strahler/main.c
===================================================================
--- grass-addons/v.strahler/main.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/main.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,356 +0,0 @@
-/****************************************************************************
- *
- * MODULE: v.strahler
- *
- * AUTHOR(S): Florian Kindl, Norbert Pfeifer, The GRASS development team
- * Modified by: Ivan Marchesini <marchesini unipg.it> and
- Annalisa Minelli <annapatri tiscali.it>
- *
- * PURPOSE: Assign Strahler order to dendritic network
- *
- * COPYRIGHT: (C) 2006 by the Authors
- *
- * This program is free software under the
- * GNU General Public License (v2).
- * Read the file COPYING that comes with GRASS
- * for details.
- *
- ****************************************************************************/
-
-/*! \file main.c
-\brief Assign Strahler order to dendritic network
-\author Florian Kindl
-
-\todo Deal with poor topology
- \li implemement sloppy mode or
- \li force clean topology
-\todo ...
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <search.h>
-#include <time.h>
-#include <grass/gis.h>
-#include <grass/Vect.h>
-#include <grass/dgl.h>
-#include <grass/dbmi.h>
-#include <grass/glocale.h>
-#include "strahler.h"
-
-#define COPYBLAH 1
-#define DDFIELDBLAH 1
-#define BADDCOLUMNBLAH 1
-#define RITEDBLAH 1
-#define EBUG 1
-#define REATEDBLAH 1
-
-int ret, type; /* return and type related to Vect_read_line and Vect_write_line */
-int fdrast; /* file descriptor for raster file is int */
-int ntrees; /* number of trees in dataset, returned by StrahForestToTrees */
-double sloppy;
-
-int main(int argc, char **argv)
-{
- int line, node;
- char buf[1024];
- int nnodes, nlines, tlines;
- DBBUF *dbbuf;
- NODEV *nodev;
- struct Option *input, *output, *dem_opt, *txout_opt, *sloppy_opt, *field_opt;
- struct GModule *module;
- char *inputset, *demset;
- struct Map_info In, Out;
- struct Map_info Trees;
- FILE *txout;
- struct Cell_head window; /* for cell sampling */
-
- extern int ntrees; /* number of trees calculated by StrahForestToTrees */
- extern int fdrast;
- extern double sloppy;
-
- /* Attribute table (from v.net.path/path.c) */
-
- struct line_cats *Cats; /* introduced and initialized the structures line_cats and line_pnts */
- static struct line_pnts *Points;
- Points = Vect_new_line_struct ();
- Cats = Vect_new_cats_struct ();
- dbString sql;
- dbDriver *driver;
- dbColumn *column;
-
- struct field_info *Fi;
- int field;
-
- /* Initialize the GIS calls */
-
- G_gisinit (argv[0]) ;
-
- input = G_define_standard_option(G_OPT_V_INPUT);
- output = G_define_standard_option(G_OPT_V_OUTPUT);
-
- module = G_define_module();
- module->description = _("Strahler order");
-
- dem_opt = G_define_standard_option(G_OPT_R_INPUT);
- dem_opt->key = "dem";
- dem_opt->description = _("Underlying DEM");
-
- txout_opt = G_define_option();
- txout_opt->key = "txout";
- txout_opt->type = TYPE_STRING;
- txout_opt->required = NO;
- txout_opt->multiple = NO;
- txout_opt->gisprompt = "new_file,file,output";
- txout_opt->description = _("Path to ASCII file where results will be written");
-
-
- /* option sloppy for bad topology */
-
- sloppy_opt = G_define_option();
- sloppy_opt->key = "sloppy";
- sloppy_opt->type = TYPE_DOUBLE;
- sloppy_opt->required = NO;
- sloppy_opt->answer = "0.0";
- sloppy_opt->multiple = NO;
- sloppy_opt->description = _("Threshold for distance within different nodes are considered the same node - may not work");
-
-
- field_opt = G_define_standard_option(G_OPT_V_FIELD);
- field = atoi (field_opt->answer);
-
- if(G_parser(argc,argv)) exit (EXIT_FAILURE);
-
- Vect_check_input_output_name( input->answer, output->answer, GV_FATAL_EXIT );
-
- inputset = G_find_vector2 (input->answer, NULL);
- if ( inputset == NULL) {
- G_fatal_error (_("Could not find input input <%s>"), input->answer);
- }
-
- /* open datasets at topology level 2 */
-
- Vect_set_open_level(2);
- Vect_open_old (&In, input->answer, inputset);
-
- G_debug(1, "Input vector opened");
-
- /* Open new vector, make 3D if input is 3D */
-
- if (1 > Vect_open_new( &Out, output->answer, Vect_is_3d( &In ))) {
- Vect_close( &In );
- G_fatal_error("Failed opening output vector file %s", output->answer);
- }
-
- /* Open DEM */
-
- demset = G_find_cell2( dem_opt->answer, NULL);
- if ( demset == NULL ) {
- G_fatal_error ("DEM file %s not found", dem_opt->answer);
- }
- if (1 > (fdrast = G_open_cell_old( dem_opt->answer, demset ))) {
- G_fatal_error("Failed opening DEM file %s", dem_opt->answer);
- }
-
- /* Open text file for results if given */
-
- if ( txout_opt->answer ) {
- txout = fopen( txout_opt->answer, "w" );
- if ( txout == NULL )
- G_fatal_error( _("Cannot open file %s"), txout_opt->answer);
- }
-
- sloppy = atof( sloppy_opt->answer );
- /* printf("blah %f", sloppy*2.4 ); */
-
-
- /* write history */
- /*Vect_copy_head_data( &In, &Out );
- Vect_hist_copy( &In, &Out );
- Vect_hist_command( &Out );*/
-
-#ifdef COPYBLAH
- /* this works, we want to continue here and ADD two columns to the &Out */
- /* Copy input to output (from v.clean/main.c) */
- /*G_debug( 1, "Copying vector lines to Output" );
-
- /* This works for both level 1 and 2 */
- /*Vect_copy_map_lines ( &In, &Out ); */
- /* 0: copy all fields, else field number */
-#endif
-
-
-#ifdef ADDFIELDBLAH /* Create table (from v.net.path/path.c) */
- /* or: add fields to existing table? */
- /* *Map, field, field_name, type */
-
- Fi = Vect_default_field_info ( &Out, field, NULL, GV_1TABLE );
- Vect_map_add_dblink ( &Out, field, NULL, Fi->table, "cat", Fi->database, Fi->driver);
-
- driver = db_start_driver_open_database ( Fi->driver, Fi->database );
- driver = db_start_driver_open_database ( Fi->driver, Vect_subst_var(Fi->database, &Out) );
- /* from v.reclass/main.c :
- */
-
- if ( driver == NULL ) {
- G_fatal_error ( "Cannot open database %s by driver %s", Fi->database, Fi->driver );
- }
-
-
- /* store the statement to create a table with category, basin ID and strahler order of the arc */
- sprintf ( buf, "CREATE TABLE %s ( cat integer, bsnid integer, sorder integer )", Fi->table );
-
- /* from v.reclass/main.c :*/
- dbString stmt;
- db_init_string ( &stmt );
- db_set_string ( &stmt, buf);
-
- if (db_execute_immediate (driver, &stmt) != DB_OK ) {
- Vect_close (&Out);
- db_close_database_shutdown_driver ( driver );
- G_fatal_error ( "Cannot create table: %s", db_get_string (&stmt) );
- }
-
- /*
- sprintf ( buf, "ALTER TABLE %s ADD COLUMN cat integer", Fi->table );
- sprintf ( buf, "ALTER TABLE %s ADD COLUMN bsnid integer", Fi->table );
- sprintf ( buf, "ALTER TABLE %s ADD COLUMN sorder integer", Fi->table );
- */
-#endif
-
-#ifdef DBADDCOLUMNBLAH
- sprintf ( buf, "bsnid");
- db_set_string ( &sql, buf );
- column->columnName = sql;
- column->sqlDataType = 1;
- sprintf ( buf, "%s", Fi->table );
- /*db_set_string ( &sql, buf );*/
- /* why such complicated arguments for db_add_column???
- dbDriver, dbString, dbColumn*/
- if ( db_add_column( driver, &sql, column ) != DB_OK ) {
- db_close_database_shutdown_driver ( driver );
- G_fatal_error ( "Cannot add column column in table %s", Fi->table );
- }
-#endif
-
-
- /* execute statement */
-#ifdef CREATEDBLAH
- sprintf ( buf, "CREATE TABLE %s ( cat integer, bsnid integer, sorder integer )", Fi->table );
- db_set_string ( &sql, buf );
- G_debug ( 2, db_get_string ( &sql ) );
-
- if (db_execute_immediate (driver, &sql) != DB_OK ) {
- db_close_database_shutdown_driver ( driver );
- G_fatal_error ( "Cannot create table: %s", db_get_string ( &sql ) );
- } else printf("table created\n");
-
- if ( db_create_index2(driver, Fi->table, "cat" ) != DB_OK ) {
- G_warning ( "Cannot create index" );
- } else printf("index created\n");
-
- if (db_grant_on_table (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP|DB_PUBLIC ) != DB_OK ) {
- G_fatal_error ( "Cannot grant privileges on table %s", Fi->table );
- } else printf("privileges granted\n");
-#endif
-
- nnodes = Vect_get_num_nodes ( &In );
- nlines = Vect_get_num_lines ( &In );
- G_debug( 1, "Number of lines: %d", nlines);
- G_debug( 1, "Number of nodes: %d", nnodes);
-
- /* Create table to store ordering */
- dbbuf = (DBBUF *) G_malloc ( (nlines + 1) * ((int)sizeof (DBBUF)) );
- nodev = (NODEV *) G_malloc ( (nnodes + 1) * ((int)sizeof (NODEV)) );
-
- /* initialize properly */
- for ( line=1; line <= nlines; line++) {
- dbbuf[line].category = dbbuf[line].line = dbbuf[line].bsnid = dbbuf[line].sorder = 0;
- }
- for ( node=1; node <= nnodes; node++) {
- nodev[node].node = nodev[node].degree = nodev[node].visited = 0;
- }
-
- ntrees = (int)StrahForestToTrees( &In, &Out, dbbuf );
- G_debug( 1, "Number of trees: %d", ntrees);
-
- StrahFindLeaves( &In, dbbuf, nodev, ntrees, fdrast );
-
- G_debug(2, "dbbuf after FindLeaves:\nline\tbsnid\tsorder");
- for ( line=1; line<=nlines; line++) {
- sprintf( buf, "%d\t%d\t%d\n",dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder);
- G_debug(2, "%s", buf);
- }
- G_debug(2, "nodev after FindLeaves:\nnode\tdegree\tvisited");
- for ( node=1; node<=nnodes; node++) {
- sprintf( buf, "%d\t%d\t%d\n",(int)nodev[node].node, nodev[node].degree, nodev[node].visited);
- G_debug(2, "%s", buf);
- }
-
- StrahOrder( &In, dbbuf, nodev );
-
- G_debug(2, "dbbuf after StrahOrder:\nline\tbsnid\tsorder\n");
- for ( line=1; line<=nlines; line++) {
- sprintf( buf, "%d\t%d\t%d\n",dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder);
- G_debug(2, "%s", buf);
- }
- G_debug(2, "nodev after StrahOrder:\nnode\tdegree\tvisited\n");
- for ( node=1; node<=nnodes; node++) {
- sprintf( buf, "%d\t%d\t%d\n",(int)nodev[node].node, nodev[node].degree, nodev[node].visited);
- G_debug(2, "%s", buf);
- }
-
-#ifdef WRITEDBLAH
- /* and write DB records */
- db_begin_transaction( driver );
-
- for ( line=1; line<=nlines; line++) {
- sprintf (buf, "insert into %s values ( %d, %d, %d )", Fi->table, line, dbbuf[line].bsnid, dbbuf[line].sorder );
- db_set_string ( &sql, buf );
- /* G_debug ( 3, db_get_string ( &sql ) ); */
- if (db_execute_immediate (driver, &sql) != DB_OK ) {
- db_close_database_shutdown_driver ( driver );
- G_fatal_error ( "Insert new row: %s", db_get_string ( &sql ) );
- }
- }
-
- db_commit_transaction( driver );
-
- db_close_database_shutdown_driver( driver );
-#endif
-
- /* writing StrahOrder instead of category in the output */
-
- for ( line=1; line<=nlines; line++) {
- type=Vect_read_line(&In, Points, Cats, line);
- ret = Vect_cat_del (Cats, field);
- G_debug(4, "deleted categories, ret=%d", ret);
- Vect_cat_set (Cats, field, dbbuf[line].sorder);
- Vect_write_line ( &Out, type, Points, Cats );
- G_debug(4, "category written for line %d", line);
- }
-
- /* assign category values */
-
- for ( line=1; line<=nlines; line++) {
- dbbuf[line].category = Vect_get_line_cat ( &In, line, field );
- G_debug(4, "added category %d for line %d", dbbuf[line].category, dbbuf[line].line);
- }
-
- /* Write text file for results if given */
-
- if ( txout_opt->answer )
- StrahWriteToFile( dbbuf, nlines, txout );
-
- Vect_close( &In );
- G_close_cell( fdrast );
-
- Vect_build( &Out, stdout );
- Vect_close( &Out );
-
- if ( txout_opt->answer )
- fclose( txout );
-
- return (EXIT_SUCCESS);
-}
Copied: grass-addons/vector/v.strahler/main.c (from rev 30226, grass-addons/v.strahler/main.c)
===================================================================
--- grass-addons/vector/v.strahler/main.c (rev 0)
+++ grass-addons/vector/v.strahler/main.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,356 @@
+/****************************************************************************
+ *
+ * MODULE: v.strahler
+ *
+ * AUTHOR(S): Florian Kindl, Norbert Pfeifer, The GRASS development team
+ * Modified by: Ivan Marchesini <marchesini unipg.it> and
+ Annalisa Minelli <annapatri tiscali.it>
+ *
+ * PURPOSE: Assign Strahler order to dendritic network
+ *
+ * COPYRIGHT: (C) 2006 by the Authors
+ *
+ * This program is free software under the
+ * GNU General Public License (v2).
+ * Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ ****************************************************************************/
+
+/*! \file main.c
+\brief Assign Strahler order to dendritic network
+\author Florian Kindl
+
+\todo Deal with poor topology
+ \li implemement sloppy mode or
+ \li force clean topology
+\todo ...
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <search.h>
+#include <time.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/dgl.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include "strahler.h"
+
+#define COPYBLAH 1
+#define DDFIELDBLAH 1
+#define BADDCOLUMNBLAH 1
+#define RITEDBLAH 1
+#define EBUG 1
+#define REATEDBLAH 1
+
+int ret, type; /* return and type related to Vect_read_line and Vect_write_line */
+int fdrast; /* file descriptor for raster file is int */
+int ntrees; /* number of trees in dataset, returned by StrahForestToTrees */
+double sloppy;
+
+int main(int argc, char **argv)
+{
+ int line, node;
+ char buf[1024];
+ int nnodes, nlines, tlines;
+ DBBUF *dbbuf;
+ NODEV *nodev;
+ struct Option *input, *output, *dem_opt, *txout_opt, *sloppy_opt, *field_opt;
+ struct GModule *module;
+ char *inputset, *demset;
+ struct Map_info In, Out;
+ struct Map_info Trees;
+ FILE *txout;
+ struct Cell_head window; /* for cell sampling */
+
+ extern int ntrees; /* number of trees calculated by StrahForestToTrees */
+ extern int fdrast;
+ extern double sloppy;
+
+ /* Attribute table (from v.net.path/path.c) */
+
+ struct line_cats *Cats; /* introduced and initialized the structures line_cats and line_pnts */
+ static struct line_pnts *Points;
+ Points = Vect_new_line_struct ();
+ Cats = Vect_new_cats_struct ();
+ dbString sql;
+ dbDriver *driver;
+ dbColumn *column;
+
+ struct field_info *Fi;
+ int field;
+
+ /* Initialize the GIS calls */
+
+ G_gisinit (argv[0]) ;
+
+ input = G_define_standard_option(G_OPT_V_INPUT);
+ output = G_define_standard_option(G_OPT_V_OUTPUT);
+
+ module = G_define_module();
+ module->description = _("Strahler order");
+
+ dem_opt = G_define_standard_option(G_OPT_R_INPUT);
+ dem_opt->key = "dem";
+ dem_opt->description = _("Underlying DEM");
+
+ txout_opt = G_define_option();
+ txout_opt->key = "txout";
+ txout_opt->type = TYPE_STRING;
+ txout_opt->required = NO;
+ txout_opt->multiple = NO;
+ txout_opt->gisprompt = "new_file,file,output";
+ txout_opt->description = _("Path to ASCII file where results will be written");
+
+
+ /* option sloppy for bad topology */
+
+ sloppy_opt = G_define_option();
+ sloppy_opt->key = "sloppy";
+ sloppy_opt->type = TYPE_DOUBLE;
+ sloppy_opt->required = NO;
+ sloppy_opt->answer = "0.0";
+ sloppy_opt->multiple = NO;
+ sloppy_opt->description = _("Threshold for distance within different nodes are considered the same node - may not work");
+
+
+ field_opt = G_define_standard_option(G_OPT_V_FIELD);
+ field = atoi (field_opt->answer);
+
+ if(G_parser(argc,argv)) exit (EXIT_FAILURE);
+
+ Vect_check_input_output_name( input->answer, output->answer, GV_FATAL_EXIT );
+
+ inputset = G_find_vector2 (input->answer, NULL);
+ if ( inputset == NULL) {
+ G_fatal_error (_("Could not find input input <%s>"), input->answer);
+ }
+
+ /* open datasets at topology level 2 */
+
+ Vect_set_open_level(2);
+ Vect_open_old (&In, input->answer, inputset);
+
+ G_debug(1, "Input vector opened");
+
+ /* Open new vector, make 3D if input is 3D */
+
+ if (1 > Vect_open_new( &Out, output->answer, Vect_is_3d( &In ))) {
+ Vect_close( &In );
+ G_fatal_error("Failed opening output vector file %s", output->answer);
+ }
+
+ /* Open DEM */
+
+ demset = G_find_cell2( dem_opt->answer, NULL);
+ if ( demset == NULL ) {
+ G_fatal_error ("DEM file %s not found", dem_opt->answer);
+ }
+ if (1 > (fdrast = G_open_cell_old( dem_opt->answer, demset ))) {
+ G_fatal_error("Failed opening DEM file %s", dem_opt->answer);
+ }
+
+ /* Open text file for results if given */
+
+ if ( txout_opt->answer ) {
+ txout = fopen( txout_opt->answer, "w" );
+ if ( txout == NULL )
+ G_fatal_error( _("Cannot open file %s"), txout_opt->answer);
+ }
+
+ sloppy = atof( sloppy_opt->answer );
+ /* printf("blah %f", sloppy*2.4 ); */
+
+
+ /* write history */
+ /*Vect_copy_head_data( &In, &Out );
+ Vect_hist_copy( &In, &Out );
+ Vect_hist_command( &Out );*/
+
+#ifdef COPYBLAH
+ /* this works, we want to continue here and ADD two columns to the &Out */
+ /* Copy input to output (from v.clean/main.c) */
+ /*G_debug( 1, "Copying vector lines to Output" );
+
+ /* This works for both level 1 and 2 */
+ /*Vect_copy_map_lines ( &In, &Out ); */
+ /* 0: copy all fields, else field number */
+#endif
+
+
+#ifdef ADDFIELDBLAH /* Create table (from v.net.path/path.c) */
+ /* or: add fields to existing table? */
+ /* *Map, field, field_name, type */
+
+ Fi = Vect_default_field_info ( &Out, field, NULL, GV_1TABLE );
+ Vect_map_add_dblink ( &Out, field, NULL, Fi->table, "cat", Fi->database, Fi->driver);
+
+ driver = db_start_driver_open_database ( Fi->driver, Fi->database );
+ driver = db_start_driver_open_database ( Fi->driver, Vect_subst_var(Fi->database, &Out) );
+ /* from v.reclass/main.c :
+ */
+
+ if ( driver == NULL ) {
+ G_fatal_error ( "Cannot open database %s by driver %s", Fi->database, Fi->driver );
+ }
+
+
+ /* store the statement to create a table with category, basin ID and strahler order of the arc */
+ sprintf ( buf, "CREATE TABLE %s ( cat integer, bsnid integer, sorder integer )", Fi->table );
+
+ /* from v.reclass/main.c :*/
+ dbString stmt;
+ db_init_string ( &stmt );
+ db_set_string ( &stmt, buf);
+
+ if (db_execute_immediate (driver, &stmt) != DB_OK ) {
+ Vect_close (&Out);
+ db_close_database_shutdown_driver ( driver );
+ G_fatal_error ( "Cannot create table: %s", db_get_string (&stmt) );
+ }
+
+ /*
+ sprintf ( buf, "ALTER TABLE %s ADD COLUMN cat integer", Fi->table );
+ sprintf ( buf, "ALTER TABLE %s ADD COLUMN bsnid integer", Fi->table );
+ sprintf ( buf, "ALTER TABLE %s ADD COLUMN sorder integer", Fi->table );
+ */
+#endif
+
+#ifdef DBADDCOLUMNBLAH
+ sprintf ( buf, "bsnid");
+ db_set_string ( &sql, buf );
+ column->columnName = sql;
+ column->sqlDataType = 1;
+ sprintf ( buf, "%s", Fi->table );
+ /*db_set_string ( &sql, buf );*/
+ /* why such complicated arguments for db_add_column???
+ dbDriver, dbString, dbColumn*/
+ if ( db_add_column( driver, &sql, column ) != DB_OK ) {
+ db_close_database_shutdown_driver ( driver );
+ G_fatal_error ( "Cannot add column column in table %s", Fi->table );
+ }
+#endif
+
+
+ /* execute statement */
+#ifdef CREATEDBLAH
+ sprintf ( buf, "CREATE TABLE %s ( cat integer, bsnid integer, sorder integer )", Fi->table );
+ db_set_string ( &sql, buf );
+ G_debug ( 2, db_get_string ( &sql ) );
+
+ if (db_execute_immediate (driver, &sql) != DB_OK ) {
+ db_close_database_shutdown_driver ( driver );
+ G_fatal_error ( "Cannot create table: %s", db_get_string ( &sql ) );
+ } else printf("table created\n");
+
+ if ( db_create_index2(driver, Fi->table, "cat" ) != DB_OK ) {
+ G_warning ( "Cannot create index" );
+ } else printf("index created\n");
+
+ if (db_grant_on_table (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP|DB_PUBLIC ) != DB_OK ) {
+ G_fatal_error ( "Cannot grant privileges on table %s", Fi->table );
+ } else printf("privileges granted\n");
+#endif
+
+ nnodes = Vect_get_num_nodes ( &In );
+ nlines = Vect_get_num_lines ( &In );
+ G_debug( 1, "Number of lines: %d", nlines);
+ G_debug( 1, "Number of nodes: %d", nnodes);
+
+ /* Create table to store ordering */
+ dbbuf = (DBBUF *) G_malloc ( (nlines + 1) * ((int)sizeof (DBBUF)) );
+ nodev = (NODEV *) G_malloc ( (nnodes + 1) * ((int)sizeof (NODEV)) );
+
+ /* initialize properly */
+ for ( line=1; line <= nlines; line++) {
+ dbbuf[line].category = dbbuf[line].line = dbbuf[line].bsnid = dbbuf[line].sorder = 0;
+ }
+ for ( node=1; node <= nnodes; node++) {
+ nodev[node].node = nodev[node].degree = nodev[node].visited = 0;
+ }
+
+ ntrees = (int)StrahForestToTrees( &In, &Out, dbbuf );
+ G_debug( 1, "Number of trees: %d", ntrees);
+
+ StrahFindLeaves( &In, dbbuf, nodev, ntrees, fdrast );
+
+ G_debug(2, "dbbuf after FindLeaves:\nline\tbsnid\tsorder");
+ for ( line=1; line<=nlines; line++) {
+ sprintf( buf, "%d\t%d\t%d\n",dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder);
+ G_debug(2, "%s", buf);
+ }
+ G_debug(2, "nodev after FindLeaves:\nnode\tdegree\tvisited");
+ for ( node=1; node<=nnodes; node++) {
+ sprintf( buf, "%d\t%d\t%d\n",(int)nodev[node].node, nodev[node].degree, nodev[node].visited);
+ G_debug(2, "%s", buf);
+ }
+
+ StrahOrder( &In, dbbuf, nodev );
+
+ G_debug(2, "dbbuf after StrahOrder:\nline\tbsnid\tsorder\n");
+ for ( line=1; line<=nlines; line++) {
+ sprintf( buf, "%d\t%d\t%d\n",dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder);
+ G_debug(2, "%s", buf);
+ }
+ G_debug(2, "nodev after StrahOrder:\nnode\tdegree\tvisited\n");
+ for ( node=1; node<=nnodes; node++) {
+ sprintf( buf, "%d\t%d\t%d\n",(int)nodev[node].node, nodev[node].degree, nodev[node].visited);
+ G_debug(2, "%s", buf);
+ }
+
+#ifdef WRITEDBLAH
+ /* and write DB records */
+ db_begin_transaction( driver );
+
+ for ( line=1; line<=nlines; line++) {
+ sprintf (buf, "insert into %s values ( %d, %d, %d )", Fi->table, line, dbbuf[line].bsnid, dbbuf[line].sorder );
+ db_set_string ( &sql, buf );
+ /* G_debug ( 3, db_get_string ( &sql ) ); */
+ if (db_execute_immediate (driver, &sql) != DB_OK ) {
+ db_close_database_shutdown_driver ( driver );
+ G_fatal_error ( "Insert new row: %s", db_get_string ( &sql ) );
+ }
+ }
+
+ db_commit_transaction( driver );
+
+ db_close_database_shutdown_driver( driver );
+#endif
+
+ /* writing StrahOrder instead of category in the output */
+
+ for ( line=1; line<=nlines; line++) {
+ type=Vect_read_line(&In, Points, Cats, line);
+ ret = Vect_cat_del (Cats, field);
+ G_debug(4, "deleted categories, ret=%d", ret);
+ Vect_cat_set (Cats, field, dbbuf[line].sorder);
+ Vect_write_line ( &Out, type, Points, Cats );
+ G_debug(4, "category written for line %d", line);
+ }
+
+ /* assign category values */
+
+ for ( line=1; line<=nlines; line++) {
+ dbbuf[line].category = Vect_get_line_cat ( &In, line, field );
+ G_debug(4, "added category %d for line %d", dbbuf[line].category, dbbuf[line].line);
+ }
+
+ /* Write text file for results if given */
+
+ if ( txout_opt->answer )
+ StrahWriteToFile( dbbuf, nlines, txout );
+
+ Vect_close( &In );
+ G_close_cell( fdrast );
+
+ Vect_build( &Out, stdout );
+ Vect_close( &Out );
+
+ if ( txout_opt->answer )
+ fclose( txout );
+
+ return (EXIT_SUCCESS);
+}
Copied: grass-addons/vector/v.strahler/r.broscoe.sh (from rev 30226, grass-addons/v.strahler/r.broscoe.sh)
===================================================================
--- grass-addons/vector/v.strahler/r.broscoe.sh (rev 0)
+++ grass-addons/vector/v.strahler/r.broscoe.sh 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,301 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.broscoe
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+#
+# PURPOSE: Calculates waerden test and t test statistics
+# for some values of threshold, on a single basin
+# according to A. J. Broscoe theory (1959)
+#
+# COPYRIGHT: (C) 2008 by 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.
+#
+# REQUIREMENTS: you need R installed with "agricolae" package
+# You can install the package by starting R as root and typing
+# install.packages("agricolae")
+# This will install the package from a CRAN mirror on the internet.
+# For more information on R see the documentation available at
+# http://www.r-project.org/
+#
+# TODO: solve stability problems for low area threshold (cf. v.strahler)
+#############################################################################
+#%Module
+#% description: Calculates waerden test and t test statistics for some values of threshold, on a single basin according to A. J. Broscoe theory (1959)
+#% keywords: waerden.test,t.test,mean stream drop
+#%End
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: thresholds
+#% type: integer
+#% description: Threshold values to calculate statistics on (separated by <space>)
+#% required : yes
+#%end
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet
+#% required : yes
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet
+#% required : yes
+#%end
+#%option
+#% key: lt
+#% type: integer
+#% description: Lesser than (in meters), the program doesn't consider stream drops lesser than this
+#% required : yes
+#%END
+#%option
+#% key: result
+#% type: string
+#% gisprompt: new_file,file,output
+#% key_desc: name
+#% description: Resultant text file to write statistics into
+#% required: yes
+#%end
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+dem=$GIS_OPT_DEM
+thresholds=$GIS_OPT_THRESHOLDS
+xcoor=$GIS_OPT_XCOOR
+ycoor=$GIS_OPT_YCOOR
+lt=$GIS_OPT_LT
+result=$GIS_OPT_RESULT
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+echo $LOCATION_NAME
+
+#g.region rast=$dem
+
+res=`g.region -p | grep res | sed 1d | cut -f2 -d':' | tr -d ' ' `
+
+#g.region vect=$input
+
+#g.region n=n+$res s=s-$res e=e+$res w=w-$res
+
+for j in $thresholds
+do
+
+g.remove vect=ordered_$j
+rm orderedtxt_$j
+
+echo "initializing statistics for threshold value=$j"
+
+r.strahler dem=$dem xcoor=$xcoor ycoor=$ycoor thr=$j output=ordered_$j textoutput=orderedtxt_$j --overwrite
+
+rm one
+rm two
+rm three
+rm maxord
+
+`cat orderedtxt_$j | sed 2d |sed 1d > one
+sort -r -k 4,4 one > two
+head -1 two > three
+awk 'NR==1{print $4}' three > "maxord"`
+
+maxord=`cat maxord`
+
+rm thrVSmaxord_$j
+
+echo "$j $maxord" > thrVSmaxord_$j
+
+ a=1
+ until [ "$a" -gt "$maxord" ]
+ do
+ g.remove vect=ord1,ord1_pl,ord1_pl_3d
+ g.remove vect=ord1_pl_3d_cat
+ rm tmp
+ rm ordcol
+ rm textout
+
+ v.extract input=ordered_$j output=ord1 type=line layer=1 new=-1 list=$a
+
+ v.build.polylines input=ord1 output=ord1_pl cats=no
+
+ v.drape input=ord1_pl type=line rast=$dem method=nearest output=ord1_pl_3d
+
+ v.category input=ord1_pl_3d output=ord1_pl_3d_cat type=line option=add cat=1 layer=1 step=1
+
+ v.db.addtable map=ord1_pl_3d_cat table=ord1_pl_3d_cat layer=1 'columns=cat integer, sx double, sy double, sz double, ex double, ey double, ez double'
+
+ v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=start units=meters 'column=sx, sy, sz'
+
+ v.to.db map=ord1_pl_3d_cat type=line layer=1 qlayer=1 option=end units=meters 'column=ex, ey, ez'
+
+ v.rast.stats vector=ord1_pl_3d_cat raster=$dem colprefix=stats percentile=90 --verbose -c
+
+
+ for i in `db.select table=ord1_pl_3d_cat database=$GISDBASE/$LOCATION_NAME/$MAPSET/dbf/ driver=dbf 'sql=select cat from ord1_pl_3d_cat' | sed 1d`
+ do
+
+ min=`echo "select stats_min from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ sz=`echo "select sz from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ ez=`echo "select ez from ord1_pl_3d_cat where cat=$i" | db.select | sed 1d`
+ if [ `echo "$ez - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ] || [ `echo "$sz - $min" | bc -l | cut -f1 -d'.'` -lt "$lt" ]
+ then
+ if [ "$ez" -gt "$sz" ]
+ then
+ echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ else
+ echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ fi
+ else
+ echo "$ez - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ echo "$sz - $min" | bc -l | cut -f1 -d'.' >> tmp
+ echo $a >> ordcol
+ fi
+ done
+ paste ordcol tmp > textout\_$a
+ #read
+
+ if [ "$a" -gt 1 ]
+ then
+ b=`echo "$a-1" | bc -l`
+ cat textout\_$b textout\_$a > textout_temp
+ rm textout\_$a textout\_$b
+ mv textout_temp textout\_$a
+ fi
+
+ a=$(($a+1))
+ done
+
+ mv textout\_$maxord textout
+ echo "plotting textout..."
+ cat textout
+
+ echo "
+ imported=read.csv('textout', sep = '\t', dec='.', header = F)
+ dframe=data.frame(imported)
+ str(dframe)
+ nrow_maxord=nrow(dframe[dframe\$V1==max(dframe\$V1),])
+ maxord=max(dframe\$V1)
+ ifelse(nrow_maxord==1, (dframe=na.omit(subset(dframe[dframe\$V1 < maxord,]))), (dframe))
+ cat('plotted dframe \n')
+ nord=c(1:max(dframe\$V1))
+ #search for negative values of drop, replacing them with 0.1 values
+ dframe\$V2[(dframe\$V2 < 0)]<-0.1
+ tmp<-dframe
+ #search for zero values of drop, replacing them with 0.1 values
+ tmp\$V2[dframe\$V2==0]<-0.1
+ dframe<-tmp
+ delta=c()
+ for(i in nord)
+ {
+ offset=(qt(0.025,nrow(dframe[dframe\$V1==i,])-1,lower.tail=FALSE)*(sd(log(dframe\$V2[dframe\$V1==i])))/sqrt(nrow(dframe[dframe\$V1==i,])))/(mean (log(dframe\$V2[dframe\$V1==i])))
+ delta=c(delta,offset)
+
+ }
+ delta=delta[delta<0.5]
+ delta<-data.frame(na.omit(delta))
+ dframe<-data.frame(na.omit(dframe[dframe\$V1 <= (nrow(delta)),]))
+ print(dframe)
+ cat('plotted input dataframe \n')
+ linear_regr=summary(lm(dframe\$V2~dframe\$V1))
+ print(linear_regr)
+ cat('plotted linear regression output \n')
+ t=linear_regr\$coefficients[2,3]
+ Pr=linear_regr\$coefficients[2,4]
+ Radj=linear_regr\$adj.r.squared
+ library(agricolae)
+ wer_t=waerden.test(dframe\$V2,dframe\$V1,group=F)
+ print(wer_t)
+ cat('plotted Van der Waerden test output \n')
+ sink(file=\"tmpfile\")
+ waerden.test(dframe\$V2,dframe\$V1,group=F)
+ sink()
+ Pval=try(system(\"cat tmpfile | grep Pvalue | cut -f2 -d' '\", intern=TRUE))
+ system(\"rm tmpfile\")
+ Pval<-as.numeric(Pval)
+ out_row<-as.numeric(c(\"$j\",t,Pr,Radj,Pval))
+ cat('plotting values: threshold, t, Pr, Radj, Pval \n')
+ print(out_row)
+ write(out_row,file='results\_$j',ncolumns=5,sep='\t')
+ " > R_temp
+
+ echo 'source ("R_temp")' | R --vanilla --slave
+
+ rm R_temp
+
+done
+
+rm $result
+rm thrVSmaxord
+
+echo "threshold t Pr Radj Pvalue">$result
+echo "threshold maxord">thrVSmaxord
+
+for k in $thresholds
+do
+ cat results_$k >> $result
+ rm results_$k
+
+ cat thrVSmaxord_$k >> thrVSmaxord
+ rm thrVSmaxord_$k
+done
+
+echo "plotting results..."
+cat $result
+
+echo "
+thrVSmaxord=read.csv(\"thrVSmaxord\", sep = '\t', dec='.', header = F)
+names(thrVSmaxord)[1]=\"Threshold\"
+names(thrVSmaxord)[2]=\"Maxord\"
+xgraph=read.csv(\"$result\", sep = '\t', dec='.', header = F)
+names(xgraph)[1]=\"Threshold\"
+names(xgraph)[2]=\"t_statistic\"
+names(xgraph)[3]=\"Pr\"
+names(xgraph)[4]=\"R_squared_adj\"
+names(xgraph)[5]=\"Pvalue\"
+
+pdf('waerden_test.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pvalue\")], type=\"l\", col=\"red\", lwd=3, lty=1, ylab=\"Pvalue\", pch=1))
+dev.off()
+
+pdf('linear_regression.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\")], type=\"l\", col=\"green\", lwd=3, lty=1, ylab=\"Pr\", pch=1))
+dev.off()
+
+pdf('all_tests.pdf')
+print(matplot(xgraph\$Threshold, xgraph[, c(\"Pr\",\"Pvalue\")], type=\"b\", lwd=3, lty=1, ylab=\"(1)Pr, (2)Pvalue\"))
+dev.off()
+
+pdf('thrVSmaxord.pdf')
+print(matplot(thrVSmaxord\$Threshold, thrVSmaxord[, c(\"Maxord\")], type=\"l\", col=\"blue\", lwd=3, lty=1, ylab=\"Maxord\", pch=1))
+dev.off()
+
+" > R_temp
+
+echo 'source ("R_temp")' | R --vanilla --slave
+
+rm R_temp
+
+
Deleted: grass-addons/vector/v.strahler/r.strahler.sh
===================================================================
--- grass-addons/v.strahler/r.strahler.sh 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/r.strahler.sh 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,174 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: r.strahler.sh
-# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
-# PURPOSE: Create a vector map of strahler ordered streems of a single
-# basin, starting from a DEM.
-# The script extracts from a DEM the network (by using
-# r.watershed), converts the network to vector lines,
-# cleans the topology, lets you decide an outlet for your basin,
-# and executes v.strahler for the upstream basin.
-#
-# COPYRIGHT: (C) 2008 by 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.
-#
-# TODO: solve stability problems for low area threshold
-#############################################################################
-#%Module
-#% description: Create a vector map of strahler ordered streems of a single basin starting from a DEM
-#% keywords: strahler, streems, vector
-#%End
-#%option
-#% key: dem
-#% type: string
-#% key_desc: dem
-#% gisprompt: old,cell,raster
-#% description: Name of DEM raster map
-#% required : yes
-#%END
-#%option
-#% key: thr
-#% type: integer
-#% description: minimum size of exterior watershed basin
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name of streems ordered map created
-#% required : yes
-#%end
-#%option
-#% key: textoutput
-#% type: string
-#% gisprompt: new_file,file,output
-#% key_desc: name
-#% description: Name for output text file
-#% required: yes
-#%end
-#%option
-#% key: bkgrmap
-#% type: string
-#% key_desc: bkgrmap
-#% gisprompt: old,cell,raster
-#% description: Name of background map to plot
-#% required : no
-#%END
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-dem=$GIS_OPT_DEM
-thr=$GIS_OPT_THR
-output=$GIS_OPT_OUTPUT
-textoutput=$GIS_OPT_TEXTOUTPUT
-bkgrmap=$GIS_OPT_BKGRMAP
-
-#check presence of raster MASK, put it aside
-MASKFOUND=0
-eval `g.findfile element=cell file=MASK`
-if [ "$file" ] ; then
- g.message "Raster MASK found, temporarily disabled"
- g.rename MASK,"${TMPNAME}_origmask" --quiet
- MASKFOUND=1
-fi
-
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
-
-#uncomment this line if you need to remove some maps remaining from previous runs
-#g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
-
-#get resolution of DEM
-res=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
-
-#find raster streems and plot them
-r.watershed elevation=$dem threshold=$thr stream=rnetwork drainage=drainage --overwrite
-r.null map=rnetwork setnull=0
-r.thin input=rnetwork output=rnetwork_th iterations=200 --overwrite
-r.mapcalc "rnetwork_thin=rnetwork_th/rnetwork_th"
-d.erase
-d.rast map=$bkgrmap
-d.rast -o map=rnetwork_thin
-
-#ask the user for zooming and choosing the outlet cell
-#-------
-echo "
-
-
-Now, please zoom to the area where you want to put the outlet cross section
-
-
-"
-
-d.zoom
-
-echo "
-
-
-Now, please click on the cell representing the cross section
-
-
-"
-
-#prevent from choosing a wrong cell outside from the streems
-cat=2
-while [ "$cat" != "1" ]
-do
-result=`d.what.rast -t -1 rnetwork_thin`
-coor=`echo $result | cut -f1 -d ' '`
-x=`echo $coor | cut -f1 -d':'`
-y=`echo $coor | cut -f2 -d':'`
-cat=`echo $result | cut -f3 -d ' ' | tr -d :`
-done
-
-#come back to the previous zoom
-d.zoom -r
-
-#find the basin and set the mask
-r.water.outlet drainage=drainage basin=basin easting=$x northing=$y
-d.rast basin
-r.null map=basin setnull=0
-r.mapcalc "MASK=basin"
-d.rast MASK
-
-#clean raster map and convert it to vector
-g.region rast=MASK
-r.to.vect input=rnetwork_thin output=vnetwork feature=line --overwrite
-d.vect map=vnetwork
-v.type input=vnetwork output=vnetwork_b type=line,boundary --overwrite
-v.centroids input=vnetwork_b output=vnetwork_b_add option=add layer=1 cat=1 step=1 --overwrite
-newres=`echo "$res/4" | bc -l`
-g.region res=$newres
-v.to.rast input=vnetwork_b_add output=rnetwork_b use=val type=area layer=1 value=1 rows=4096 --overwrite
-v.to.rast input=vnetwork output=rnetwork_l use=val type=line layer=1 value=1 rows=4096 --overwrite
-r.patch input=rnetwork_b,rnetwork_l output=rnetwork_patch --overwrite
-r.thin input=rnetwork_patch output=rnetwork_thin2 iterations=200 --overwrite
-r.to.vect input=rnetwork_thin2 output=vnetwork2 feature=line --overwrite
-
-#remove dangles and create polylines
-soglia=`echo "sqrt(2*($res^2))+1" | bc -l`
-v.clean input=vnetwork2 output=vnetwork_dangle type=line tool=rmdangle thresh=$soglia --overwrite
-v.build.polylines input=vnetwork_dangle output=vnetwork_poly cats=first --overwrite
-d.vect vnetwork_poly
-
-#enlarge region to ensure v.strahler get raster values at streems ends and use v.strahler
-g.region -a res=$res n=n+$soglia s=s-$soglia e=e+$soglia w=w-$soglia
-v.strahler input=vnetwork_poly output=$output dem=$dem sloppy=0 layer=1 txout=$textoutput
-
-#remove temporary files
-g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
-
Copied: grass-addons/vector/v.strahler/r.strahler.sh (from rev 30226, grass-addons/v.strahler/r.strahler.sh)
===================================================================
--- grass-addons/vector/v.strahler/r.strahler.sh (rev 0)
+++ grass-addons/vector/v.strahler/r.strahler.sh 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,221 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.strahler.sh
+# AUTHOR(S): Annalisa Minelli, Ivan Marchesini
+# PURPOSE: Create a vector map of strahler ordered streems of a single
+# basin, starting from a DEM.
+# The script extracts from a DEM the network (by using
+# r.watershed), converts the network to vector lines,
+# cleans the topology, lets you decide an outlet for your basin,
+# and executes v.strahler for the upstream basin.
+#
+# COPYRIGHT: (C) 2008 by 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.
+#
+# TODO: solve stability problems for low area threshold
+#############################################################################
+#%Module
+#% description: Create a vector map of strahler ordered streams of a single basin starting from a DEM
+#% keywords: strahler, streams, vector
+#%End
+#%flag
+#% key: i
+#% description: find interactively the outlet coords
+#%END
+#%option
+#% key: dem
+#% type: string
+#% key_desc: dem
+#% gisprompt: old,cell,raster
+#% description: Name of DEM raster map
+#% required : yes
+#%END
+#%option
+#% key: xcoor
+#% type: double
+#% description: x coord of outlet
+#% required : no
+#%end
+#%option
+#% key: ycoor
+#% type: double
+#% description: y coord of outlet
+#% required : no
+#%end
+#%option
+#% key: thr
+#% type: integer
+#% description: minimum size of exterior watershed basin
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Name of streams ordered map created
+#% required : yes
+#%end
+#%option
+#% key: textoutput
+#% type: string
+#% gisprompt: new_file,file,output
+#% key_desc: name
+#% description: Name for output text file
+#% required: yes
+#%end
+#%option
+#% key: bkgrmap
+#% type: string
+#% key_desc: bkgrmap
+#% gisprompt: old,cell,raster
+#% description: Name of background map to plot
+#% required : no
+#%END
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+dem=$GIS_OPT_DEM
+xcoor=$GIS_OPT_XCOOR
+ycoor=$GIS_OPT_YCOOR
+thr=$GIS_OPT_THR
+output=$GIS_OPT_OUTPUT
+textoutput=$GIS_OPT_TEXTOUTPUT
+bkgrmap=$GIS_OPT_BKGRMAP
+
+#check presence of x and y coordinates if you want to use the un-interactive mode
+
+if [ $GIS_FLAG_I -eq 0 ] && [ -z $xcoor ] && [ -z $ycoor ]; then
+echo "
+
+
+If you don't want to use the interactive-mode you must select the outlets' coodinates
+
+
+"
+exit 0
+
+fi
+
+#check presence of raster MASK, put it aside
+MASKFOUND=0
+eval `g.findfile element=cell file=MASK`
+if [ "$file" ] ; then
+ g.message "Raster MASK found, temporarily disabled"
+ g.rename MASK,"${TMPNAME}_origmask" --quiet
+ MASKFOUND=1
+fi
+
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
+
+#uncomment this line if you need to remove some maps remaining from previous runs
+#g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
+#get resolution of DEM
+res=`g.region -p | grep nsres | cut -f2 -d':' | tr -d ' '`
+
+#find raster streams and plot them
+r.watershed elevation=$dem threshold=$thr stream=rnetwork drainage=drainage --overwrite
+r.null map=rnetwork setnull=0
+r.thin input=rnetwork output=rnetwork_th iterations=200 --overwrite
+r.mapcalc "rnetwork_thin=rnetwork_th/rnetwork_th"
+d.erase
+
+if [ "$bkgrmap" ]; then
+d.rast map=$bkgrmap
+fi
+
+
+d.rast -o map=rnetwork_thin
+
+#use the interactive mode
+if [ $GIS_FLAG_I -eq 1 ] ; then
+
+ #ask the user for zooming and choosing the outlet cell
+ #-------
+ echo "
+
+
+ Now, please zoom to the area where you want to put the outlet cross section
+
+
+ "
+
+ d.zoom
+
+ echo "
+
+
+ Now, please click on the cell representing the cross section
+
+
+ "
+
+ #prevent from choosing a wrong cell outside from the streams
+ cat=2
+ while [ "$cat" != "1" ]
+ do
+ result=`d.what.rast -t -1 rnetwork_thin`
+ coor=`echo $result | cut -f1 -d ' '`
+ x=`echo $coor | cut -f1 -d':'`
+ y=`echo $coor | cut -f2 -d':'`
+ cat=`echo $result | cut -f3 -d ' ' | tr -d :`
+ done
+
+ #come back to the previous zoom
+ d.zoom -r
+
+#do not use the interactive mode
+else
+ x=$xcoor
+ y=$ycoor
+fi
+
+
+#find the basin and set the mask
+r.water.outlet drainage=drainage basin=basin easting=$x northing=$y
+d.rast basin
+r.null map=basin setnull=0
+r.mapcalc "MASK=basin"
+d.rast MASK
+
+#clean raster map and convert it to vector
+g.region rast=MASK
+r.to.vect input=rnetwork_thin output=vnetwork feature=line --overwrite
+d.vect map=vnetwork
+v.type input=vnetwork output=vnetwork_b type=line,boundary --overwrite
+v.centroids input=vnetwork_b output=vnetwork_b_add option=add layer=1 cat=1 step=1 --overwrite
+newres=`echo "$res/4" | bc -l`
+g.region res=$newres
+v.to.rast input=vnetwork_b_add output=rnetwork_b use=val type=area layer=1 value=1 rows=4096 --overwrite
+v.to.rast input=vnetwork output=rnetwork_l use=val type=line layer=1 value=1 rows=4096 --overwrite
+r.patch input=rnetwork_b,rnetwork_l output=rnetwork_patch --overwrite
+r.thin input=rnetwork_patch output=rnetwork_thin2 iterations=200 --overwrite
+r.to.vect input=rnetwork_thin2 output=vnetwork2 feature=line --overwrite
+
+#remove dangles and create polylines
+soglia=`echo "sqrt(2*($res^2))+1" | bc -l`
+v.clean input=vnetwork2 output=vnetwork_dangle type=line tool=rmdangle thresh=$soglia --overwrite
+v.build.polylines input=vnetwork_dangle output=vnetwork_poly cats=first --overwrite
+d.vect vnetwork_poly
+
+#enlarge region to ensure v.strahler get raster values at streams ends and use v.strahler
+g.region -a res=$res n=n+$soglia s=s-$soglia e=e+$soglia w=w-$soglia
+v.strahler input=vnetwork_poly output=$output dem=$dem sloppy=0 layer=1 txout=$textoutput
+
+#remove temporary files
+g.remove rast=rnetwork,rnetwork_thin,rnetwork_b,rnetwork_l,rnetwork_patch,rnetwork_thin2,drainage,MASK,basin vect=vnetwork,vnetwork_b,vnetwork_b_add,vnetwork2,vnetwork_dangle
+
Deleted: grass-addons/vector/v.strahler/strahler.c
===================================================================
--- grass-addons/v.strahler/strahler.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/strahler.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,210 +0,0 @@
-#include "strahler.h"
-
-/* Find Leaves of tree and distinguish outlet from sources */
-int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast )
-{
- int nnodes, degr, aline, node, unode, dnode, tree, outlet;
- double x, y, z, z_init;
-
- int n;
-
- struct Cell_head window;
- INTERP_TYPE method = NEAREST;
-
- OUTLETS *outlets;
-
- nnodes = Vect_get_num_nodes( In );
-
- z = z_init = 27000000.0; /* is it safe to initialize lowest z with the height of Olympus Mons above Martian Datum in millimeters? */
-
- outlets = (OUTLETS *) G_malloc ( (ntrees + 1) * ((int)sizeof (OUTLETS)) );
-
- /* initialize outlet table properly */
- for ( n=1; n <= ntrees; n++ ) {
- outlets[n].z = z_init;
- outlets[n].leaf = 0;
- }
-
- G_debug( 1, "Reached StrahFindLeaves with %d trees", ntrees);
-
- G_get_window( &window );
- G_debug( 2, "window: N %f S %f W %f E %f\n", window.north,window.south,window.west,window.east);
-
- G_debug( 2, "%d nodes in map", nnodes);
- for ( node = 1; node <= nnodes; node++ ) {
- /*degr = Vect_get_node_n_lines( In, node );*/
- degr = StrahGetDegr( In, node );
- nodev[node].node = node;
- nodev[node].degree = degr;
-
- G_debug( 4, "node %d: degr=%d", node, degr );
-
- if ( degr == 1 ) {
- unode = node;
- nodev[node].visited = 1; /*set the "visited" parameter to 1... can't be less*/
- /*aline = abs( Vect_get_node_line( In, unode, 0 ) );*/
- aline = abs( StrahGetNodeLine( In, unode, 0 ) );
- /*Vect_get_line_nodes( In, aline, &unode, &dnode);*/
- /* these are not necessarily up-node and down-node but we need not distinguish here... */
- dbbuf[aline].sorder = 1;
- /* no, do not visit yet
- nodev[node].visited += 1;
- */
- /*
- nodev[unode].visited += 1;
- nodev[dnode].visited += 1;
- */
- /* get z value or DEM value under node here, remember lowest */
- if ( Vect_is_3d( In ) ) {
- Vect_get_node_coor( In, node, &x, &y, &z );
- } else {
- Vect_get_node_coor( In, node, &x, &y, NULL );
- z = (double)G_get_raster_sample( fdrast, &window, NULL, y, x, 0, method );
- }
-
- G_debug( 5, "fdrast=%d node=%d y=%f x=%f z=%f\n", fdrast, node, y, x, z );
-
- tree = dbbuf[aline].bsnid;
- if ( z < outlets[tree].z ) {
- outlets[tree].z = z;
- outlets[tree].leaf = aline;
- }
- }
- }
- /* reset outlets to 0 */
- for ( n=1; n<=ntrees; n++) {
- G_debug( 2, "outlet tree %d: line %d with %f", n, outlets[n].leaf, outlets[n].z );
- outlet = outlets[n].leaf;
- dbbuf[outlet].sorder = 0;
- /*
- Vect_get_line_nodes( In, outlet, &unode, &dnode);
- nodev[unode].visited = 0;
- nodev[dnode].visited = 0;
- */
- }
- return 1;
-}
-
-/* assign strahler order */
-int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev )
-{
- int nlines, line;
- int fnode, tnode, unode, dnode; /* from-, to-, up-, down- node */
- int degr, d, aline, aorder; /* find adjacent lines */
- int corder, norder, dorder, dline; /* assign order */
- int cline, rfinish; /* follow one line until this stop-condition is 1 */
-
- struct dbbuf;
- struct nodev;
-
- nlines = Vect_get_num_lines ( In );
-
- G_debug( 1, "reached StrahOrder");
-
- for (line = 1; line <= nlines; line++) {
- if (dbbuf[line].sorder == 1) { /* get lines of order 1 */
-
- cline = line; /* and start run downwards */
- rfinish = 0;
- while (rfinish == 0) {
- G_debug( 3, "reached line %d", cline);
-
- Vect_get_line_nodes( In, cline, &fnode, &tnode); /* and their nodes */
-
- /**** BAUSTELLE ****/
- /* check nodes - we do not rely on flow direction */
- if ( nodev[fnode].visited == 1 ) {
- unode = fnode;
- dnode = tnode;
- } else {
- unode = tnode;
- dnode = fnode;
- }
-
- /*
- if (nodev[fnode].visited == nodev[tnode].visited) {
- printf("Line %d: both nodes are visited or not - what now?\n", line);
- } else if (nodev[fnode].visited == 0) {
- node = fnode;
- } else {
- node = tnode;
- }
- */
-
- /*
- if ( (nodev[fnode].degree - nodev[fnode].visited) < (nodev[tnode].degree - nodev[tnode].visited) ) {
- unode = fnode;
- dnode = tnode;
- } else if ( (nodev[fnode].degree - nodev[fnode].visited) > (nodev[tnode].degree - nodev[tnode].visited) ) {
- unode = tnode;
- dnode = fnode;
- } else {
- G_message("Line %d: nodes have same diff(degree,visited) - what now?\n", cline);
- */
- /* ok, what now?
- - occurs when leaf meets node with norder==0 -> visit nodes and continue
- - or: is outlet -> assign order
- */
- /* continue? */
- /*
- printf("visiting tnode %d and fnode %d for line %d and break\n\n", tnode, fnode, cline);
- nodev[tnode].visited += 1;
- nodev[fnode].visited += 1;
- break;
- }
- */
-
- G_debug( 3, "unode is %d and dnode is %d", unode, dnode);
-
- /*degr = Vect_get_node_n_lines ( In, dnode ); use downward node */
- degr = StrahGetDegr( In, dnode );
-
- G_debug( 4, "deg for dnode %d is %d", dnode, degr );
-
- corder = dbbuf[cline].sorder; /* current order - result won't be less than that */
- norder = dorder = 0; /* no order - how many lines have none?, highest order of others */
- for (d = 0; d < degr; d++) {
- /*aline = abs( Vect_get_node_line ( In, dnode, d ) );*/
- aline = abs( StrahGetNodeLine( In, dnode, d ) );
- G_debug( 4, "line %d at dnode %d is %d", d, dnode, aline);
- if (aline != cline) { /* we don't need the line we come from */
- aorder = dbbuf[aline].sorder;
- if (aorder == 0) {
- norder++;
- dline = aline; /* downward line */
- } else if (aorder > dorder) {
- dorder = aorder; /* find highest order */
- }
- }
- }
-
- if (norder > 1 || norder == 0) { /* if (norder > 1: node indeterminate OR norder == 0 : subtree finished) */
- G_debug( 3, "Finished run at line %d because norder=%d", cline, norder);
- G_debug( 4, "visiting unode %d for line %d\n", unode, cline);
- nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
- rfinish = 1; /* finish this run */
-
- } else { /* or */
- if (dorder > corder) {
- ; /* assign highest order of alines */
- } else if (dorder < corder) {
- dorder=corder; /* or keep order of cline */
- } else {
- dorder++; /* or raise order by one */
- }
-
- dbbuf[dline].sorder = dorder; /* assign order */
-
- G_debug( 4, "visiting dnode %d and unode %d for line %d", dnode, unode, cline);
- nodev[dnode].visited += 1; /* visit dnode */
- nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
-
- cline = dline; /* and continue with this line */
-
- G_debug( 3, "StrahOrder for dline %d: %d", dline, dorder);
- }
- } /* while */
- }
- }
- return 1;
-}
Copied: grass-addons/vector/v.strahler/strahler.c (from rev 30226, grass-addons/v.strahler/strahler.c)
===================================================================
--- grass-addons/vector/v.strahler/strahler.c (rev 0)
+++ grass-addons/vector/v.strahler/strahler.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,210 @@
+#include "strahler.h"
+
+/* Find Leaves of tree and distinguish outlet from sources */
+int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast )
+{
+ int nnodes, degr, aline, node, unode, dnode, tree, outlet;
+ double x, y, z, z_init;
+
+ int n;
+
+ struct Cell_head window;
+ INTERP_TYPE method = NEAREST;
+
+ OUTLETS *outlets;
+
+ nnodes = Vect_get_num_nodes( In );
+
+ z = z_init = 27000000.0; /* is it safe to initialize lowest z with the height of Olympus Mons above Martian Datum in millimeters? */
+
+ outlets = (OUTLETS *) G_malloc ( (ntrees + 1) * ((int)sizeof (OUTLETS)) );
+
+ /* initialize outlet table properly */
+ for ( n=1; n <= ntrees; n++ ) {
+ outlets[n].z = z_init;
+ outlets[n].leaf = 0;
+ }
+
+ G_debug( 1, "Reached StrahFindLeaves with %d trees", ntrees);
+
+ G_get_window( &window );
+ G_debug( 2, "window: N %f S %f W %f E %f\n", window.north,window.south,window.west,window.east);
+
+ G_debug( 2, "%d nodes in map", nnodes);
+ for ( node = 1; node <= nnodes; node++ ) {
+ /*degr = Vect_get_node_n_lines( In, node );*/
+ degr = StrahGetDegr( In, node );
+ nodev[node].node = node;
+ nodev[node].degree = degr;
+
+ G_debug( 4, "node %d: degr=%d", node, degr );
+
+ if ( degr == 1 ) {
+ unode = node;
+ nodev[node].visited = 1; /*set the "visited" parameter to 1... can't be less*/
+ /*aline = abs( Vect_get_node_line( In, unode, 0 ) );*/
+ aline = abs( StrahGetNodeLine( In, unode, 0 ) );
+ /*Vect_get_line_nodes( In, aline, &unode, &dnode);*/
+ /* these are not necessarily up-node and down-node but we need not distinguish here... */
+ dbbuf[aline].sorder = 1;
+ /* no, do not visit yet
+ nodev[node].visited += 1;
+ */
+ /*
+ nodev[unode].visited += 1;
+ nodev[dnode].visited += 1;
+ */
+ /* get z value or DEM value under node here, remember lowest */
+ if ( Vect_is_3d( In ) ) {
+ Vect_get_node_coor( In, node, &x, &y, &z );
+ } else {
+ Vect_get_node_coor( In, node, &x, &y, NULL );
+ z = (double)G_get_raster_sample( fdrast, &window, NULL, y, x, 0, method );
+ }
+
+ G_debug( 5, "fdrast=%d node=%d y=%f x=%f z=%f\n", fdrast, node, y, x, z );
+
+ tree = dbbuf[aline].bsnid;
+ if ( z < outlets[tree].z ) {
+ outlets[tree].z = z;
+ outlets[tree].leaf = aline;
+ }
+ }
+ }
+ /* reset outlets to 0 */
+ for ( n=1; n<=ntrees; n++) {
+ G_debug( 2, "outlet tree %d: line %d with %f", n, outlets[n].leaf, outlets[n].z );
+ outlet = outlets[n].leaf;
+ dbbuf[outlet].sorder = 0;
+ /*
+ Vect_get_line_nodes( In, outlet, &unode, &dnode);
+ nodev[unode].visited = 0;
+ nodev[dnode].visited = 0;
+ */
+ }
+ return 1;
+}
+
+/* assign strahler order */
+int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev )
+{
+ int nlines, line;
+ int fnode, tnode, unode, dnode; /* from-, to-, up-, down- node */
+ int degr, d, aline, aorder; /* find adjacent lines */
+ int corder, norder, dorder, dline; /* assign order */
+ int cline, rfinish; /* follow one line until this stop-condition is 1 */
+
+ struct dbbuf;
+ struct nodev;
+
+ nlines = Vect_get_num_lines ( In );
+
+ G_debug( 1, "reached StrahOrder");
+
+ for (line = 1; line <= nlines; line++) {
+ if (dbbuf[line].sorder == 1) { /* get lines of order 1 */
+
+ cline = line; /* and start run downwards */
+ rfinish = 0;
+ while (rfinish == 0) {
+ G_debug( 3, "reached line %d", cline);
+
+ Vect_get_line_nodes( In, cline, &fnode, &tnode); /* and their nodes */
+
+ /**** BAUSTELLE ****/
+ /* check nodes - we do not rely on flow direction */
+ if ( nodev[fnode].visited == 1 ) {
+ unode = fnode;
+ dnode = tnode;
+ } else {
+ unode = tnode;
+ dnode = fnode;
+ }
+
+ /*
+ if (nodev[fnode].visited == nodev[tnode].visited) {
+ printf("Line %d: both nodes are visited or not - what now?\n", line);
+ } else if (nodev[fnode].visited == 0) {
+ node = fnode;
+ } else {
+ node = tnode;
+ }
+ */
+
+ /*
+ if ( (nodev[fnode].degree - nodev[fnode].visited) < (nodev[tnode].degree - nodev[tnode].visited) ) {
+ unode = fnode;
+ dnode = tnode;
+ } else if ( (nodev[fnode].degree - nodev[fnode].visited) > (nodev[tnode].degree - nodev[tnode].visited) ) {
+ unode = tnode;
+ dnode = fnode;
+ } else {
+ G_message("Line %d: nodes have same diff(degree,visited) - what now?\n", cline);
+ */
+ /* ok, what now?
+ - occurs when leaf meets node with norder==0 -> visit nodes and continue
+ - or: is outlet -> assign order
+ */
+ /* continue? */
+ /*
+ printf("visiting tnode %d and fnode %d for line %d and break\n\n", tnode, fnode, cline);
+ nodev[tnode].visited += 1;
+ nodev[fnode].visited += 1;
+ break;
+ }
+ */
+
+ G_debug( 3, "unode is %d and dnode is %d", unode, dnode);
+
+ /*degr = Vect_get_node_n_lines ( In, dnode ); use downward node */
+ degr = StrahGetDegr( In, dnode );
+
+ G_debug( 4, "deg for dnode %d is %d", dnode, degr );
+
+ corder = dbbuf[cline].sorder; /* current order - result won't be less than that */
+ norder = dorder = 0; /* no order - how many lines have none?, highest order of others */
+ for (d = 0; d < degr; d++) {
+ /*aline = abs( Vect_get_node_line ( In, dnode, d ) );*/
+ aline = abs( StrahGetNodeLine( In, dnode, d ) );
+ G_debug( 4, "line %d at dnode %d is %d", d, dnode, aline);
+ if (aline != cline) { /* we don't need the line we come from */
+ aorder = dbbuf[aline].sorder;
+ if (aorder == 0) {
+ norder++;
+ dline = aline; /* downward line */
+ } else if (aorder > dorder) {
+ dorder = aorder; /* find highest order */
+ }
+ }
+ }
+
+ if (norder > 1 || norder == 0) { /* if (norder > 1: node indeterminate OR norder == 0 : subtree finished) */
+ G_debug( 3, "Finished run at line %d because norder=%d", cline, norder);
+ G_debug( 4, "visiting unode %d for line %d\n", unode, cline);
+ nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
+ rfinish = 1; /* finish this run */
+
+ } else { /* or */
+ if (dorder > corder) {
+ ; /* assign highest order of alines */
+ } else if (dorder < corder) {
+ dorder=corder; /* or keep order of cline */
+ } else {
+ dorder++; /* or raise order by one */
+ }
+
+ dbbuf[dline].sorder = dorder; /* assign order */
+
+ G_debug( 4, "visiting dnode %d and unode %d for line %d", dnode, unode, cline);
+ nodev[dnode].visited += 1; /* visit dnode */
+ nodev[unode].visited += 1; /* visit unode, for the sake of completeness */
+
+ cline = dline; /* and continue with this line */
+
+ G_debug( 3, "StrahOrder for dline %d: %d", dline, dorder);
+ }
+ } /* while */
+ }
+ }
+ return 1;
+}
Deleted: grass-addons/vector/v.strahler/strahler.h
===================================================================
--- grass-addons/v.strahler/strahler.h 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/strahler.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,95 +0,0 @@
-#include <stdlib.h>
-#include <string.h>
-#include <search.h>
-#include <time.h>
-#include "grass/gis.h"
-#include "grass/Vect.h"
-#include "grass/dgl.h"
-#include "grass/dbmi.h"
-
-/*! \typedef DBBUF
- \brief Buffer to keep track of assigned tree and Strahler order for each line
-*/
-
-/*! \typedef NODEV
- \brief Buffer to keep track of valency for each node
-*/
-
-/*! \typedef OUTLETS
- \brief Buffer to determine lowest leaf of each tree
-*/
-
-typedef struct { /* collect lines in Strahler order */
- int category;
- int line; /* line number */
- int bsnid; /* basin ID */
- int sorder; /* Strahler order */
-} DBBUF;
-
-typedef struct { /* visited nodes */
- int node; /* node number */
- int degree; /* valency */
- int visited; /* is between 0 and degr(node) */
-} NODEV;
-
-typedef struct { /* keep track of lowest leaf (=outlet) of each tree */
- double z; /* z-value of leaf */
- int leaf; /* line-id of outlet leaf */
-} OUTLETS;
-
-/* in forest2tree.c */
-/*! \fn int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf );
- \brief Returns the number of trees in *In
- \param *In The input map
- \param *Out The output map
- \param *dbbuf The buffer table for line orders
-*/
-
-
-/* in strahler.c */
-/*! \fn int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast );
- \brief Identifies all leaves of each tree and the one that lies lowest
- \param *In The input map
- \param *dbbuf The buffer table for line orders
- \param *nodev The buffer table for node valency
- \param ntrees Number of trees in map
- \param fdrast File descriptor of DEM raster map
-*/
-
-/*! \fn int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev );
- \brief Calculates the Strahler order of each line
- \param *In The input map
- \param *dbbuf The buffer table for line orders
- \param *nodev The buffer table for node valency
-*/
-
-
-/* in write.c */
-/*! \fn int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout );
- \brief Writes ASCII representation of calculated order to file
- \param *dbbuf The buffer table for line orders
- \param nlines Number of lines in *In map
- \param *txout ASCII output file name
-*/
-
-
-/* in helper.c */
-/*! \fn int StrahGetDegr( struct Map_info *In, int node );
- \brief Get degree of node (for sloppy mode)
- \param *In The input map
- \param node Node number in *In
-*/
-
-/*! \fn int StrahNodeLine( struct Map_info *In, int node, int d );
- \brief Get lines connected in node (for sloppy mode)
- \param *In The input map
- \param node Node number in *In
- \param d n-th line connected to node
-*/
-
-int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf );
-int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast );
-int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev );
-int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout );
-int StrahGetDegr( struct Map_info *In, int node );
-int StrahGetNodeLine( struct Map_info *In, int node, int d );
Copied: grass-addons/vector/v.strahler/strahler.h (from rev 30226, grass-addons/v.strahler/strahler.h)
===================================================================
--- grass-addons/vector/v.strahler/strahler.h (rev 0)
+++ grass-addons/vector/v.strahler/strahler.h 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,95 @@
+#include <stdlib.h>
+#include <string.h>
+#include <search.h>
+#include <time.h>
+#include "grass/gis.h"
+#include "grass/Vect.h"
+#include "grass/dgl.h"
+#include "grass/dbmi.h"
+
+/*! \typedef DBBUF
+ \brief Buffer to keep track of assigned tree and Strahler order for each line
+*/
+
+/*! \typedef NODEV
+ \brief Buffer to keep track of valency for each node
+*/
+
+/*! \typedef OUTLETS
+ \brief Buffer to determine lowest leaf of each tree
+*/
+
+typedef struct { /* collect lines in Strahler order */
+ int category;
+ int line; /* line number */
+ int bsnid; /* basin ID */
+ int sorder; /* Strahler order */
+} DBBUF;
+
+typedef struct { /* visited nodes */
+ int node; /* node number */
+ int degree; /* valency */
+ int visited; /* is between 0 and degr(node) */
+} NODEV;
+
+typedef struct { /* keep track of lowest leaf (=outlet) of each tree */
+ double z; /* z-value of leaf */
+ int leaf; /* line-id of outlet leaf */
+} OUTLETS;
+
+/* in forest2tree.c */
+/*! \fn int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf );
+ \brief Returns the number of trees in *In
+ \param *In The input map
+ \param *Out The output map
+ \param *dbbuf The buffer table for line orders
+*/
+
+
+/* in strahler.c */
+/*! \fn int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast );
+ \brief Identifies all leaves of each tree and the one that lies lowest
+ \param *In The input map
+ \param *dbbuf The buffer table for line orders
+ \param *nodev The buffer table for node valency
+ \param ntrees Number of trees in map
+ \param fdrast File descriptor of DEM raster map
+*/
+
+/*! \fn int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev );
+ \brief Calculates the Strahler order of each line
+ \param *In The input map
+ \param *dbbuf The buffer table for line orders
+ \param *nodev The buffer table for node valency
+*/
+
+
+/* in write.c */
+/*! \fn int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout );
+ \brief Writes ASCII representation of calculated order to file
+ \param *dbbuf The buffer table for line orders
+ \param nlines Number of lines in *In map
+ \param *txout ASCII output file name
+*/
+
+
+/* in helper.c */
+/*! \fn int StrahGetDegr( struct Map_info *In, int node );
+ \brief Get degree of node (for sloppy mode)
+ \param *In The input map
+ \param node Node number in *In
+*/
+
+/*! \fn int StrahNodeLine( struct Map_info *In, int node, int d );
+ \brief Get lines connected in node (for sloppy mode)
+ \param *In The input map
+ \param node Node number in *In
+ \param d n-th line connected to node
+*/
+
+int StrahForestToTrees( struct Map_info *In, struct Map_info *Out, DBBUF *dbbuf );
+int StrahFindLeaves( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev, int ntrees, int fdrast );
+int StrahOrder( struct Map_info *In, DBBUF *dbbuf, NODEV *nodev );
+int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout );
+int StrahGetDegr( struct Map_info *In, int node );
+int StrahGetNodeLine( struct Map_info *In, int node, int d );
Deleted: grass-addons/vector/v.strahler/write.c
===================================================================
--- grass-addons/v.strahler/write.c 2008-02-18 07:54:48 UTC (rev 30225)
+++ grass-addons/vector/v.strahler/write.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -1,14 +0,0 @@
-#include "strahler.h"
-
-int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout )
-{
- int line;
-
- G_debug( 1, "Reached WriteToFile - writing %d lines", nlines);
- fprintf( txout, "== Result of Strahler Order ==\n" );
- fprintf( txout, " Category: Line: Basin: Order:\n" ); /*added cat column*/
- for ( line=1; line<=nlines; line++) {
- fprintf( txout, "%9d%8d%8d%8d\n", dbbuf[line].category, dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder );
- }
- return 1;
-}
Copied: grass-addons/vector/v.strahler/write.c (from rev 30226, grass-addons/v.strahler/write.c)
===================================================================
--- grass-addons/vector/v.strahler/write.c (rev 0)
+++ grass-addons/vector/v.strahler/write.c 2008-02-18 10:12:35 UTC (rev 30227)
@@ -0,0 +1,14 @@
+#include "strahler.h"
+
+int StrahWriteToFile( DBBUF *dbbuf, int nlines, FILE *txout )
+{
+ int line;
+
+ G_debug( 1, "Reached WriteToFile - writing %d lines", nlines);
+ fprintf( txout, "== Result of Strahler Order ==\n" );
+ fprintf( txout, " Category: Line: Basin: Order:\n" ); /*added cat column*/
+ for ( line=1; line<=nlines; line++) {
+ fprintf( txout, "%9d%8d%8d%8d\n", dbbuf[line].category, dbbuf[line].line, dbbuf[line].bsnid, dbbuf[line].sorder );
+ }
+ return 1;
+}
More information about the grass-commit
mailing list