[GRASS-SVN] r39313 - in grass-addons/raster: . r.xtent
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 28 10:02:13 EDT 2009
Author: benducke
Date: 2009-09-28 10:02:12 -0400 (Mon, 28 Sep 2009)
New Revision: 39313
Added:
grass-addons/raster/r.xtent/
grass-addons/raster/r.xtent/Makefile
grass-addons/raster/r.xtent/at_exit_funcs.c
grass-addons/raster/r.xtent/at_exit_funcs.h
grass-addons/raster/r.xtent/description.html
grass-addons/raster/r.xtent/geom.c
grass-addons/raster/r.xtent/geom.h
grass-addons/raster/r.xtent/globals.h
grass-addons/raster/r.xtent/gt_vector.c
grass-addons/raster/r.xtent/gt_vector.h
grass-addons/raster/r.xtent/main.c
grass-addons/raster/r.xtent/r.xtent.001.png
grass-addons/raster/r.xtent/r.xtent.002.png
grass-addons/raster/r.xtent/r.xtent.003.png
grass-addons/raster/r.xtent/r.xtent.004.png
grass-addons/raster/r.xtent/r.xtent.005.png
grass-addons/raster/r.xtent/r.xtent.006.png
grass-addons/raster/r.xtent/r.xtent.007.png
grass-addons/raster/r.xtent/tools.c
grass-addons/raster/r.xtent/tools.h
grass-addons/raster/r.xtent/xtent.c
grass-addons/raster/r.xtent/xtent.h
Log:
Added r.xtent
Added: grass-addons/raster/r.xtent/Makefile
===================================================================
--- grass-addons/raster/r.xtent/Makefile (rev 0)
+++ grass-addons/raster/r.xtent/Makefile 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.xtent
+
+LIBES = $(GISLIB) $(DATETIMELIB) $(VECTLIB) $(DBMILIB)
+
+DEPENDENCIES = $(RASTERDEP) $(VECTDEP) $(DBMIDEP) $(GISDEP)
+
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+EXTRA_INC = $(VECT_INC)
+
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/raster/r.xtent/at_exit_funcs.c
===================================================================
--- grass-addons/raster/r.xtent/at_exit_funcs.c (rev 0)
+++ grass-addons/raster/r.xtent/at_exit_funcs.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,74 @@
+/***************************************************************************
+ * at_exit_funcs.c
+ *
+ * Mon Apr 18 14:52:20 2005
+ * Copyright 2005 Benjamin Ducke
+ ****************************************************************************/
+
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "globals.h"
+#include "tools.h"
+
+
+void clean_tmpfiles ( void ) {
+
+ int i;
+ char tmp[5000];
+ static int done = 0;
+
+
+ /* make sure this gets only run once */
+ if ( done > 0 ) {
+ return;
+ }
+
+ done ++;
+
+ if ( PROGRESS )
+ fprintf ( stdout, "Deleting temporary maps.\n");
+
+ /* clean up temp files: delete temporary maps */
+ if ( ( DISTANCE == COST ) && (COSTMODE == ADHOC) ) {
+ for ( i=0; i<num_centers; i ++ ) {
+ sprintf ( tmp, "g.remove rast=%s", tmapnames[i] );
+ if ( !VERBOSE ) {
+ strcat ( tmp," --q" );
+ }
+ run_cmd ( tmp );
+ }
+ if ( make_error ) {
+ sprintf ( tmp, "g.remove rast=%s", temapname );
+ if ( !VERBOSE ) {
+ strcat ( tmp," --q" );
+ }
+ run_cmd ( tmp );
+ }
+ /* TODO: check, if all maps were successfully removed */
+ }
+
+}
+
Added: grass-addons/raster/r.xtent/at_exit_funcs.h
===================================================================
--- grass-addons/raster/r.xtent/at_exit_funcs.h (rev 0)
+++ grass-addons/raster/r.xtent/at_exit_funcs.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,6 @@
+#ifndef AT_EXIT_FUNCS_H_
+#define AT_EXIT_FUNCS_H_
+
+void clean_tmpfiles ( void );
+
+#endif /*AT_EXIT_FUNCS_H_*/
Added: grass-addons/raster/r.xtent/description.html
===================================================================
--- grass-addons/raster/r.xtent/description.html (rev 0)
+++ grass-addons/raster/r.xtent/description.html 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,265 @@
+<H2>DESCRIPTION</H2>
+
+<em>r.xtent</em> outputs a raster map layer representing the Voronoi diagram,
+weighted Voronoi diagram or a more complex territorial partitioning of space around points (centers)
+in a vector input map (option <em>centers=</em>), based on the XTENT formula.
+
+<H2>BACKGROUND</H2>
+<p>
+The algorithm is based on the simple XTENT formula, first published by Renfrew and Level in 1979:
+<p>
+<em>I = C^a - k*d</em>
+<p>
+In the formula above, <em>I</em> is the strength of influence that each centers has on any given location (i.e. raster map cell) in the current region.
+The basic idea is that each cell will be allocated to the center that scores the highest <em>I</em> for that cell.
+The magnitude of <em>I</em> at a cell <em>x</em> for each center is determined by two terms that are weighted against each other:
+center weight (or size) <em>C</em> and distance <em>d</em>. Obviously, a large center in close proximity will have the best chance to score the
+highest <em>I</em> (i.e. "dominate" a cell). But a very large center can also be dominant, even if it is farther away.
+The two coefficients <em>a</em> and <em>k</em> determine the balance between center size and distance.
+The importance of distance increases in a linear manner while the importance of size increases exponentially.
+Thus, larger centers will compete stronger in relation to smaller ones, even at an increased distance.
+<p>
+The result map (option <em>output=</em>) is of type CELL (integer) raster with
+each cell representing the unique ID (category value) of the center to whose territory it belongs.<br>
+<p>
+Under certain parameters (specifically: <em>k=1</em>, no <em>-s</em> flag given and straight-line distance used), <em>r.xtent</em>
+will produce the equivalent of a simple (constant <em>C=1</em> and <em>a=0.5</em>) or weighted (variable <em>C</em>) Voronoi diagram.
+Using a distance measurement derived from a cost surface raster will improve result accuracy.
+
+<h2>USAGE</h2>
+
+In the most simple case, using the default values <em>a</em>=0.5 and <em>k</em>=1.0 and keeping <em>C</em>=1.0
+(see below), the result will equal that of a voronoi diagram calculated with <em>v.voronoi</em>. But <em>r.xtent</em> offers many options for
+inclusion of center size (weight), movement cost, barriers, pathways and a priori constraints to provide
+more detailed, realistic models of territorial allocation processes.
+<p>
+By default, <em>r.xtent</em> behaves differently than the original model published by Renfrew and Level in that it will
+always produce a complete partitioning of space where each cell is assigned to the territory of exactly one center.
+<br>
+The original definition of the XTENT formula demands that <em>I>=0.0</em>. This can be enforced using the "strict" flag
+(<em>-s</em>) in which case cells for which no center produces an <em>I>=0.0</em> will be set to NULL in the result map.
+This allows modeling incomplete partitionings of space where one ore more cells may not be assigned to the territory of
+any center. Another way of achieving an incomplete partitioning, which offers a bit more
+control, is to specify a "reach" attribute for each center (see section on "CONSTRAINTS" below).
+<p>
+
+The <em>a</em> and <em>k</em> coefficients are set globally, but center weight <e>C</e> must be set individually for
+each center. This is done by specifying the name of a double type attribute in the input vector points map (option <em>c=</em>).
+The program will read each center's weight from that attribute. If no attribute name is given, then C will be set
+to constant "1.0" for every center in the input map.
+
+<h3>CATEGORIES, LABELS, COLORS AND REPORTS</h3>
+
+The basic output map is an integer (CELL) type raster map. Each cell stores the unique number (ID) that represents
+the most dominant center at that location. By default, the integer type attribute "cat", which should be present
+in any GRASS vector map will be used to assign category numbers to each center. As an alternative, any attribute
+in the input points vector map of type integer may be chosen to represent the category number (option <em>categories=</em>)
+, but the user must make sure that every value is unique!<br>
+If no category value is specified and "cat" (the default category attribute for GRASS vector maps) is not present,
+the new IDs will be generated automaticall starting at "1" and being increased by one for each center in the
+input map.<br>
+Output map analysis may be facilitated if in addition to this a string (text) type attribute is specified to
+label the output map (i.e. for providing center names; option <em>labels=</em>).
+<p>
+The color of centers' territories in the output map can be controlled by suppling an attribute of type string that
+is set to the format RRR:GGG:BBB to specify red, green and blue intensity, respectively (option <em>rgb_column=</em>).
+Intensities can be in the range 0 to 255 for each component. E.g. 125:125:125 is a neutral gray.
+By default, <em>r.xtent</em> looks for a coloring attribute named "GRASSRGB" in the input map, which is in
+accordance with the behavior of <em>d.vect</em>.
+<p>
+Category and color attributes have default names ("cat" and "GRASSRGB"), which will be used, if present.
+If the respective attributes are not found or not of the rigth type, a warning will be issued and defaults used
+(i.e., a new index starting at "1" and a random color table).
+<p>
+A summary of the model output can be written to an ASCII file using the <em>report=</em> option. To get an output
+that is better suited for import into a spreadsheet, set the "tabular" flag (<em>-t</em>).
+<p>
+The report will contain at most the following fields for each center
+(some of the fields will not be shown or set to "n.a.", depending on the model option and results):
+
+<ul>
+ <li>id: the category value of the center</li>
+ <li>name: the name of the current center</li>
+ <li>dominates: a list of other centers (ids and names) that lie within the territory of the current site</li>
+ <li>dominated by: id and name of the center in whose territory the current center lies</li>
+ <li>cells: the number of map cells within the current site's territory</li>
+ <li>area: current site's territory in square meters (does not work for non-world coordinate systems)</li>
+ <li>reach cost/dist: (cost) distance between current center and farthest cell within its territory</li>
+ <li>aggressor: number of times that the current center is second in influence to any other center</li>
+ <li>competitor: id and name of the most frequent competitor for this center (see later)</li>
+ <li>competitor strength avg/var/std: statistics for the competitor's strength within the current center's territory</li>
+</ul>
+
+<p>
+The "dominates" line is skipped in the tabular report output as it will potentially create very long fields.
+The "dominates" relationship can be vizualized using a vector map symbology based on the "dominated by" attribute.
+
+
+<h3>DISTANCE AND MOVEMENT COSTS</h3>
+
+One straight-forward step towards more realism is to replace the straight line measure (which is actually a geodetic distance
+in the case of a lat/long region) used in the distance term <em>k*d</em> with a cost of movement measure that considers the actual, physical effort of moving
+on the ground from one location (cell) to another. For this purpose, <em>r.xtent</em> uses a cost surface map for each center <em>c</em>
+in the input vector points map that records the cost of movement from each cell in the map to the cell containing <em>c</em>.
+For each center, the user must first calculate a cost surface using an appropriate tools such as GRASS' <em>r.cost</em>
+and <em>r.walk</em> modules. The name of the corresponding raster map must then be added as a text type attribute field in the input map.
+The <em>costs_att=</em> attribute is used to tell <em>r.xtent</em> which attribute to read the cost map names from.
+
+
+<h3>BOUNDARIES AND PATHWAYS</h3>
+
+Boundaries and pathways are topographic features (usually of line geometry) that block or facilitate movement. Examples for
+boundaries may be broad rivers, mountain ridges or built walls. Pathways could be natural passes or built roads. Exactly what qualifies
+as a boundary or pathway depends on the type of movement modeled in the territorial allocation process.
+<p>
+If boundaries and/or pathways are to be added to the cost surfaces then they should be "burned" into a friction surfaces for
+use with <em>r.walk</em>. The <em>r.burn.frict</em> module should be used for this purpose.
+<p>
+Pathways should be "burned" after boundaries if pathways across boundaries (such as bridges) need to be modeled.
+
+
+<h3>COMPETITORS</h3>
+
+Territorial allocation is not always a clear-cut case. Many different combinations of variables may lead to situations where a
+center <em>A</em> "wins" highest <em>I</em> only by a relatively small margin over another center <em>B</em>. In this context, we will refer to center <em>B</em>
+as the "competitor" and to the inverse of the margin between <em>A</em> and <em>B</em> as "competition strength". Thus, whenever there is reason to say
+that "center <em>B</em> is almost as influential as center <em>A</em>" at a location <em>x</em>, then the strength of competition will be high and <em>B</em> will
+be classified as the competitor at that location (cell). Strength of competition is normalized to the range 0..1 for the entire region,
+with "1.0" being the highest strength of competition achievable.
+<p>
+Both the competitor's ID and the competition strength can be mapped (options <em>comp_id=</em> and <em>comp_strength=</em>)
+and may form the basis for interpretations.
+It is possible, e.g., to map "conflict" zones between center <em>A</em> and any number of other centers by isolating the cells within the
+territory of <em>A</em> whose competition exceeds a threshold value (e.g. 0.5) and optionally filtering them by the IDs of one or more
+possible competitors.
+<p>
+Competitor ID and competitor's strength can only be calculated for models in which at least two territories share at least one boundary.
+<p>
+Alternatively, competition strength can be interpreted as an error measure, with higher competition in locations where the territorial
+allocation may be erroneous due to flawed input data.
+
+<h3>CONSTRAINTS</h3>
+
+Further modeling flexibility may be gained by including <em>a priori</em> information as model constraints.
+It is possible to define a maximum distance (or costs) for each center, beyond which it cannot exercise and influence,
+no matter how marginal. Use the <em>reach=</em> option to specify a double type attribute that stores a maximum
+"reach" for each center. Set it to "-1.0" for centers with unrestricted reach.
+<p>
+It is also possible to define a center <em>A</em> as being dominated (ruled) by another center <em>B</em>, which means that it cannot compete against
+its "ruler" and any territory dominated by <em>A</em> will be allocated to <em>B</em>. This is done by specifying an integer attribute
+in the input vector map which should point to the ID of the ruling center or to that of the center itself in case there is no "ruler" (option <em>ruler=</em>).
+
+<h2>HINTS</h2>
+
+Model parameters (<em>C</em>,<em>a</em> and <em>k</em>) will most likely have to be adjusted if running in strict mode,
+otherwise the output may be a map with only NULL cells (<em>r.xtent</em> should issue a warning if this is the case).
+This can happen easily if the site weights are quantified in a smaller range than the distances. Unfortunately,
+the two may happen to use completely different value ranges. E.g. population size might give consistently much
+smaller values than distance in meters.
+<br>
+The model's parameters can be hard to control. It is advisable to go with the defaults for <em>a</em> and <em>k</em> check if there are centers
+with an extreme difference in weight.
+<p>
+The competitor's ID map in combination with a threshold error value can be used to map border regions from the
+perspective of any center. The extent of border zones could indicate potential "trouble zones". Zones with high error are zones
+that are potentially hard to control for the dominant center. There may also be a close third, fourth etc. contestant, so care
+must be exercised with the interpretation.
+
+
+<h3>CAVEATS AND LIMITATIONS</h3>
+
+Attributes in the input vector map that are to be used in the model (e.g. center weight, category attribute, labeling attribute)
+must be complete. There has to be a value for each center. The program's behavior for missing values is currently undefined.
+<p>
+Keep in mind that the calculation of costs between any cell in the region and a center is done by calculating the globally optimal
+shortest path from that cell to the center on a cost surface (or straight line). This assumes total spatial information being available,
+to the entities "on the move", an assumption that may or may not hold in a real-world scenario.
+<p>
+Territories should not get fragmented but form contiguous spatial phenomena. There are extreme cases, however, in which some small degree
+of fragmentation may occur, e.g. in topographically complex situations.
+<p>
+
+
+<h3>COMPUTATIONAL EFFICIENCY</h3>
+
+The memory footprint of <em>r.xtent</em> is very low and CPU time will grow linearly with the number of centers.
+If "ruler" constraints are given, however, processing time will increase exponentially.
+
+<H2>EXAMPLES</h2>
+
+The following are examples of territorial models with increasing complexity.
+
+<p>
+
+A basic model with center weight approximated by population size and store in the attribute "size" in the input
+map "towns"
+<div class="code"><pre>
+r.xtent centers=towns output=territories c=size
+</pre></div>
+
+<p align="center"><img src="r.xtent.001.png" alt="BPN chart"></p>
+
+<p>
+
+The same model, but this time combined with an elevation map to get a more realistic cost of movement distance
+measure (attribute table field "map_names" stores the names of the cost surface maps):
+
+<div class="code"><pre>
+r.xtent centers=towns output=territories c=size costs_att=map_names
+</pre></div>
+
+<p align="center"><img src="r.xtent.002.png" alt="BPN chart"><img src="r.xtent.003.png" alt="BPN chart"></p>
+
+<p>
+
+Here, an additional vector input map of a river representing a natural boundary has been added to
+a friction map using <em>r.burn.frict</em> and used in the cost surface calculation with <em>r.walk</em>:
+
+
+<p align="center"><img src="r.xtent.004.png" alt="BPN chart"><img src="r.xtent.005.png" alt="BPN chart"></p>
+
+<p>
+
+Another variation, with a smaller boundary weight (read from an attribute stored in the vector map "river" when performing the <em>r.burn.frict</em> step)
+to model a river that is easier to cross:
+
+
+<p align="center"><img src="r.xtent.006.png" alt="BPN chart"></p>
+
+<p>
+
+Another basic model, but this time with the additional <em>a priori</em> information that two of the towns are ruled by
+the third. This relationship is set in the integer attribute "ruler" in the input map "towns". The value of
+"ruler" is set to the category value of the ruling town (or to that of the town itself):
+
+<div class="code"><pre>
+r.xtent centers=towns output=territories c=size ruler=ruler --o
+</pre></div>
+
+<p align="center"><img src="r.xtent.007.png" alt="BPN chart"></p>
+
+
+<h2>BUGS</h2>
+
+The program may create temporary maps in the current mapset. If it crashes unexpectedly, it might leave
+those (prefixed with "XTENT") in the user's mapset. Use <em>g.mremove</em> to delete them.
+
+<H2>SEE ALSO</H2>
+<p>
+<EM><A HREF="v.voronoi">v.voronoi</A></EM>,
+<EM><A HREF="d.vect.html">d.vect</A></EM>,
+<EM><A HREF="r.cost.html">r.walk</A></EM>,
+<EM><A HREF="r.walk.html">r.walk</A></EM>,
+<EM><A HREF="r.burn.frict.html">r.burn.frict</A></EM>
+<EM><A HREF="g.mremove.html">g.mremove</A></EM>
+
+<H2>REFERENCES</H2>
+
+Renfrew, C. and Level, E. V. 1979: Exploring Dominance: Predicting Polities from Centers. In Transformations:
+Mathematical Approaches to Culture Change, edited by C. Renfrew and K. L. Cooke (Academic Press, New York), pp. 145-166.
+
+<H2>AUTHORS</h2>
+
+Benjamin Ducke (benjamin.ducke AT oadigital.net)
+
+<p><i>Last changed: 2009/09/28</i>
Added: grass-addons/raster/r.xtent/geom.c
===================================================================
--- grass-addons/raster/r.xtent/geom.c (rev 0)
+++ grass-addons/raster/r.xtent/geom.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,127 @@
+/*
+ * PRIMITIVE GEOMETRY FUNCTIONS UTILIZED BY FUNCTIONS IN xtent.c
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+
+#include "globals.h"
+#include "gt_vector.h"
+
+
+/* returns the record number of the vector point record that is spatially closest
+ to the point specified by x and y
+*/
+int get_closest_point ( double x, double y, GVT_map_s *map ) {
+ long int i;
+ long int closest;
+ double d;
+ double min = 0.0;
+ int min_set = 0;
+
+ GVT_rewind ( map );
+ i = 0;
+ closest = 0;
+ while ( GVT_next ( map ) ) {
+ i ++;
+ d = GVT_get_2D_point_distance ( x, y, map );
+ if ( !min_set ) {
+ min = d;
+ min_set = 1;
+ } else {
+ if ( d < min ) {
+ min = d;
+ closest = i-1;
+ }
+ }
+ }
+
+ return ( closest );
+
+}
+
+
+/* returns distance of point x,y to current point in *map */
+double get_distance ( double x, double y, double *costs, GVT_map_s *map ) {
+
+ long int i;
+
+ if ( DISTANCE == STRAIGHT ) {
+ return ( GVT_get_2D_point_distance ( x, y, map ) );
+ }
+
+ if ( DISTANCE == COST ) {
+ i = GVT_get_current (map);
+ return ( costs[i] );
+ }
+
+ return ( 0.0 );
+}
+
+
+/* get maximum distance to sites in *map for all raster cells in current region */
+double get_max_distance_region ( GVT_map_s *map ) {
+ struct Cell_head window;
+
+ int row, col;
+ int nrows, ncols;
+ double x, y;
+
+ double dist, max_dist;
+
+ G_get_window (&window);
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+
+ if ( PROGRESS ) {
+ fprintf ( stdout, " Finding maximum distance for all centers.\n" );
+ fflush ( stdout );
+ }
+
+ max_dist = 0.0;
+ for (row=0; row < nrows; row ++) {
+ y = G_row_to_northing ( (double) row, &window );
+ for (col=0; col < ncols; col ++) {
+ x = G_col_to_easting ( (double) col, &window );
+ /* TODO: bottle neck: implement coordinate caching in gt_vector.c */
+ GVT_rewind ( map );
+ while ( GVT_next ( map ) ) {
+ dist = GVT_get_2D_point_distance ( x, y, map );
+ if ( dist > max_dist ) {
+ max_dist = dist;
+ }
+ }
+ }
+ if ( PROGRESS ) {
+ G_percent (row, (nrows-1), 2);
+ fflush ( stdout );
+ }
+ }
+
+ return ( max_dist );
+}
+
+
+/* find biggest C in all centers */
+double get_max_size ( double *C, GVT_map_s *map ) {
+
+ long int i;
+ double max_C;
+
+ /* TODO: get rid of *C when internal attribute caching is done in vect tools lib */
+ max_C = 0.0;
+ GVT_rewind ( map );
+ i = 0;
+ while ( GVT_next ( map ) ) {
+ i ++;
+ if ( C[i] > max_C ) {
+ max_C = C[i];
+ }
+ }
+
+ return ( max_C );
+}
+
Added: grass-addons/raster/r.xtent/geom.h
===================================================================
--- grass-addons/raster/r.xtent/geom.h (rev 0)
+++ grass-addons/raster/r.xtent/geom.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,11 @@
+#ifndef GEOM_H_
+#define GEOM_H_
+
+/* distance calculation helpers */
+int get_closest_point ( double x, double y, GVT_map_s *map );
+double get_distance ( double x, double y, double *costs, GVT_map_s *map );
+/* normalization helpers */
+double get_max_distance_region ( GVT_map_s *map );
+double get_max_size ( double *C, GVT_map_s *map );
+
+#endif /*GEOM_H_*/
Added: grass-addons/raster/r.xtent/globals.h
===================================================================
--- grass-addons/raster/r.xtent/globals.h (rev 0)
+++ grass-addons/raster/r.xtent/globals.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,179 @@
+/***************************************************************************
+ * globals.h
+ *
+ * Mon Apr 18 15:04:11 2005
+ * Copyright 2005 Benjamin Ducke
+ ****************************************************************************/
+
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+
+/*
+#define EXPERIMENTAL
+*/
+
+#ifndef _GLOBALS_H
+#define _GLOBALS_H
+
+/* put a
+#define LOCAL
+into main.c ! */
+
+#ifdef LOCAL
+#define EXTERN
+#else
+#define EXTERN extern
+#endif
+
+/* uncomment for latest GRASS 6.3.CVS features */
+#define GRASS64
+
+#define PROGVERSION 0.99
+#define PROGNAME "r.xtent"
+
+/* if DEBUG > 1 we will dump even more details */
+#define DEBUG 0
+
+/* enable progress display? */
+EXTERN int PROGRESS;
+
+/* additional verbosity? */
+EXTERN int VERBOSE;
+
+/* overwrite existing output map? */
+EXTERN int OVERWRITE;
+
+/* GEOGRAPHIC OR XY COORDINATES? */
+#define XY 0
+#define GEOGRAPHIC 1
+EXTERN int SYSTEM;
+
+/* STUFF FOR MOVEMENT COST CALCULATIONS */
+#define STRAIGHT 0
+#define COST 1
+EXTERN int DISTANCE;
+
+#define NONE 0
+#define BASENAME 1
+#define ATTNAME 2
+#define MAPNAME 3
+#define ADHOC 4
+#define PSEUDO 5
+EXTERN int COSTMODE;
+
+#define MAXFPVAL 2000000000.0 /* maximum FP value (e.g. for absolute boundaries) */
+
+#define LAMBDA_DEFAULT 1.0
+#define WEIGHT_DEFAULT 1000.0 /* default weight of boundaries */
+#define PATH_DEFAULT 0.0 /* default weight of pathways */
+
+
+/* ================================================================================
+
+ STRUCTURES AND TYPES
+
+===================================================================================*/
+
+
+
+/* stores report results for one center */
+struct report_struct {
+ int id;
+ char *name;
+ long int out_count;
+ double area;
+ double percentage;
+ double max_cost;
+ long int sec_count;
+ long int err_count;
+ double err_sum;
+ double err_avg;
+ double err_var;
+ double err_std;
+ int competitor;
+ int competitor_id;
+ double competitor_p;
+ long int aggressor;
+ double aggressor_p;
+ /* the following point to other centers */
+ int boss;
+ int boss_id;
+ char *boss_name;
+ int num_subjects;
+ int *subject;
+ int *subject_id;
+ char **subject_name;
+};
+
+
+/* ================================================================================
+
+ GLOBAL MODULE OPTIONS AND FLAGS
+
+===================================================================================*/
+
+
+/* by making this global, all func's will have access to module options and flags */
+EXTERN struct GModule *module;
+
+EXTERN struct
+{
+ struct Option *centers; /* known political centers */
+ struct Option *output; /* name of new raster map */
+ struct Option *report; /* name of ASCII report file */
+ struct Option *costs_att; /* att name for cost surfaces */
+ struct Option *labels; /* name of attribute in centers map used for raster category labels */
+ struct Option *rgbcol; /* column with color specs in format RRR:GGG:BBB */
+ struct Option *cats; /* name of indexing attribute for writing result map values */
+ struct Option *errors; /* raster map/attribute in to store normalized error term */
+ struct Option *second; /* raster map to store ID of center with second-highest I */
+ struct Option *maxdist; /* maximum dist/cost attribute for each center */
+ struct Option *ruler; /* int attribute that points to a dominating center */
+ struct Option *ally; /* str attribut with comma-separated indices to allied centers */
+ struct Option *c; /* read C from this attribute in centers map */
+ struct Option *k; /* global k */
+ struct Option *a; /* global a */
+}
+parm;
+EXTERN struct
+{
+ struct Flag *tabular; /* tabular output format for report */
+ struct Flag *strict; /* imposes I >= 0 restriction (original formula) */
+}
+flag;
+
+/* WE NEED TO KEEP THE FOLLOWING GLOBAL, IN ORDER TO BE ABLE
+ * TO EASILY SPLIT UP THE COMPLEX PROGRAM LOGICS IN THE MAIN
+ * ROUTINE AND TO PASS INFORMATION TO THE AT_EXIT CLEAN UP
+ * FUNCTION
+ */
+
+/* basic XTENT parameters */
+EXTERN double a,k;
+EXTERN double *C;
+EXTERN double *C_norm; /* normalized C */
+EXTERN int I, I_second;
+EXTERN int num_centers;
+
+/* for adhoc generation of cost surface maps */
+EXTERN char **costcmds;
+EXTERN char **tmapnames; /* array of temp map names (costs) */
+EXTERN char *temapname; /* temp map name (unnormalized error map) */
+EXTERN int make_pseudo;
+EXTERN int make_error;
+EXTERN int make_colors;
+
+#endif /* _GLOBALS_H */
Added: grass-addons/raster/r.xtent/gt_vector.c
===================================================================
--- grass-addons/raster/r.xtent/gt_vector.c (rev 0)
+++ grass-addons/raster/r.xtent/gt_vector.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,1103 @@
+/*
+
+
+TODO:
+
+NULL value handling:
+ GVT_get_no_cache currently returns NULL for values that are NULL or missing -- maybe differentiate
+ between these two cases`?
+ use GRASS NULL values for CELL and DCELL types as int and double NULL representations !
+ -> look at value.c for hints !
+
+ATTRIBUTE CHECKING:
+- a function to check numerical attributes for validity, range and NULL data
+- the same for text attributes
+- a function to read all values for an attribute into an array and return it with NULL-termination
+- a function that returns any numerical attribute as DOUBLE
+- a function that returns any numerical attribute as INT
+
+HANLDE OTHER TYPE OF ATTRIBUTES:
+ DATE type
+ ... (?)
+
+
+ATTRIBUTE CACHING:
+ All functions reading/checking attributes call stub funcs GVT_get_attr() and GVT_attr_exists().
+ Those two need to be modified ...
+
+from d.what.vect:
+
+fprintf( stdout, _("Layer: %d\ncategory: %d\n"), Cats->field[j], Cats->cat[j] );
+
+currently, map->field seems to point to the currently layer's db connection, after:
+map->field = Vect_get_field (&map->in_vect_map, map->vect_cats->field[ map->current_layer - 1 ]);
+
+For advanced analytical operations, we need to able to open a map in level 2!
+Any analytical function requiring this should check if that has already been done,
+ if not attempt to open level 2, if that fails: warn or abort.
+
+Speed up attribute operations: copy all attribute values into an internal array, upon first
+query for that attribute (cache status=dirty). Serve all queries from that cache in the future.
+Set cache status back to dirty, if sth. has changed in the attribute table or another record
+has been added, so that the cache will be rebuilt on next access.
+
+Use Vect_get_num_updated_lines() and Vect_get_updated_line() to determine if cache needs to
+be set to dirty.
+
+A function that will return the numerical value of an attribute as type double, as long as
+that attribute is any number format (integer, double, ...)
+
+*/
+
+/* split into low and high level functions (e.g. export new vector map from SQL selection */
+
+/* if layer number = -1, parse all layers (see Martin Landa's recent posting; do other GRASS
+ modules handle it this way, too? Should this be default when opening a map? */
+
+/* in addition to the global
+
+ int GVT_NO_REGION_CONSTRAINTS;
+ int GVT_NO_TYPE_CONSTRAINTS;
+ int GVT_NO_SQL_CONSTRAINTS;
+
+ we need to be able to set these for each map, individually, too!
+
+ */
+
+
+#define GVT_DEBUG 1
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <grass/glocale.h>
+
+#include <grass/gis.h>
+#include <grass/dbmi.h>
+#include <grass/Vect.h>
+
+#include "globals.h"
+#include "gt_vector.h"
+
+
+/*
+ Check if GVT has been initialized and abort, if not.
+*/
+int GVT_check_init ( void ) {
+ return (0);
+}
+
+
+/* Allocates memory for a vector map, keeping track of the amount
+ of mem allocated in a counter in the struct.
+*/
+void * GVT_malloc ( size_t size, GVT_map_s *map ) {
+ void *p;
+ char tmp [2048];
+
+ p = G_malloc ( size );
+ if ( p == NULL ) {
+ /* TODO: %i is not the right type */
+ sprintf ( tmp, "GVT: Failed to allocate %i bytes of memory.\n", size );
+ if ( GVT_FATAL_MEM_ERRORS )
+ G_fatal_error (_( tmp ));
+ G_warning (_( tmp ));
+ return (p);
+ }
+ map->mem_size = map->mem_size + size;
+ return (p);
+}
+
+
+/*
+
+ PUBLIC FUNCTIONS
+
+*/
+
+
+/*
+ Initializes a bunch of global variables to default values.
+ Always call this function once in your program before using
+ any other GVT functions!
+*/
+void GVT_init ( void ) {
+ GVT_NO_REGION_CONSTRAINTS = FALSE;
+ GVT_NO_TYPE_CONSTRAINTS = FALSE;
+ GVT_NO_SQL_CONSTRAINTS = FALSE;
+
+ GVT_NO_ATTRIBUTE_CACHING = FALSE;
+
+ GVT_FATAL_FILE_ERRORS = TRUE;
+ GVT_FATAL_TYPE_ERRORS = TRUE;
+ GVT_FATAL_MEM_ERRORS = TRUE;
+
+ GVT_INIT_DONE = TRUE;
+}
+
+
+/* returns a pointer to a new GVT_map_s structure */
+/* caller must not pre-allocate that pointer! */
+GVT_map_s *GVT_new_map ( void ) {
+ GVT_map_s *map;
+ char tmp [2048];
+
+ struct Cell_head window;
+
+ map = G_malloc ( sizeof (GVT_map_s));
+ if ( map == NULL ) {
+ sprintf ( tmp, "GVT: Failed to allocate %i bytes of memory for new map.\n",
+ sizeof (GVT_map_s) );
+ if ( GVT_FATAL_MEM_ERRORS )
+ G_fatal_error (_( tmp ));
+ G_warning (_( tmp ));
+ return (map);
+ }
+ map->mem_size = map->mem_size + sizeof (GVT_map_s);
+ map->open_level = -1;
+
+ /* TODO: keep track of number of bytes alloc'd */
+ map->vect_points = Vect_new_line_struct ();
+ map->vect_cats = Vect_new_cats_struct ();
+
+ /* -2 indicates that no records exist in this map yet */
+ map->current_record = -2;
+
+ /* DEFAULT CONSTRAINTS */
+ /* filter vector objects to current region */
+ G_get_window (&window);
+ Vect_set_constraint_region (&map->in_vect_map, window.north, window.south,
+ window.east, window.west, window.top, window.bottom);
+ map->north = window.north;
+ map->south = window.south;
+ map->east = window.east;
+ map->west = window.west;
+ map->top = window.top;
+ map->bottom = window.bottom;
+
+ /* initialize SQL query */
+ db_init_string ( &map->sql );
+ db_set_string ( &map->sql, GVT_EMPTY_SQL );
+ /* initialize string buffer */
+ db_init_string ( &map->str );
+
+ /*lat-long location? */
+ if ( G_projection () == 3 ) {
+ map->isLatLong = TRUE;
+ } else {
+ map->isLatLong = FALSE;
+ }
+
+ /* distance calculations will be needed for lots of analytical applications, so might as
+ well initialize them right now ...
+ */
+ G_begin_distance_calculations ();
+
+ /* TODO: set all attribute caches to 'dirty' */
+
+ return ( map );
+}
+
+
+/* returns number of free'd bytes */
+int GVT_free_map (GVT_map_s *map) {
+
+ int freed = 0;
+
+ /* TODO: keep track of number of bytes dealloc'd */
+
+ if ( map == NULL ) {
+ return (0);
+ }
+
+ /* TODO: free:
+ Vect_new_line_struct ();
+ Vect_new_cats_struct ();
+ */
+
+ free ( map );
+
+ /* TODO: compare amount free'd with map->mem_size */
+
+ return ( freed );
+};
+
+
+/* constraints can only be set or unset when openening a new map! */
+/* returns the open level of the map */
+/* TODO: open in level2 ? */
+/* TODO: let user choose open level and abort/warn if it cannot be reached */
+/* TODO: handle layers: let user choose layer, if 0: look for the first layer
+ * that has objects of the specified type (issue a warning) or abort,
+ * if none found.
+ */
+int GVT_open_map ( char *mapname, int type, GVT_map_s *map ) {
+ char *mapset;
+ char tmp [2048];
+ int i, dbcol;
+
+ dbColumn *column;
+ dbValue *dbvalue;
+ int ctype, sqltype, more;
+ char *colname;
+
+ /* DO NOT ALLOCATE ANYTHING IN THE *MAP STRUCT BETWEEN HERE ... */
+
+ /* TODO: check if this object has already been created (?) */
+ if ( map->open_level > -1 ) {
+ sprintf (tmp, "GVT: Vector map already loaded into this object.\n" );
+ if ( GVT_FATAL_FILE_ERRORS )
+ G_fatal_error (_(tmp));
+ G_warning (_( tmp ));
+ return (-1);
+ }
+
+ /* TODO: check if a valid type has been specified */
+ /* and exit, if GVT_FATAL_TYPE_ERRORS is set */
+
+ /* check if map exists and if there are > 0 vector points in
+ the specified layer */
+ mapset = G_find_vector2 (mapname, "");
+ if (mapset == NULL) {
+ sprintf (tmp, "GVT: Could not open vector map: %s.\n", mapname );
+ if ( GVT_FATAL_FILE_ERRORS )
+ G_fatal_error (_(tmp));
+ G_warning (_( tmp ));
+ return (-1);
+ }
+
+ /* if success opening the vector map: */
+ map->open_level = 0;
+
+ /* attempt to open input map at level 1 */
+ Vect_set_open_level (1);
+ if (1 > Vect_open_old (&map->in_vect_map, mapname, mapset)) {
+ Vect_close (&map->in_vect_map);
+ free (mapset);
+ map->open_level = -1; /* reset open level */
+ sprintf (tmp, "GVT: Could not open vector map for level 1 access: %s.\n", mapname );
+ if ( GVT_FATAL_FILE_ERRORS )
+ G_fatal_error (_(tmp));
+ G_warning (_( tmp ));
+ return (-1);
+ }
+ map->open_level = 1;
+
+
+ /* ... AND HERE ! */
+
+ /* TODO: override constraints, if appropriate global vars have been set */
+
+ /* DEFAULT CONSTRAINTS */
+ /* filter vector objects to current region and point types only */
+ Vect_set_constraint_region (&map->in_vect_map, map->north, map->south,
+ map->east, map->west, map->top, map->bottom);
+ Vect_set_constraint_type (&map->in_vect_map, type);
+ map->type_constraints = type;
+
+ map->type = type;
+
+ if ( Vect_is_3d ( &map->in_vect_map ) ) {
+ map->is3d = TRUE;
+ } else {
+ map->is3d = FALSE;
+ }
+
+ sprintf ( map->name, mapname );
+
+ /* calculate number of vector points in map and current region */
+ map->num_records = 0;
+ /* TODO: map->type = : this should not be necessary if type constraints are in effect! */
+ while ((map->type = Vect_read_next_line (&map->in_vect_map, map->vect_points, map->vect_cats)) > 0) {
+ map->num_records ++;
+ }
+ /* TODO: determine number of layers and repeat the object counting for all layers */
+ /* exit, if no layer has the type requested and GVT_FATAL_TYPE_ERRORS is set */
+ /* save number of layers with valid type and associated layer numbers */
+
+ /* check for attributes */
+ /* TODO: this needs to be done for all layers, as given by n_cats (?)
+ so make sure there are no mem leaks in this loop
+ */
+ map->num_layers = 0;
+ map->current_layer = 0;
+ /*
+
+ TODO: currently, this only connects to the first layer in the map!
+ also, there is no abort condition, if map->vect_cats->n_cats < 1, i.e. there is no DB connection !!!
+ -> QUESTION: is there always at least on layer (field) with cats int attribute?
+
+ The loop should look like:
+
+ for ( i = 0; i < map->vect_cats->n_cats; i++ ) {
+ ...
+ }
+
+ because n_cats has the number of DB connections for this map.
+
+ */
+ for ( i = 0; i < 1; i++ ) {
+ map->field = Vect_get_field (&map->in_vect_map, map->vect_cats->field[i]);
+ /* TODO: in which cases could field == NULL? what to do then? */
+ /* TODO: catch cases without a DB connection */
+ if ( map->field != NULL ) {
+ /* start DB connection for the current layer (field) */
+ map->driver = db_start_driver ( map->field->driver );
+ db_init_handle ( &map->handle );
+ db_set_handle ( &map->handle, map->field->database, NULL);
+ if (db_open_database(map->driver, &map->handle) != DB_OK){
+ /* TODO: DB access failed, anything else we need to do? */
+ db_shutdown_driver( map->driver );
+ G_fatal_error (_("GVT: Error accessing database layer.\n"));
+ } else {
+ /* the default SQL statement is: "select all records" (in current layer) */
+ sprintf ( map->sql_buf, "select * from %s where %s = %d ",
+ map->field->table,
+ map->field->key,
+ map->vect_cats->cat[i]);
+ db_set_string ( &map->sql, map->sql_buf );
+ /* now open a cursor to fetch data from the DB */
+ /* other GVT funcs can use this cursor directly to access the attribute data for this map */
+ db_open_select_cursor(map->driver, &map->sql, &map->cursor, DB_SEQUENTIAL);
+ map->table = db_get_cursor_table (&map->cursor);
+ db_fetch ( &map->cursor, DB_NEXT, &more );
+ map->dbncols = db_get_table_number_of_columns (map->table);
+ /* TODO: lots of the vars in here should not be part of GVT_map_s struct
+ many of them need to be tracked for each layer separately
+ */
+ for ( dbcol = 0; dbcol < map->dbncols; dbcol ++ ) {
+ /* TODO: save this info in GVT_map_s, so that other funcs that need to check
+ attribute table types and structure can access it !
+ This info includes the name and type of each attribute in the table.
+ */
+ column = db_get_table_column(map->table, dbcol);
+ sqltype = db_get_column_sqltype (column);
+ ctype = db_sqltype_to_Ctype(sqltype);
+ dbvalue = db_get_column_value(column);
+ db_convert_value_to_string( dbvalue, sqltype, &map->str);
+ colname = (char*) db_get_column_name (column);
+
+ if ( DEBUG > 1 ) {
+ fprintf (stderr, "ATT = %s\n", colname);
+ }
+ }
+
+ db_close_cursor(&map->cursor);
+
+ }
+ }
+ map->num_layers ++;
+ }
+
+ /* the DB connection will remain open for other GVT funcs to use !!! */
+
+ /* TODO: activate first layer of given type */
+ map->current_layer = 1;
+
+ /* TODO: */
+ /* init attribute caches for current layer (?) */
+
+ Vect_rewind (&map->in_vect_map);
+ map->current_record = -1;
+ //free (mapset);
+
+ return ( map->open_level );
+}
+
+
+int GVT_open_map_points ( char *mapname, int type, GVT_map_s *map ) {
+ return ( GVT_open_map ( mapname, GV_POINT, map ) );
+}
+
+
+void GVT_close_map ( GVT_map_s *map ) {
+ Vect_close (&map->in_vect_map);
+ map->current_record = -2;
+
+ /* close DBMI driver and table(s) and cursor(s) */
+ /* TODO: this needs to be update for multilayered maps !!!
+ AND CHECK if a table was opened, at all !!!
+ */
+ db_close_database( map->driver);
+ db_shutdown_driver( map->driver);
+
+ /* TODO: what else needs to be done? */
+
+}
+
+
+/* jump to first record */
+void GVT_first ( GVT_map_s *map ) {
+ Vect_rewind ( &map->in_vect_map );
+ Vect_read_next_line (&map->in_vect_map, map->vect_points, map->vect_cats);
+ map->current_record = 0;
+}
+
+
+/* jump to next record */
+/* returns:
+ 1 (TRUE) SUCCESS
+ 0 (FALSE) NO MORE RECORDS (EOF)
+*/
+int GVT_next ( GVT_map_s *map ) {
+ /* TODO: does this honor current type constraints? */
+ if ( Vect_read_next_line (&map->in_vect_map, map->vect_points, map->vect_cats) > 0 ) {
+ map->current_record ++;
+ return ( TRUE );
+ } else {
+ return ( FALSE );
+ }
+}
+
+
+/* jump to last record */
+void GVT_last ( GVT_map_s *map ) {
+ GVT_seek ( map->num_records-1, map );
+ map->current_record = map->num_records-1;
+}
+
+
+/* rewind database: jump to position -1 (immediately before first record)! */
+void GVT_rewind ( GVT_map_s *map ) {
+ Vect_rewind ( &map->in_vect_map );
+ map->current_record = -1;
+}
+
+
+/* rewind and refilter data using SQL string */
+/* first record has position 0 ! */
+/* last record has position map->num_records-1 ! */
+/* returns:
+ 1 (TRUE) SUCCESS
+ 0 (FALSE attempting to seek to illegal position
+*/
+int GVT_seek ( long int pos, GVT_map_s *map ) {
+ long int i;
+
+ /* TODO: check, if this work also with layers of mixed types! */
+ if ( pos < 0 ) {
+ return ( FALSE );
+ }
+ if ( pos > map->num_records-1 ) {
+ return ( FALSE );
+ }
+
+ Vect_rewind ( &map->in_vect_map );
+ for ( i = 0; i <= pos; i ++ ) {
+ Vect_read_next_line (&map->in_vect_map, map->vect_points, map->vect_cats);
+ }
+ map->current_record = i;
+ return ( TRUE );
+}
+
+
+/* returns ID (cat - 1) of currently active record */
+/* TODO: maybe return cat (current_record + 1)
+ is ID always = cat-1?
+*/
+long int GVT_get_current ( GVT_map_s *map ) {
+ return ( map->current_record );
+}
+
+
+void GVT_query () {
+}
+
+
+
+/*
+ This function checks for the existence of attribute *name value in the current *map record's attributes
+ in the currently active layers.
+ It bypasses the GVT attribute caching scheme, using SQL calls to get the attribute data
+ from the DB which *map is connected to
+ If an attribute *name exists, TRUE will be returned, FALSE otherwise.
+ Type checking can be done by setting type to GVT_DOUBLE|INT|STRING ...
+ Set type to -1 to accept any type of attribute.
+*/
+int GVT_attr_exists_no_cache ( char *name, int type, GVT_map_s *map ) {
+ int dbcol;
+ int ctype, sqltype, more;
+ /* these two are pointers to members of the attribute table
+ they need not be free'd ! */
+ dbColumn *column;
+ char *colname;
+
+ if ( map->field != NULL ) {
+ /* Run SQL statement to refresh attribute information */
+ /* TODO: Do we really have to go through all this in order to read the current record's attribute values? */
+ /* START */
+ sprintf ( map->sql_buf, "select * from %s where %s = %d ",
+ map->field->table,
+ map->field->key,
+ map->vect_cats->cat[ map->current_layer - 1 ]);
+ db_set_string ( &map->sql, map->sql_buf );
+ db_open_select_cursor(map->driver, &map->sql, &map->cursor, DB_SEQUENTIAL); /* this does not alloc any memory, so no need to clean up */
+ map->table = db_get_cursor_table (&map->cursor);
+ db_fetch ( &map->cursor, DB_NEXT, &more );
+ /* END of stuff to do for reading attribute values */
+
+ for ( dbcol = 0; dbcol < map->dbncols; dbcol ++ ) {
+ column = db_get_table_column(map->table, dbcol);
+ sqltype = db_get_column_sqltype (column);
+ ctype = db_sqltype_to_Ctype(sqltype);
+ colname = (char*) db_get_column_name (column);
+ if (!strcmp (colname,name)){
+ /* found an attribute of that name! */
+ if ( type == -1 ) {
+ /* no type checking required! */
+ db_close_cursor(&map->cursor);
+ return ( TRUE );
+ }
+ if (ctype == type) {
+ /* now we can break out of this loop! */
+ db_close_cursor(&map->cursor);
+ return ( TRUE );
+ }
+ }
+ }
+ db_close_cursor(&map->cursor);
+ }
+ return ( FALSE );
+}
+
+/*
+ STUB FUNCTION
+*/
+int GVT_attr_exists ( char *name, int type, GVT_map_s *map ) {
+ return GVT_attr_exists_no_cache ( name, type, map );
+}
+
+
+/*
+ The following functions check if an attribute of a specific name and type exists.
+*/
+int GVT_double_exists ( char *name, GVT_map_s *map ) {
+ return ( GVT_attr_exists (name, GVT_DOUBLE, map));
+}
+
+int GVT_int_exists ( char *name, GVT_map_s *map ) {
+ return ( GVT_attr_exists (name, GVT_INT, map));
+}
+
+int GVT_numeric_exists ( char *name, GVT_map_s *map ) {
+ if ( (GVT_attr_exists (name, GVT_INT, map)) || (GVT_attr_exists (name, GVT_DOUBLE, map)) ) {
+ return ( TRUE );
+ }
+ return ( FALSE );
+}
+
+int GVT_string_exists ( char *name, GVT_map_s *map ) {
+ return ( GVT_attr_exists (name, GVT_STRING, map));
+}
+
+int GVT_any_exists ( char *name, GVT_map_s *map ) {
+ return ( GVT_attr_exists (name, -1, map));
+}
+
+
+/*
+ This function returns an attribute value from the current *map record's attributes
+ in the currently active layers.
+ It bypasses the GVT attribute caching scheme, using SQL calls to get the attribute data
+ from the DB which *map is connected to
+ If an attribute *name exists, it's string value will be returned. New memory will be
+ allocated for the return value and the caller is responsible for releasing it.
+
+ If the attribute does not exist, NULL will be returned. NULL will also be returned, if the
+ attribute exists, but does not have a value.
+ In addition, *ctype will be set to the attribute type or to -1, if the attribute was
+ not found. This allows you to differentiate between cases where NULL is returned because
+ of a missing attribute and where it is returned because of a missing or NULL value!
+
+*/
+char *GVT_get_attr_no_cache ( char *name, int *ctype, GVT_map_s *map ) {
+ int dbcol;
+ char *val = NULL;
+
+ /* these three are pointers to members of the attribute table
+ they need not be free'd ! */
+ dbColumn *column;
+ dbValue *dbvalue;
+ char *colname;
+
+ int sqltype, more;
+
+ if ( map->field != NULL ) {
+ *ctype = -1;
+ /* TODO: Do we really have to go through all this in order to read the current record's attribute values? */
+ /* START */
+ sprintf ( map->sql_buf, "select * from %s where %s = %d ",
+ map->field->table,
+ map->field->key,
+ map->vect_cats->cat[ map->current_layer - 1 ]);
+ db_set_string ( &map->sql, map->sql_buf );
+ db_open_select_cursor(map->driver, &map->sql, &map->cursor, DB_SEQUENTIAL); /* this does not alloc any memory, so no need to clean up */
+ map->table = db_get_cursor_table (&map->cursor);
+ db_fetch ( &map->cursor, DB_NEXT, &more );
+ /* END of stuff to do for reading attribute values */
+
+ /* now find the right attribute col and read the value in there */
+ /* TODO: need to update dbncols = db_get_table_number_of_columns (map->table) */
+ for ( dbcol = 0; dbcol < map->dbncols; dbcol ++ ) {
+ column = db_get_table_column(map->table, dbcol);
+ sqltype = db_get_column_sqltype (column);
+ *ctype = db_sqltype_to_Ctype(sqltype);
+ colname = (char*) db_get_column_name (column);
+ if (!strcmp (colname,name)){
+ dbvalue = db_get_column_value(column);
+ if ( db_test_value_isnull (dbvalue) ) {
+ /* got a null value, so return NULL! */
+ db_close_cursor(&map->cursor);
+ return ( NULL );
+ }
+ db_convert_value_to_string( dbvalue, sqltype, &map->str);
+ val = strdup ( db_get_string ( &map->str ) );
+
+ if ( GVT_DEBUG > 1 ) {
+ fprintf ( stderr, "GVT: ATT TXT = %s\n", val );
+ }
+ db_close_cursor(&map->cursor);
+ return (val);
+ }
+ }
+ db_close_cursor(&map->cursor);
+ }
+
+ /* return value will be NULL, if nothing was read or attribute missing !!!! */
+ return ( NULL );
+}
+
+/*
+ STUB FUNCTION
+*/
+char *GVT_get_attr ( char *name, int *ctype, GVT_map_s *map ) {
+ return GVT_get_attr_no_cache ( name, ctype, map );
+}
+
+
+
+/*
+ This function returns a DOUBLE attribute value from the current *map record's attributes
+ in the currently active layers.
+
+ If the attribute does not exist or is not of the right type, program execution will be aborted.
+ In all other cases, the return value will be a double value.
+*/
+double GVT_get_double ( char *name, GVT_map_s *map) {
+ double val;
+ char *att;
+ int ctype;
+
+ /* TODO: set this to standardized NULL value for DOUBLE attributes !!! */
+ val = 0.0;
+
+ /* get raw attribute text */
+ att = GVT_get_attr ( name, &ctype, map );
+
+ if ( att == NULL ) {
+ /* something is wrong */
+ if ( ctype == -1 ) {
+ G_fatal_error (_( "GVT: Missing attribute '%s' in map '%s'.\n"), name, map->name );
+ }
+ if ( ctype != GVT_DOUBLE ) {
+ G_fatal_error (_( "GVT: Attribute '%s' in map '%s' not of type DOUBLE.\n"), name, map->name);
+ }
+ /* double attribute exists, but has no value */
+ return ( val );
+ }
+
+ val = atof ( att );
+ G_free ( att );
+
+ if ( GVT_DEBUG > 1 ) {
+ fprintf ( stderr, "GVT: DBL VAL = %.f\n", val );
+ }
+
+ return (val);
+}
+
+
+/*
+ This function returns an INTEGER attribute value from the current *map record's attributes
+ in the currently active layers.
+
+ If the attribute does not exist or is not of the right type, program execution will be aborted.
+ In all other cases, the return value will be a double value.
+*/
+int GVT_get_int ( char *name, GVT_map_s *map) {
+ int val;
+ char *att;
+ int ctype;
+
+ /* TODO: set this to standardized NULL value for INTEGER attributes !!! */
+ val = 0;
+
+ /* get raw attribute text */
+ att = GVT_get_attr ( name, &ctype, map );
+
+ if ( att == NULL ) {
+ /* something is wrong */
+ if ( ctype == -1 ) {
+ G_fatal_error (_("GVT: Missing attribute '%s' in map '%s'.\n"), name, map->name );
+ }
+ if ( ctype != GVT_INT ) {
+ G_fatal_error (_("GVT: Attribute '%s' in map '%s' not of type INTEGER.\n"), name, map->name );
+ }
+ /* integer attribute exists, but has no value */
+ return ( val );
+ }
+
+ val = atoi ( att );
+ G_free ( att );
+
+ if ( GVT_DEBUG > 1 ) {
+ fprintf ( stderr, "GVT: INT VAL = %.i\n", val );
+ }
+
+ return (val);
+}
+
+
+/*
+ This will read an attribute that's either double or integer type and always
+ returns the value as a double type.
+
+ The return value is the same as for function GVT_get_double ().
+*/
+double GVT_get_numeric ( char *name, GVT_map_s *map) {
+ double val;
+ char *att;
+ int ctype;
+
+ /* TODO: set this to standardized NULL value for DOUBLE attributes !!! */
+ val = 0.0;
+
+ /* get raw attribute text */
+ att = GVT_get_attr ( name, &ctype, map );
+
+ if ( att == NULL ) {
+ /* something is wrong */
+ if ( ctype == -1 ) {
+ G_fatal_error (_("GVT: Missing attribute '%s' in map '%s'.\n"), name, map->name );
+ }
+ if ( ctype != (GVT_DOUBLE || GVT_INT) ) {
+ G_fatal_error (_("GVT: Attribute '%s' in map '%s' not of type DOUBLE or INTEGER.\n"), name, map->name );
+ }
+ /* attribute exists, but has no value */
+ return ( val );
+ }
+
+ val = atof ( att );
+ G_free ( att );
+
+ if ( GVT_DEBUG > 1 ) {
+ fprintf ( stderr, "GVT: NUM VAL = %.f\n", val );
+ }
+
+ return (val);
+}
+
+
+/*
+ This function returns a STRING attribute value from the current *map record's attributes
+ in the currently active layers.
+
+ If the attribute does not exist or is not of the right type, program execution will be aborted.
+ In all other cases, the return value will be a string value (NULL if the attribute did not have a value).
+
+ New memory will be allocated for the string returned and must be free'd by the caller.
+*/
+char *GVT_get_string ( char *name, GVT_map_s *map) {
+ char *att;
+ int ctype;
+
+
+ /* get raw attribute text */
+ att = GVT_get_attr ( name, &ctype, map );
+
+ if ( att == NULL ) {
+ /* something is wrong */
+ if ( ctype == -1 ) {
+ G_fatal_error (_( "GVT: Missing attribute '%s' in map '%s'.\n"), name, map->name );
+ }
+ if ( ctype != GVT_STRING ) {
+ G_fatal_error (_("GVT: Attribute '%s' in map '%s' not of type STRING.\n"), name, map->name );
+ }
+ /* double attribute exists, but has no value */
+ return ( NULL );
+ }
+
+ if ( GVT_DEBUG > 1 ) {
+ fprintf ( stderr, "GVT: STR VAL = %s\n", att );
+ }
+
+ return ( att );
+}
+
+
+/* TODO: implement C API attribute writes in GVT_set_double2 and delete this version ! */
+/* this writes a new value into a double type attribute of the current record */
+int GVT_set_double ( char *name, double dvalue, GVT_map_s *map ) {
+
+ char tmp [GVT_MAX_STRING_SIZE];
+
+ map->field = Vect_get_field (&map->in_vect_map, map->vect_cats->field[ map->current_layer - 1 ]);
+ /* TODO: in which cases could field == NULL? what to do then? */
+ /* TODO: catch cases without a DB connection */
+ if ( map->field != NULL ) {
+ sprintf ( tmp, "echo \"UPDATE %s SET %s=%f WHERE cat=%d\" | db.execute", map->name, name, dvalue, map->vect_cats->cat[map->current_layer - 1]);
+ system ( tmp );
+ }
+
+ return ( TRUE );
+}
+
+
+
+/* TODO: this is far too much overhead for reading a single attribute! */
+/* (same as GVT_get_double) */
+/* write a value into a double attribute */
+/* returns:
+ 1 (TRUE) SUCCESS: VALUE SET
+ 0 (FALSE) FAILURE
+*/
+int GVT_set_double2 ( char *name, double dvalue, GVT_map_s *map ) {
+ int dbcol;
+ int found;
+
+ dbColumn *column;
+ dbValue *dbvalue;
+ int ctype, sqltype, more;
+ char *colname;
+
+ map->field = Vect_get_field (&map->in_vect_map, map->vect_cats->field[ map->current_layer - 1 ]);
+ /* TODO: in which cases could field == NULL? what to do then? */
+ /* TODO: catch cases without a DB connection */
+ if ( map->field != NULL ) {
+ /* TODO: we should not need to re-open everything, if GVT_open_map () does not close it */
+ map->driver = db_start_driver ( map->field->driver );
+ db_init_handle ( &map->handle );
+ db_set_handle ( &map->handle, map->field->database, NULL);
+ if (db_open_database(map->driver, &map->handle) != DB_OK){
+ /* TODO: DB access failed */
+ db_shutdown_driver( map->driver );
+ } else {
+ sprintf ( map->sql_buf, "select * from %s where %s = %d ",
+ map->field->table,
+ map->field->key,
+ map->vect_cats->cat[ map->current_layer - 1 ]);
+ db_set_string ( &map->sql, map->sql_buf );
+ db_open_select_cursor(map->driver, &map->sql, &map->cursor, DB_SEQUENTIAL);
+
+ map->table = db_get_cursor_table (&map->cursor);
+ db_fetch ( &map->cursor, DB_NEXT, &more );
+ map->dbncols = db_get_table_number_of_columns (map->table);
+ /* TODO: lots of the vars in here should not be part of GVT_map_s struct
+ many of them need to be tracked for each layer separately
+ */
+ found = FALSE;
+ /* TODO: too much overhead: better to keep a list of attributes in GVT_map_s and search that! */
+ for ( dbcol = 0; dbcol < map->dbncols; dbcol ++ ) {
+ column = db_get_table_column(map->table, dbcol);
+ sqltype = db_get_column_sqltype (column);
+ ctype = db_sqltype_to_Ctype(sqltype);
+ colname = (char*) db_get_column_name (column);
+ if (!strcmp (colname,name)){
+ if (ctype == GVT_DOUBLE) {
+ found = TRUE;
+ dbvalue = db_get_column_value(column); /* dbvalue now points to attribute */
+ db_convert_value_to_string( dbvalue, sqltype, &map->str);
+ if ( db_test_value_isnull ( dbvalue ) ) {
+ fprintf ( stderr, " NULL\n" );
+ }
+ db_set_value_double (dbvalue, dvalue);
+ db_convert_value_to_string( dbvalue, sqltype, &map->str);
+ fprintf ( stderr, "SETTING %s TO: %f\n", db_get_string ( &map->str ), dvalue );
+
+ /* TODO: now we could actually skip the other runs of this loop! */
+ } else {
+ /* TODO: abort if not double attribute or return something strange to signify this */
+ return (FALSE);
+ }
+ }
+ }
+ db_close_cursor(&map->cursor);
+ db_close_database( map->driver);
+ db_shutdown_driver( map->driver);
+ if ( found == FALSE ) {
+ /* TODO: check if attribute exists */
+ return (FALSE);
+ }
+ }
+ }
+
+ return ( TRUE );
+}
+
+
+/* check if an attribute is NULL */
+int GVT_is_null () {
+ /* look at functions in lib/db/dbmi_base:
+ column.c: db_test_column_null_allowed ();
+ value.c: db_test_value_isnull (dbValue *value)
+ */
+ return ( FALSE );
+}
+
+
+int GVT_set_null () {
+ return ( TRUE );
+}
+
+
+/* append new record
+*/
+int GVT_append ( ) {
+ return ( TRUE );
+}
+
+
+/* delete record
+*/
+int GVT_delete ( ) {
+ return ( TRUE );
+}
+
+
+/* seek to first record that matches
+ a search string
+*/
+int GVT_find () {
+ return ( TRUE );
+}
+
+
+int GVT_find_next () {
+ return ( TRUE );
+}
+
+
+/* create a new double attribute
+ returns:
+ 1 (TRUE) SUCCESS
+ 0 (FALSE) FAILURE
+*/
+int GVT_add_double ( char *name, GVT_map_s *map ) {
+ char buf [GVT_MAX_STRING_SIZE];
+
+ /* TODO: replace system call with appropriate C API calls */
+ /* TODO: check if column with same name (but different type) exists */
+
+ sprintf ( buf, "echo \"ALTER TABLE %s ADD COLUMN %s double precision\" | db.execute", map->name, name );
+ fprintf ( stderr, "%s\n", buf );
+ system ( buf );
+
+ return ( TRUE );
+}
+
+
+/* delete an attribute from a table
+*/
+int GVT_drop_double ( char *name, GVT_map_s *map ) {
+ /* TODO: replace system call with appropriate C API calls */
+ /* NOTE: db.execute seems to not support DROP COLUMN statements! */
+ return ( TRUE );
+}
+
+
+long int GVT_get_num_records ( GVT_map_s *map ) {
+ return ( map->num_records );
+}
+
+
+int GVT_is_lat_long ( GVT_map_s *map ) {
+ return ( map->isLatLong );
+}
+
+
+/* pass NULL for constraints that you don't want to touch */
+int GVT_set_constraints ( struct Cell_head *window, char *SQL, GVT_map_s *map ) {
+
+ if ( window != NULL ) {
+ map->north = window->north;
+ map->south = window->south;
+ map->east = window->east;
+ map->west = window->west;
+ map->top = window->top;
+ map->bottom = window->bottom;
+ }
+
+ if ( SQL != NULL ) {
+ }
+ return (0);
+}
+
+
+
+int GVT_get_constraints ( void ) {
+ return (0);
+}
+
+
+int GVT_constraints ( int type, int SQL, int region ) {
+ return (0);
+}
+
+
+/* returns the file name of the map */
+char *GVT_get_mapname ( GVT_map_s *map ) {
+ /* first check, if opened */
+ return (NULL);
+}
+
+
+/* returns the mapset of the map */
+char *GVT_get_mapset ( GVT_map_s *map ) {
+ /* first check, if opened */
+ return (NULL);
+}
+
+
+/* ANALYTICAL FUNCTIONS */
+
+/* gets the 2D distance from the point in the map's current record
+ to the point defined by x and y
+ This also works for lat long locations
+*/
+double GVT_get_2D_point_distance ( double x, double y, GVT_map_s *map ) {
+ double d = 0.0;
+ double px, py;
+ double a,b,c;
+
+ /* TODO: check if it's a points map! */
+
+ if ( map->isLatLong ) {
+ px = *map->vect_points->x;
+ py = *map->vect_points->y;
+ d = G_geodesic_distance ( x, y, px, py );
+ } else {
+ px = *map->vect_points->x;
+ py = *map->vect_points->y;
+ a = fabs (px - x);
+ b = fabs (py - y);
+ c = sqrt ( pow (a,2.0) + pow (b,2.0) );
+ d = c;
+ }
+
+ return (d);
+}
+
+double GVT_get_point_x ( GVT_map_s *map ) {
+ /* TODO: check if this is a point map */
+ return ( *map->vect_points->x );
+}
+
+double GVT_get_point_y ( GVT_map_s *map ) {
+ /* TODO: check if this is a point map */
+ return ( *map->vect_points->y );
+}
+
+
+
+
Added: grass-addons/raster/r.xtent/gt_vector.h
===================================================================
--- grass-addons/raster/r.xtent/gt_vector.h (rev 0)
+++ grass-addons/raster/r.xtent/gt_vector.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,148 @@
+/***************************************************************************
+ * gt_vector.h
+ *
+ * Mon Apr 18 15:23:06 2005
+ * Copyright 2005 Benjamin Ducke
+ ****************************************************************************/
+
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ */
+
+#ifndef _GT_VECTOR_H
+#define _GT_VECTOR_H
+
+#include <grass/gis.h>
+#include <grass/dbmi.h>
+#include <grass/Vect.h>
+
+#define CLEAN 0
+#define DIRTY 1
+
+#define GVT_MAX_STRING_SIZE 4096
+#define GVT_EMPTY_SQL "EMPTY"
+
+/* GVT data types match column types in grass/dbmi.h */
+#define GVT_DOUBLE DB_C_TYPE_DOUBLE
+#define GVT_INT DB_C_TYPE_INT
+#define GVT_STRING DB_C_TYPE_STRING
+
+typedef struct GVT_map {
+
+ char name[GVT_MAX_STRING_SIZE];
+ char mapset[GVT_MAX_STRING_SIZE];
+
+ int num_layers; /* layers are called fields in the vect API */
+ int current_layer; /* -1 means: parse ALL layers! */
+ long int current_record;
+ long int num_records;
+ int num_attributes;
+
+ struct Map_info in_vect_map;
+ struct line_pnts *vect_points;
+ struct line_cats *vect_cats;
+ struct field_info *field; /* a pointer to the current layer (=field) */
+ int dbncols;
+ char sql_buf[GVT_MAX_STRING_SIZE];
+ dbString sql, str;
+ dbDriver *driver;
+ dbHandle handle;
+ dbCursor cursor;
+ dbTable *table;
+
+ /* topology properties */
+ short int open_level;
+
+ /* coordinate data properties */
+ short int is3d;
+ short int isLatLong;
+ short int type;
+
+ /* constraints */
+ double north, south, east, west, top, bottom;
+ int type_constraints;
+
+ /* attribute cache and memory ops */
+ short int *status;
+ double **cache_double;
+ short int **cache_short;
+ int **cache_int;
+ long int **cache_long;
+ char **cache_char;
+
+ size_t mem_size;
+
+} GVT_map_s;
+
+int GVT_NO_REGION_CONSTRAINTS;
+int GVT_NO_TYPE_CONSTRAINTS;
+int GVT_NO_SQL_CONSTRAINTS;
+
+int GVT_NO_ATTRIBUTE_CACHING;
+
+/* if set, GRASS module will terminate on file I/O errors */
+int GVT_FATAL_FILE_ERRORS;
+int GVT_FATAL_TYPE_ERRORS;
+int GVT_FATAL_MEM_ERRORS;
+
+int GVT_INIT_DONE;
+
+/**** FUNCTIONS *****/
+
+void GVT_init ( void );
+
+GVT_map_s *GVT_new_map ( void );
+int GVT_free_map ( GVT_map_s *map );
+
+int GVT_open_map ( char *mapname, int type, GVT_map_s *map );
+int GVT_open_map_points ( char *mapname, int type, GVT_map_s *map );
+void GVT_close_map ( GVT_map_s *map );
+
+int GVT_next ( GVT_map_s *map );
+int GVT_seek ( long int pos, GVT_map_s *map );
+void GVT_first ( GVT_map_s *map );
+void GVT_last ( GVT_map_s *map );
+void GVT_rewind ( GVT_map_s *map );
+long int GVT_get_current ( GVT_map_s *map );
+
+/* Attribute functions */
+
+/* check attributes */
+int GVT_double_exists ( char *name, GVT_map_s *map );
+int GVT_int_exists ( char *name, GVT_map_s *map );
+int GVT_numeric_exists ( char *name, GVT_map_s *map );
+int GVT_string_exists ( char *name, GVT_map_s *map );
+int GVT_any_exists ( char *name, GVT_map_s *map );
+
+/* read attributes */
+double GVT_get_double ( char *name, GVT_map_s *map);
+int GVT_get_int ( char *name, GVT_map_s *map);
+double GVT_get_numeric ( char *name, GVT_map_s *map);
+char *GVT_get_string ( char *name, GVT_map_s *map);
+
+/* modify attributes */
+int GVT_set_double ( char *name, double dvalue, GVT_map_s *map );
+int GVT_add_double ( char *name, GVT_map_s *map );
+
+/* Info functions */
+long int GVT_get_num_objects ( GVT_map_s *map );
+
+/* Geometry functions */
+double GVT_get_2D_point_distance ( double x, double y, GVT_map_s *map );
+double GVT_get_point_x ( GVT_map_s *map );
+double GVT_get_point_y ( GVT_map_s *map );
+
+
+#endif /* _GT_VECTOR_H */
Added: grass-addons/raster/r.xtent/main.c
===================================================================
--- grass-addons/raster/r.xtent/main.c (rev 0)
+++ grass-addons/raster/r.xtent/main.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,1625 @@
+/* r.xtent */
+
+
+/*
+
+ TODO Version 1.0
+
+ BUGS:
+ - Win32:
+ Complains about not being able to remove fcell element for tmapnames[i]. But there is no
+ fcell element for those maps, anyway !!! Maybe this is just a bug in current CVS
+
+ DBF TROUBLES:
+ - column names are limited to 10 chars
+
+ TEST:
+ - what happens, if some attribute values are bad in the input centers map?
+ - does this work with CENTROIDS instead of POINTS, as well ?
+ - do line breaks work for the protocol on Win32 (both regular and tabular)?
+ - NULL cells in input cost surfaces
+
+ General coding
+ - VALGRIND
+ - replace all G_system() with G_spawn(). G_spawn returns -1 on error, so this can be checked!
+ -> adjust tools.c run_cmd() to use G_spawn() instead of G_system()
+ - improve HTML documentation: math typesetting, better formatting. Pictures. Cross references.
+ - make status display work on Win32. [currently impossible?]
+ - ATTRIBUTE TYPE CHECKING: SOME COULD BE DOUBLE OR INT, BUT SHOULD ALWAYS BE HANDLED AS DOUBLE
+ -> use functions GVT_numeric_exists and GVT_get_numeric (), that always returns a double !!!
+ - Currently, only C attribute supports double and int, all other do strict type-checking. Should this
+ be relaxed?
+ - replace Posix funcs like strcpy with GRASS funcs, where available
+ - make sure that memory is dealloc'd properly for things like cats, labels, colors ...
+
+ Program functionality
+ - a switch to output areal calculations in ha
+ - let user define sqkm per cell for areal calculations in non-world systems
+
+ Documentation updates:
+ - C gets normalized automatically!
+ - territorial partitioning will always be complete, unless 'strict' flag given or reach
+ constraint set
+ - paragraph breaks between figures
+ - add more examples, with strict mode and reach constraints
+ - reach constraint together with unrestricted formula is the recommended way of defining
+ territorial reach, as it is easier to control than the original I>=0 contraint!
+
+
+ TODO VERSION 2.X
+
+ A version 1.x should be done if there is enough user demand
+
+ Political attributes:
+ - Model constraints
+ * find a way to store for each site:
+ - allies: cannot claim each other's territory: if second highest I is an ally: do not touch it!
+ - if any constraints are given, then there should be a check for validity
+ to see if they hold in the model result (?!)
+
+ Misc.:
+ - user should be able to provide a map for spatially variable k (?)
+ - in very complex situations, the model may produce fragmented territories. E.g.:
+ r.xtent cent=sites output=test comp_id=comp comp_str=strength elev=dem -g --o
+ (Oaxaca dataset, with region set to extents of 'dem' and res=0:00:59)
+ -> Should this be catered for? If so, if a center gets captured, then all of
+ its cells in the output map (basic I) should be reclassified to the ID of
+ the new "boss" (this will automatically be taken care of in cases where
+ there is an a priori boss for that center, specified in the ruler= attribute).
+ - latex, html, rtf output style for protocols
+
+
+
+
+******************************
+
+REVISON 0.99, OCT 2007
+
+******************************
+
+
+FIX BUGS:
+
+UPDATE DOCUMENTATION!
+
+geom.c: in function get_distance():
+ delete line distance measure in case where a cost map is used to derive
+ distance: only take the cost value (costs[i]) and return that! [DONE]
+
+xent.c: the debugging output for "I=..." should be set to %.10f to allow floating
+ point values in much smaller ranges [DONE]
+
+
+Improvements for the next version of r.xtent:
+
+1. Normalize center weights, so that largest center's weight=1.0
+ -> Divide all weights by largest weight. This won't do
+ any harm to data that is already normalized, as it just
+ divides everything by "1"! This is also what R and L did
+ in their original study. [DONE]
+
+2. Use standard weight of "k=0.0001" for regions that use projected
+meter systems. This is the same as R and L did (although they used km as
+their unit of distance measurements which meant "k=0.1").
+Use standard weight of "k=0.000001?" for cost surface maps generated with
+r.walk. [DONE]
+
+
+3. Provide different functions for distance decay:
+
+a. Linear (already DONE)
+b. Exponential
+c. Sigmoid-function: P(t)=(1)/(1+e^t)
+ --> see Wikipedia entry on Sigmoid
+ For our purposes : P(t)=( 1- (1)/(1+e^t) )
+
+
+*** Open questions:
+
+
+a) How to model the hierarchical relations of sites that form alliances?
+
+b) How to choose "k" for regions that use lat-long and plain x/y systems?
+
+c) what to do with "residual territories", i.e. territories still claimed by a site A that has been
+taken over by another site B? Maybe allocate all territory that once belonged to A to B in a
+second "clean-up" pass?
+
+
+
+
+
+
+*/
+
+#define LOCAL
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <time.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "globals.h"
+#include "at_exit_funcs.h"
+#include "gt_vector.h"
+#include "xtent.h"
+#include "tools.h"
+
+
+/* variables global in main.c */
+GVT_map_s *pmap;
+double b_weight; /* weight for boundaries */
+double p_weight; /* weight for pathways */
+int *costf;
+CELL *cell;
+CELL *cell_second;
+DCELL *dcell;
+DCELL *dcell_boundaries;
+DCELL **costdcell;
+double *costs;
+double *reach;
+int *ruler;
+double *max_cost;
+
+char *boundaries_map;
+char *pathways_map;
+
+/* temporary map name strings */
+char *cmapname; /* name of a cost map */
+char *pseudomap; /* a pseudo friction map (all costs set to 1), used when only boundaries or pathways are supplied */
+
+
+
+/* DEFINE FLAGS AND OPTIONS */
+void define_opts ( int argc, char *argv[] ) {
+
+ module = G_define_module ();
+ module->description = "Calculate variations or original version of Renfrew and Level's Xtent model for territorial modeling.";
+ module->label = "r.xtent: territorial modelling with boundaries and movement costs";
+ module->keywords = _("raster,Xtent,territories,territorial analysis");
+
+ /* vector map with centers (sites only or network) */
+ parm.centers = G_define_standard_option (G_OPT_V_INPUT);
+ parm.centers->key = "centers";
+ parm.centers->type = TYPE_STRING;
+ parm.centers->required = YES;
+ parm.centers->description = _("Input map with vector points representing centers");
+ parm.centers->guisection = _("Basic");
+
+ /* Output map name */
+ parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.output->key = "territories";
+ parm.output->description = _("Output raster map showing territories");
+ parm.output->guisection = _("Basic");
+
+ /* attribute that stores C for each center */
+ parm.c = G_define_option ();
+ parm.c->key = "c";
+ parm.c->type = TYPE_STRING;
+ parm.c->required = NO;
+ parm.c->description = _("Double or integer attribute in 'centers' map with individual 'C' values");
+ parm.c->guisection = _("Basic");
+
+ /* global a */
+ parm.a = G_define_option ();
+ parm.a->key = "a";
+ parm.a->type = TYPE_DOUBLE;
+ parm.a->required = YES;
+ parm.a->answer = "0.5";
+ parm.a->description = _("Global setting for 'a'");
+ parm.a->guisection = _("Basic");
+
+ /* global k */
+ parm.k = G_define_option ();
+ parm.k->key = "k";
+ parm.k->type = TYPE_DOUBLE;
+ parm.k->required = NO;
+ parm.k->description = _("Global setting for 'k'");
+ parm.k->guisection = _("Basic");
+
+ /* cost surface attribute */
+ parm.costs_att = G_define_option ();
+ parm.costs_att->key = "costs_att";
+ parm.costs_att->type = TYPE_STRING;
+ parm.costs_att->required = NO;
+ parm.costs_att->description = _("Text attribute in 'centers' map that stores names of raster cost maps");
+ parm.costs_att->guisection = _("Basic");
+
+ /* write ASCII report file? */
+ parm.report = G_define_standard_option (G_OPT_F_OUTPUT);
+ parm.report->key = "report";
+ parm.report->type = TYPE_STRING;
+ parm.report->required = NO;
+ parm.report->description = _("Name of ASCII file to write report to");
+ parm.report->guisection = _("Basic");
+
+ /* impose I >= 0 restriction? */
+ flag.strict = G_define_flag ();
+ flag.strict->key = 's';
+ flag.strict->description = _("Impose original restriction I >= 0");
+ flag.strict->guisection = _("Basic");
+
+ /* tabular format report file? */
+ flag.tabular = G_define_flag ();
+ flag.tabular->key = 't';
+ flag.tabular->description = _("Tabular output format for report file");
+ flag.tabular->guisection = _("Basic");
+
+ parm.second = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.second->key = "comp_id";
+ parm.second->type = TYPE_STRING;
+ parm.second->required = NO;
+ parm.second->description = _("New raster map to store ID of strongest competitor");
+ parm.second->guisection = _("Enhancements");
+
+ /* Error map/attribute name */
+ parm.errors = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.errors->key = "comp_strength";
+ parm.errors->type = TYPE_STRING;
+ parm.errors->required = NO;
+ parm.errors->description = _("New raster map to store strength of competitor");
+ parm.errors->guisection = _("Enhancements");
+
+ parm.maxdist = G_define_option ();
+ parm.maxdist->key = "reach";
+ parm.maxdist->type = TYPE_STRING;
+ parm.maxdist->required = NO;
+ parm.maxdist->description = _("Double attribute in 'centers' map for maximum distance/cost reach");
+ parm.maxdist->guisection = _("Enhancements");
+
+ parm.ruler = G_define_option ();
+ parm.ruler->key = "ruler";
+ parm.ruler->type = TYPE_STRING;
+ parm.ruler->required = NO;
+ parm.ruler->description = _("Integer attribute in 'centers' map that points to ID of ruling center");
+ parm.ruler->guisection = _("Enhancements");
+
+ #ifdef EXPERIMENTAL
+ parm.ally = G_define_option ();
+ parm.ally->key = "ally";
+ parm.ally->type = TYPE_STRING;
+ parm.ally->required = NO;
+ parm.ally->description = _("Text attribute in 'centers' map with comma-separated list of ally IDs");
+ parm.ally->guisection = _("Enhancements");
+ #endif
+
+ /* indexing attribute */
+ parm.cats = G_define_option ();
+ parm.cats->key = "categories";
+ parm.cats->type = TYPE_STRING;
+ parm.cats->required = NO;
+ parm.cats->answer = "cat";
+ parm.cats->description = _("Integer attribute in 'centers' map to use as output raster map categories");
+ parm.cats->guisection = _("Mapping");
+
+ /* can be used to label raster categories (mode 'areas') or add a labeling attribute in mode 'sites'*/
+ parm.labels = G_define_option ();
+ parm.labels->key = "labels";
+ parm.labels->type = TYPE_STRING;
+ parm.labels->required = NO;
+ parm.labels->description = _("Text attribute in 'centers' map for labeling output raster map categories");
+ parm.labels->guisection = _("Mapping");
+
+ /* RGB colouring attribute */
+ parm.rgbcol = G_define_option ();
+ parm.rgbcol->key = "rgb_column";
+ parm.rgbcol->type = TYPE_STRING;
+ parm.rgbcol->required = NO;
+ parm.rgbcol->answer = "GRASSRGB";
+ parm.rgbcol->description = _("String attribute (RRR:GGG:BBB) in 'centers' map for category colors");
+ parm.rgbcol->guisection = _("Mapping");
+
+ /* setup some basic GIS stuff */
+ G_gisinit (argv[0]);
+ if (G_parser(argc, argv))
+ exit(-1);
+
+ /* do not pause after a warning message was displayed */
+ G_sleep_on_error (0);
+
+}
+
+
+/* check for valid combinations of parameters/options */
+/* check for valid map and attribute names and types */
+void check_params ( void ) {
+
+ int num_specs;
+ double k_val;
+
+
+ /* enable verbosity, overwrite, progess display? */
+ VERBOSE = 0;
+ OVERWRITE = 0;
+ PROGRESS = 1;
+ if ( module->verbose ) {
+ VERBOSE = 1;
+ }
+ if ( module->overwrite ) {
+ OVERWRITE = 1;
+ }
+
+ /* stuff relating to report file writing */
+ if ( flag.tabular->answer ) {
+ if ( parm.report->answer == NULL ) {
+ G_warning (_("-t flag only makes sense if option 'report=' is specified. Ignoring.\n"));
+ }
+ }
+
+ /* 1. FORMULA PROPERTIES */
+ if ( (parm.errors->answer) && (parm.second->answer==NULL) ) {
+ G_fatal_error (_("Mapping strength of competition requires providing a name for the 'competitor ID' map.\n"));
+ }
+
+ /* set appropriate value for k, if user did not provide one! */
+ if ( parm.k->answer == NULL ) {
+
+ G_warning (_("No value for k given. Will set a default value."));
+
+ parm.k->answer = G_malloc ( sizeof (char) * 255 );
+ k_val = 0.0001;
+
+ /* lat-long system */
+ if ( G_projection () == 3 ) {
+ G_warning (_("You are working in a lat-long location. Appropriate value for 'k' depends on latitude."));
+ } else {
+ /* we assume a system based on meters */
+ G_warning (_("Assuming meter-based coordinate system."));
+ }
+
+ if ( parm.costs_att->answer != NULL ) {
+ G_warning (_("Using cost maps. Appropriate value for k depends on cost units."));
+ }
+
+ sprintf ( parm.k->answer, "%f", k_val);
+ G_warning (_("Guessing k = %f. No guarantee this will work"), k_val);
+ }
+
+
+ /* 3. DISTANCE MODE */
+
+ DISTANCE = STRAIGHT; /* default is to use a straight-line distance measure */
+ COSTMODE = NONE;
+
+ /* check if user specified more than one way of naming cost maps */
+ if ( parm.costs_att->answer != NULL ) {
+ num_specs ++;
+ COSTMODE = ATTNAME;
+ DISTANCE = COST;
+ }
+}
+
+
+/* write a report to an ASCII file on disk */
+void write_report ( int *cats, char **labels, int argc, char *argv[] ) {
+ struct report_struct *report;
+ FILE *fp;
+ int error;
+ int i,j;
+ int fd_output, fd_second, fd_errors;
+ long int total_cells;
+ long int *competitor;
+ long int max;
+
+ char tmp [5000];
+
+ time_t systime;
+ clock_t proctime;
+ unsigned long timeused;
+ unsigned int days, hours, mins, secs;
+
+ double sqdiff;
+
+ struct Cell_head window;
+ int row, col;
+ int nrows, ncols;
+ double east, north;
+
+ int boss;
+
+
+ G_get_window (&window);
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+
+ fp = fopen ( parm.report->answer, "w+" );
+ if ( fp == NULL ) {
+ G_fatal_error ("Failed to open file '%s' for report writing.\n", parm.report->answer );
+ }
+
+ report = G_malloc ( num_centers * sizeof ( struct report_struct ) );
+
+ fd_output = G_open_cell_old ( parm.output->answer, G_find_cell ( parm.output->answer, "" ) );
+
+ fd_second = -1;
+ if ( parm.second->answer ) {
+ fd_second = G_open_cell_old ( parm.second->answer, G_find_cell ( parm.second->answer, "" ) );
+ }
+ fd_errors = -1;
+ if ( parm.errors->answer ) {
+ fd_errors = G_open_cell_old ( parm.errors->answer, G_find_cell ( parm.errors->answer, "" ) );
+ }
+
+ if ( PROGRESS )
+ G_message (_( "Writing report file to '%s':\n"), parm.report->answer);
+
+ /* alloc some counter arrays */
+ competitor = G_malloc ( sizeof ( long int) * num_centers );
+
+ /* make room for storing subjects */
+ for ( i=0; i < num_centers; i ++ ) {
+ report[i].num_subjects = 0;
+ report[i].subject = G_malloc ( sizeof (int) * num_centers );
+ report[i].subject_id = G_malloc ( sizeof (int) * num_centers );
+ report[i].subject_name = G_malloc ( sizeof (char*) * num_centers );
+ }
+
+ /* init area calculations */
+ error = G_begin_cell_area_calculations ();
+
+ /* get map statistics */
+ if ( PROGRESS )
+ G_message (_(" Getting map statistics:\n" ));
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ if ( error == 0 ) {
+ /* area cannot be calculated in x/y locations */
+ report[i].area = -1.0;
+ } else {
+ report[i].area = 0.0;
+ }
+
+ report[i].id = cats[i];
+ if ( parm.labels->answer ) {
+ report[i].name = strdup ( labels[i] );
+ }
+
+ /* get basic raster stats */
+ report[i].out_count = 0;
+ total_cells = 0;
+ for (row=0; row < nrows; row ++) {
+ G_get_c_raster_row ( fd_output, cell, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_c_null_value ( &cell[col]) ) {
+ if ( cell[col] == cats[i] ) {
+ report[i].out_count ++;
+ if ( error > 0 ) {
+ /* update area stats */
+ report[i].area = report[i].area + G_area_of_cell_at_row ( row );
+ }
+ }
+ total_cells ++;
+ }
+ }
+ }
+
+ /* calculate area (km^2), percentage and maximum reach of this center's territory */
+ report[i].area = report[i].area / 1000000;
+ report[i].percentage = ( ( (double) report[i].out_count / (double) total_cells) * 100.0 );
+ report[i].max_cost = max_cost[i];
+
+ /* get error stats */
+ if ( parm.errors->answer ) {
+ report[i].err_count = 0;
+ report[i].err_sum = 0.0;
+ for (row=0; row < nrows; row ++) {
+ G_get_c_raster_row ( fd_output, cell, row );
+ G_get_d_raster_row ( fd_errors, dcell, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_c_null_value ( &cell[col]) ) {
+ if ( cell[col] == cats[i] ) {
+ report[i].err_count ++;
+ report[i].err_sum = report[i].err_sum + (double) dcell[col];
+ }
+ }
+ }
+ }
+ report[i].err_avg = ((double) report[i].err_sum / report[i].err_count);
+
+ /* now get error var and std */
+ report[i].err_var = 0;
+ report[i].err_std = 0.0;
+ sqdiff = 0.0;
+ for (row=0; row < nrows; row ++) {
+ G_get_c_raster_row ( fd_output, cell, row );
+ G_get_d_raster_row ( fd_errors, dcell, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_c_null_value ( &cell[col]) ) {
+ if ( cell[col] == cats[i] ) {
+ sqdiff = sqdiff + pow ( ((double) dcell[col] - (double) report[i].err_avg), 2.0 );
+ }
+ }
+ }
+ }
+ report[i].err_var = (1.0 / (double) report[i].err_count) * sqdiff;
+ report[i].err_std = sqrt ( report[i].err_var );
+ }
+
+ /* get information about competitors (second highest I) */
+ if ( parm.second->answer ) {
+
+ /* initialize counters */
+ for ( j=0; j < num_centers; j ++ ) {
+ competitor[j] = 0;
+ }
+ report[i].aggressor = 0;
+
+ for (row=0; row < nrows; row ++) {
+ G_get_c_raster_row ( fd_output, cell, row );
+ G_get_c_raster_row ( fd_second, cell_second, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_c_null_value ( &cell[col]) ) {
+ if ( cell[col] == cats[i] ) {
+ /* we have a cell that belongs to the current center */
+
+ /* which on is the most frequent competitor in current center's territory?
+
+ * count frequencies for all centers in array of i counters
+ * we'll find the index of the largest counter when the raster map loop is done.
+ *
+ */
+ for (j=0; j < num_centers; j++) {
+ if ( cats[j] == cell_second[col] ) {
+ competitor [j] ++;
+ }
+ }
+ } else {
+ /* we have a cell that is outside the current center's territory */
+
+ /* how frequent is the current center second highest I in the territory of another? */
+ if ( cell_second[col] == cats[i] ) {
+ /* just sum up and also calculate percentage (against number of cells in foreing territories) */
+ report[i].aggressor ++;
+ }
+ }
+ }
+ }
+ }
+
+ /* find most frequent competitor */
+ max = -1;
+ report[i].competitor = 0;
+ report[i].competitor_id = cats[0];
+ for ( j=1; j< num_centers; j ++ ) {
+ if ( (competitor[j] > 0) && ( competitor[j] > max ) ) {
+ max = competitor[j];
+ report[i].competitor = j;
+ report[i].competitor_id = cats[j];
+ }
+ }
+
+ /* no competitor found! */
+ if ( max == -1 ) {
+ report[i].competitor = -1;
+ } else {
+ /* get percentage */
+ report[i].competitor_p = ((double) max / (double) report[i].out_count) * 100.0;
+ }
+ report[i].aggressor_p = ( (double)report[i].aggressor / ((double)total_cells - (double)report[i].out_count) ) * 100.0;
+ }
+
+ if ( PROGRESS ) {
+ G_percent (i, num_centers-1, 1); }
+ i ++;
+ }
+
+ /* now check for boss and/or subjects! */
+
+ /* check the cell under each center. If it's not owned by the center, then that
+ * center is dominated by another and we need to store its ID.
+ */
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ report[i].boss = -1;
+ east = G_easting_to_col ( GVT_get_point_x ( pmap ), &window );
+ north = G_northing_to_row ( GVT_get_point_y ( pmap ), &window );
+ G_get_c_raster_row ( fd_output, cell, (int) north );
+ boss = (CELL) cell[(int) east];
+ if ( boss != cats[i] ) {
+ /* this center is dominated by another center! Store that center's ID */
+ for ( j=0; j < num_centers; j ++ ) {
+ if ( cats[j] == boss ) {
+ report[i].boss = j;
+ report[i].boss_id = boss;
+ if ( parm.labels->answer ) {
+ report[i].boss_name = strdup ( labels[j] );
+ }
+ /* also: append ID and name of this center to the bosse's list */
+ report[j].subject[report[j].num_subjects] = i;
+ report[j].subject_id[report[j].num_subjects] = cats[i];
+ if ( parm.labels->answer ) {
+ report[j].subject_name[report[j].num_subjects] = strdup ( labels[i] );
+ }
+ report[j].num_subjects ++;
+ }
+ }
+ }
+ i++;
+ }
+
+ /* write records to ASCII file */
+ if ( PROGRESS )
+ G_message (_( " Writing individual records:\n" ));
+ GVT_rewind ( pmap );
+ if ( flag.tabular->answer ) {
+ /* header line for tabular report output */
+ fprintf ( fp, "id\tname\trulerid\trname\tcells\tcellsp\tarea\treach" );
+ if ( parm.second->answer ) {
+ fprintf ( fp, "\taggressor\taggperc\tcompid\tcompname\tcompperc" );
+ }
+ if ( parm.errors->answer ) {
+ fprintf ( fp, "\tstravg\tstrvar\tstrdev" );
+ }
+ fprintf ( fp, "\n" );
+ }
+
+ if ( !flag.tabular->answer ) {
+ /* make a nice header */
+ fprintf ( fp, "This is %s version %.2f\n", argv[0], PROGVERSION);
+ systime = time (NULL);
+ fprintf ( fp, "Calculation started on %s", ctime ( &systime ) );
+ fprintf ( fp, "\tlocation = %s\n", G_location() );
+ fprintf ( fp, "\tmapset = %s\n", G_mapset() );
+ sprintf ( tmp, "%s ", argv[0]);
+ for ( j=1; j < argc; j ++ ) {
+ strcat ( tmp, argv[j] );
+ if ( j < (argc - 1) ) {
+ strcat ( tmp, " " );
+ }
+ }
+ fprintf ( fp, "\tcommand = %s\n\n", tmp );
+ fprintf ( fp, "===\n\n" );
+ }
+
+ for ( i=0; i < num_centers; i++ ) {
+ if ( flag.tabular->answer ) {
+ /* tabular output: one line per center and always with complete fields */
+ fprintf ( fp, "%i", report[i].id );
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "\t%s", report[i].name );
+ } else {
+ fprintf ( fp, "\tn.a." );
+ }
+ if ( report[i].boss > -1 ) {
+ fprintf ( fp, "\t%i", report[i].boss_id );
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "\t%s", report[i].boss_name );
+ } else {
+ fprintf ( fp, "\tn.a." );
+ }
+ fprintf ( fp, "\t0" );
+ fprintf ( fp, "\t0.0" );
+ fprintf ( fp, "\t0.0" );
+ fprintf ( fp, "\t0.0" );
+ } else {
+ fprintf ( fp, "\t%i", report[i].id );
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "\t%s", report[i].name );
+ } else {
+ fprintf ( fp, "\tn.a." );
+ }
+ fprintf ( fp, "\t%li", report[i].out_count );
+ fprintf ( fp, "\t%.3f", report[i].percentage );
+ fprintf ( fp, "\t%.3f", report[i].area );
+ fprintf ( fp, "\t%.3f", report[i].max_cost );
+ }
+ if ( parm.second->answer ) {
+ if (( report[i].aggressor > 0 ) && ( report[i].boss < 0 )) {
+ fprintf ( fp, "\t%li", report[i].aggressor);
+ fprintf ( fp, "\t%.3f", report[i].aggressor_p );
+ } else {
+ fprintf ( fp, "\t0" );
+ fprintf ( fp, "\t0.0" );
+ }
+ if (( report[i].competitor > -1 ) && ( report[report[i].competitor].boss != i ) && ( report[i].boss < 0 )) {
+ fprintf ( fp, "\t%i", report[i].competitor_id);
+ fprintf ( fp, "\t%s", labels[report[i].competitor]);
+ fprintf ( fp, "\t%.3f", report[i].competitor_p );
+ } else {
+ fprintf ( fp, "\t0" );
+ fprintf ( fp, "\tn.a." );
+ fprintf ( fp, "\t0.0" );
+ }
+ }
+ if ( parm.errors->answer ) {
+ if (( report[i].competitor > -1 ) && ( report[report[i].competitor].boss != i ) && ( report[i].boss < 0 )) {
+ fprintf ( fp, "\t%.3f", report[i].err_avg );
+ fprintf ( fp, "\t%.3f", report[i].err_var );
+ fprintf ( fp, "\t%.3f", report[i].err_std );
+ } else {
+ fprintf ( fp, "\t0.0" );
+ fprintf ( fp, "\t0.0" );
+ fprintf ( fp, "\t0.0" );
+ }
+ }
+ fprintf ( fp, "\n" ) ;
+ } else {
+ /* this is regular report format: one paragraph per center */
+
+ /* write records to report file */
+ fprintf ( fp, "(center %i of %i)\n", i + 1, num_centers );
+ fprintf ( fp, "id:\t\t%i\n", report[i].id );
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "name:\t\t%s\n", report[i].name );
+ } else {
+ fprintf ( fp, "name:\t\tn.a.\n" );
+ }
+ if ( report[i].num_subjects > 0 ) {
+ fprintf ( fp, "dominates:\t" );
+ for ( j=0; j<report[i].num_subjects; j ++ ) {
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "%i (%s)", report[i].subject_id[j], report[i].subject_name[j] );
+ } else {
+ fprintf ( fp, "%i", report[i].subject_id[j] );
+ }
+ if ( j < (report[i].num_subjects-1) ) {
+ fprintf ( fp, ", " );
+ }
+ }
+ fprintf ( fp, "\n");
+ } else {
+ fprintf ( fp, "dominates:\tn.a.\n" );
+ }
+ if ( report[i].boss > -1 ) {
+ if ( parm.labels->answer ) {
+ fprintf ( fp, "dominated by:\t%i, %s\n", report[i].boss_id, report[i].boss_name );
+ } else {
+ fprintf ( fp, "dominated by:\t%i\n", report[i].boss_id );
+ }
+ } else {
+ fprintf ( fp, "dominated by:\tn.a.\n" );
+ fprintf ( fp, "cells:\t\t%li (%.3f%%)\n", report[i].out_count, report[i].percentage );
+ if ( report[i].area > -1.0 ) {
+ fprintf ( fp, "area:\t\t%.3f sqkm\n", report[i].area );
+ }
+ if ( DISTANCE == COST ) {
+ fprintf ( fp, "reach cost:\t%.3f cells\n", report[i].max_cost );
+ } else {
+ fprintf ( fp, "reach dist:\t%.3f km\n", report[i].max_cost / 1000 );
+ }
+ }
+ if ( parm.second->answer ) {
+ if (( report[i].aggressor > 0 ) && ( report[i].boss < 0 )) {
+ fprintf ( fp,"aggressor:\t%li (%.3f%%)\n", report[i].aggressor, report[i].aggressor_p );
+ } else {
+ fprintf ( fp,"aggressor:\tn.a.\n" );
+ }
+ if (( report[i].competitor > -1 ) && ( report[report[i].competitor].boss != i ) && ( report[i].boss < 0 )) {
+ fprintf ( fp,"competitor:\t%i, %s (%.3f%%)\n", report[i].competitor_id, labels[report[i].competitor], report[i].competitor_p );
+ } else {
+ fprintf ( fp,"competitor:\tn.a.\n" );
+ }
+ }
+ if ( parm.errors->answer ) {
+ if (( report[i].competitor > -1 ) && ( report[report[i].competitor].boss != i ) && ( report[i].boss < 0 )) {
+ fprintf ( fp, "strength avg:\t%.3f\n", report[i].err_avg );
+ fprintf ( fp, "strength var:\t%.3f\n", report[i].err_var );
+ fprintf ( fp, "strength std:\t%.3f\n", report[i].err_std );
+ } else {
+ fprintf ( fp, "strength avg:\tn.a.\n" );
+ fprintf ( fp, "strength var:\tn.a.\n" );
+ fprintf ( fp, "strength std:\tn.a.\n" );
+ }
+ }
+ fprintf ( fp, "\n\n" );
+ if ( PROGRESS ) {
+ G_percent ( i, num_centers-1, 1 );
+ }
+ }
+ }
+
+ if ( !flag.tabular->answer ) {
+ /* write a nice footer */
+ /* write processing time to logfile */
+ proctime = clock ();
+ timeused = (unsigned long) proctime / CLOCKS_PER_SEC;
+ days = timeused / 86400;
+ hours = (timeused - (days * 86400)) / 3600;
+ mins = (timeused - (days * 86400) - (hours * 3600)) / 60;
+ secs = (timeused - (days * 86400) - (hours * 3600) - (mins * 60));
+ systime = time (NULL);
+ fprintf ( fp, "===" );
+ fprintf ( fp, "\n\nCalculation finished %s",ctime(&systime));
+ fprintf ( fp, "Processing time: %id, %ih, %im, %is\n",
+ days, hours, mins, secs );
+ fflush ( fp );
+ }
+
+ fclose ( fp );
+
+ /* TODO: also release invidual members of the report struct */
+ G_free ( report );
+
+ /* close output map(s) */
+ G_close_cell ( fd_output );
+ if ( parm.second->answer ) {
+ G_close_cell ( fd_second );
+ }
+ if ( parm.errors->answer ) {
+ G_close_cell ( fd_errors );
+ }
+
+}
+
+/* ================================================================================
+
+ MAIN
+
+===================================================================================*/
+int main (int argc, char *argv[]) {
+
+ int i,j;
+ char tmp[5000];
+
+ struct Cell_head window;
+ int row, col;
+ int nrows, ncols;
+ double east, north;
+
+ int fd, fd_temp, fd_error;
+ int fd_second;
+ int fd_boundaries;
+
+ double x, y;
+ double max_d;
+ double max_c;
+ double C_total;
+
+ /* for raster map categories, labels and colours */
+ int *cats;
+ char **labels;
+ struct Categories pcats;
+ char **rgb;
+ int r,g,b;
+
+ /* map history */
+ struct History hist;
+
+ /* color map rules */
+ struct Colors *colors;
+ DCELL *v1, *v2;
+
+ int error;
+ double diff;
+ double max_diff;
+ double dist;
+
+ int valid_I; /* this var is set to "1", if at least one cell was found to belong to the territory of any
+ any center in the input map. Its purpose is to guard against cases, where strict mode
+ is used (-s) and the model parameters are set badly, so that the user is presented
+ with an all NULL output map
+ */
+ int matched;
+ int boss;
+
+
+ define_opts ( argc, argv );
+
+ G_get_window (&window);
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+
+ check_params ();
+
+ GVT_init ();
+ pmap = GVT_new_map ( );
+ GVT_open_map ( parm.centers->answer, GV_POINT , pmap );
+
+ a = atof ( parm.a->answer );
+ k = atof ( parm.k->answer );
+
+ C = G_malloc ( sizeof (double) * pmap->num_records );
+
+ /* init file descriptors */
+ fd_temp = -1;
+ fd_second = -1;
+ fd_boundaries = -1;
+
+ if ( parm.c->answer != NULL ) {
+ /* check if C attribute exists */
+ if ( GVT_any_exists ( parm.c->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading C attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.c->answer, parm.centers->answer );
+ }
+ if ( GVT_numeric_exists ( parm.c->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading C attribute: '%s' is not of type double or integer in input vector map '%s'.\n"),
+ parm.c->answer, parm.centers->answer );
+ }
+ i = 0;
+ GVT_rewind ( pmap );
+ C_total = 0;
+ while ( GVT_next ( pmap ) ) {
+ C[i] = GVT_get_numeric ( parm.c->answer, pmap );
+ if ( DEBUG > 1 ) {
+ G_message (_( " %i: C=%.f\n"), i, C[i]);
+ }
+ C_total = C_total + C[i];
+ i ++;
+ }
+ num_centers = i;
+
+ /* normalize C attributes */
+ C_norm = G_malloc ( sizeof (double) * num_centers );
+ for ( i = 0; i < num_centers; i ++ ) {
+ C_norm [i] = (double) C[i] / (double) C_total;
+ }
+
+ } else {
+ /* we still need to know at least the number of centers (sites) in the input map! */
+ num_centers = 0;
+ GVT_rewind ( pmap );
+ while ( GVT_next ( pmap ) ) {
+ num_centers ++;
+ }
+ }
+
+
+ /* if no attribute for C given: set it to 1.0 for each center */
+ if (parm.c->answer == NULL ) {
+ G_warning (_("No C attribute specified. Assuming constant C=1.0.\n"));
+ for ( i=0; i<pmap->num_records; i++) {
+ C[i] = 1.0;
+ }
+ }
+
+ cats = NULL;
+ /* no categories attribute specified? create a new index starting at 1! */
+ if ( parm.cats->answer == NULL ) {
+ cats = G_malloc ( sizeof (int) * num_centers );
+ for ( i = 0; i < num_centers; i ++ ) {
+ cats[i] = i + 1;
+ }
+ }
+
+ /* check if cats attribute exists and is of type integer */
+ if ( parm.cats->answer != NULL ) {
+ if ( GVT_any_exists ( parm.cats->answer, pmap) == FALSE ) {
+ G_warning (_("Error reading categories attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.cats->answer, parm.centers->answer );
+ cats = G_malloc ( sizeof (int) * num_centers );
+ for ( i = 0; i < num_centers; i ++ ) {
+ cats[i] = i + 1;
+ }
+ } else {
+ if ( GVT_int_exists ( parm.cats->answer, pmap) == FALSE ) {
+ G_warning (_("Error reading categories attribute: '%s' is not of integer type in input vector map '%s'.\n"),
+ parm.cats->answer, parm.centers->answer );
+ cats = G_malloc ( sizeof (int) * num_centers );
+ for ( i = 0; i < num_centers; i ++ ) {
+ cats[i] = i + 1;
+ }
+ }
+ }
+ if ( GVT_int_exists ( parm.cats->answer, pmap) == TRUE ) {
+ /* categories label is cool, so let's read in all category values in input vector map */
+ cats = G_malloc ( sizeof (int) * num_centers );
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ cats[i] = GVT_get_int ( parm.cats->answer, pmap);
+ i ++;
+ }
+ }
+ }
+
+ labels = NULL;
+ /* labeling attribute? */
+ if ( parm.labels->answer != NULL ) {
+ if ( GVT_any_exists ( parm.labels->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading label attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.labels->answer, parm.centers->answer );
+ }
+ if ( GVT_string_exists ( parm.labels->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading label attribute: '%s' is not of text type in input vector map '%s'.\n"),
+ parm.labels->answer, parm.centers->answer );
+ }
+ }
+
+ colors = NULL;
+ /* colouring attribute? */
+ make_colors = 0;
+ if ( parm.rgbcol->answer ) {
+ if ( GVT_any_exists ( parm.rgbcol->answer, pmap) == FALSE ) {
+ if ( strcmp ( parm.rgbcol->answer, "GRASSRGB" ) ) {
+ G_warning (_("Colors attribute: '%s' does not exist in input vector map '%s'. Using default colors.\n"),
+ parm.rgbcol->answer, parm.centers->answer );
+ }
+ } else {
+ if ( GVT_string_exists ( parm.rgbcol->answer, pmap) == FALSE ) {
+ G_warning (_("Colors attribute: '%s' is not of text type in input vector map '%s'. Using default colors.\n"),
+ parm.rgbcol->answer, parm.centers->answer );
+ }
+ }
+
+ if ( GVT_string_exists ( parm.rgbcol->answer, pmap) == TRUE ) {
+ /* colouring attribute is cool, so read in colours */
+ rgb = G_malloc ( sizeof (char*) * num_centers );
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ rgb[i] = GVT_get_string ( parm.rgbcol->answer, pmap);
+ i ++;
+ }
+ make_colors = 1;
+ }
+ }
+
+ /* allocate an array of center reaches (max cost/dist from a center to claim a cell) */
+ reach = G_malloc ( sizeof (double) * pmap->num_records );
+ for ( i=0; i < num_centers; i ++ ) {
+ reach[i] = -1.0; /* set all centers to unconstrained reach by default */
+ }
+
+ /* maximum cost/distance attribute specified ? */
+ if ( parm.maxdist->answer ) {
+ /* check for valid attribute */
+ if ( GVT_any_exists ( parm.maxdist->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading max reach attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.maxdist->answer, parm.centers->answer );
+ }
+ if ( GVT_double_exists ( parm.maxdist->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading max reach attribute: '%s' exists in input vector map '%s' but is not of type double.\n"),
+ parm.maxdist->answer, parm.centers->answer );
+ }
+ /* read in all reach values */
+ i = 0;
+ GVT_rewind ( pmap );
+ while ( GVT_next ( pmap ) ) {
+ reach[i] = GVT_get_double ( parm.maxdist->answer, pmap );
+ i ++;
+ }
+ }
+
+ /* allocate an array of ruler pointers */
+ ruler = G_malloc ( sizeof (int) * pmap->num_records );
+ for ( i=0; i < num_centers; i ++ ) {
+ ruler[i] = i; /* set all centers to be their own rulers by default */
+ }
+
+ /* ruler attribute specified ? */
+ if ( parm.ruler->answer ) {
+ /* check for valid attribute */
+ if ( GVT_any_exists ( parm.ruler->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading ruler attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.ruler->answer, parm.centers->answer );
+ }
+ if ( GVT_int_exists ( parm.ruler->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading ruler attribute: '%s' exists in input vector map '%s' but is not of type integer.\n"),
+ parm.ruler->answer, parm.centers->answer );
+ }
+ /* read in all reach values */
+ i = 0;
+ GVT_rewind ( pmap );
+ while ( GVT_next ( pmap ) ) {
+ ruler[i] = GVT_get_int ( parm.ruler->answer, pmap );
+ i ++;
+ }
+ /* now translate the category IDs to the sequential IDs in pmap */
+ for ( i=0; i <= num_centers; i ++ ) {
+ matched = 0;
+ for ( j=0; j <= num_centers; j ++ ) {
+ if ( cats[j] == ruler[i] ) {
+ ruler[i] = j;
+ matched = 1;
+ }
+ }
+ if ( matched == 0 ) {
+ G_fatal_error (_("Ruler ID '%i' in record no. %i in input map '%s' points to a non-existing record.\n"),
+ ruler[i], i, parm.output->answer );
+ }
+ }
+ }
+
+ /* register an AT EXIT function that will make sure temp maps
+ * get deleted
+ */
+ atexit ( clean_tmpfiles );
+
+ /* cost surface mode? check for valid maps or create them ad hoc !!! */
+ if ( DISTANCE == COST ) {
+
+ make_pseudo = 0;
+
+ if ( COSTMODE == ATTNAME ) { /* User has specified an attribute name for cost surfaces */
+ /* check if map name attribute exists */
+ if ( GVT_any_exists ( parm.costs_att->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading map name attribute: '%s' does not exist in input vector map '%s'.\n"),
+ parm.costs_att->answer, parm.centers->answer );
+ }
+ if ( GVT_string_exists ( parm.costs_att->answer, pmap) == FALSE ) {
+ G_fatal_error (_("Error reading map name attribute: '%s' is not of string type in input vector map '%s'.\n"),
+ parm.costs_att->answer, parm.centers->answer );
+ }
+ /* allocate mem for pointers to cost maps and open all maps */
+ costf = G_calloc ( (unsigned) num_centers, sizeof (int));
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ cmapname = GVT_get_string ( parm.costs_att->answer, pmap );
+ costf[i] = G_open_cell_old ( cmapname, G_find_cell ( cmapname, "" ) );
+ if ( costf[i] < 0 ) {
+ G_fatal_error (_("Error opening cost surface raster '%s').\n"),cmapname);
+ }
+ i ++;
+ G_free ( cmapname );
+ }
+ /* allocate mem for raster map buffers to read data from cost maps */
+ costdcell = G_calloc ( (unsigned) num_centers, sizeof (DCELL*) );
+ for ( i = 0; i < num_centers; i++ ) {
+ costdcell[i] = G_allocate_d_raster_buf();
+ }
+ /* allocate an array of doubles for storing cost values */
+ costs = G_calloc ( (unsigned) num_centers, sizeof (double));
+ }
+ }
+
+ /*
+ * CALCULATE XTENT MODEL
+ */
+
+ if ( PROGRESS )
+ G_message (_("Calculating Xtent model for %i centers:\n"), num_centers);
+
+ valid_I = 0;
+
+ /* this globally tracks the maximum distance to each center from within
+ * any cell in its territory. We may need this for report writing later
+ */
+ max_cost = G_malloc ( num_centers * sizeof ( double ) );
+ for ( i = 0; i < num_centers; i ++ ) {
+ max_cost[i] = -1.0;
+ }
+
+ cell = G_allocate_c_raster_buf();
+ fd = G_open_cell_new ( parm.output->answer );
+ if ( fd < 0 ) {
+ G_fatal_error (_("Could not open output map %s for writing.\n"), parm.output->answer);
+ }
+
+ /* also output a map of "second best" centers?*/
+ if ( parm.second->answer != NULL ) {
+ cell_second = G_allocate_c_raster_buf();
+ fd_second = G_open_cell_new ( parm.second->answer );
+ if ( fd_second < 0 ) {
+ G_fatal_error (_("Could not open 'second highest I' map %s for writing.\n"), parm.second->answer);
+ }
+ }
+
+ /* also output a map of error probability ? */
+ if ( parm.errors->answer != NULL ) {
+ /* for now, we'll just create a temporary map to store the
+ * unnormalized error
+ */
+ temapname = mktmpmap ( "XTENTtmp" );
+ dcell = G_allocate_d_raster_buf();
+ G_set_fp_type(DCELL_TYPE);
+ fd_temp = G_open_fp_cell_new ( temapname );
+ if ( fd_temp < 0 ) {
+ G_fatal_error (_("Could not open temporary error map %s for writing.\n"), temapname );
+ }
+ make_error = 1;
+ }
+
+ max_c = 1;
+ max_d = 1;
+
+ /* write a raster map with an I value for each cell */
+ if ( PROGRESS )
+ G_message (_(" Calculating I values for all cells in region:\n"));
+ for (row=0; row < nrows; row ++) {
+ G_set_c_null_value ( cell, ncols );
+ if ( parm.errors->answer != NULL ) {
+ G_set_d_null_value ( dcell, ncols );
+ }
+ if ( DISTANCE == COST ) {
+ /* if we have cost surfaces, we need to read one row from each map */
+ for ( i = 0; i < num_centers; i++ ) {
+ G_get_d_raster_row ( costf[i], costdcell[i], row );
+ }
+ }
+ y = G_row_to_northing ( (double) row, &window );
+ for (col=0; col < ncols; col ++) {
+ if ( DISTANCE == COST ) {
+ /* if we have cost surfaces, we need to read the cost value for each center in this column */
+ for ( i = 0; i < num_centers; i++ ) {
+ costs[i] = (double) costdcell[i][col];
+ }
+ }
+ x = G_col_to_easting ( (double) col, &window );
+
+ /* calculate XTENT value I at current coordinates */
+ if ( parm.c->answer != NULL ) {
+ I = get_center_xtent_original ( x, y, C_norm, a, k, costs, reach, ruler, pmap );
+ } else {
+ I = get_center_xtent_original ( x, y, C, a, k, costs, reach, ruler, pmap );
+ }
+
+ /* update max cost value for this center if necessary */
+ if ( I >= 0 ) {
+ if ( DISTANCE == COST ) {
+ if ( max_cost[I] < costs[I] ) {
+ if ( !G_is_d_null_value ( (DCELL*) &costs[I] ) ) {
+ max_cost[I] = costs[I];
+ }
+ }
+ }
+ if ( DISTANCE == STRAIGHT ) {
+ dist = GVT_get_2D_point_distance ( x, y, pmap );
+ if ( dist > max_cost[I] ) {
+ max_cost[I] = dist;
+ }
+ }
+ }
+
+ if ( I < 0 ) {
+ /* this can only happen if there is an error calculating I or strict mode is turned on */
+ /* otherwise, the return value of I will be the INDEX of the dominant center ! */
+ G_set_c_null_value ( &cell[col], 1 );
+ } else {
+ cell[col] = (CELL) cats[I];
+ }
+
+ /* this signals that the model has managed to produce at least one center index ;) */
+ if ( (!G_is_c_null_value ( &cell[col] )) && ( I >= 0 ) ) {
+ valid_I = 1;
+ }
+
+ /* if one of the cost surfaces has a NULL value here, we will copy it thru, as well! */
+ if ( DISTANCE == COST ) {
+ for ( i = 0; i < num_centers; i++ ) {
+ if ( G_is_d_null_value (&costdcell[i][col])) {
+ G_set_c_null_value ( &cell[col], 1 );
+ }
+ }
+ }
+ }
+
+ G_put_map_row ( fd, cell );
+
+ if ( PROGRESS ) {
+ G_percent (row, (nrows-1), 2);
+ }
+ }
+
+ /* close basic output map for now */
+ G_close_cell ( fd );
+
+ /* Now we make a second pass to find the strongest competitor for each cell.
+ * This code is full of redundant terminology. The following are all
+ * basically the same: second I, second biggest I, competitor
+ *
+ * The competitor is the center which would have claimed a cell if another,
+ * more dominant center had not taken it. The difference between the influences
+ * of the dominant center and the competitor is stored in the 'errors' map.
+ * It can be interpreted both as the strength of competition and the likelihood
+ * of making an error when claiming that a specific cell is in the territory of
+ * center A rather then B.
+ *
+ * The error (= competition strength) measure is normalized in another pass after this loop.
+ *
+ */
+ max_diff = 0;
+ if ( parm.second->answer ) {
+
+ /* re-open basic output map */
+ fd = G_open_cell_old ( parm.output->answer, G_find_cell ( parm.output->answer, "" ) );
+
+ if ( PROGRESS )
+ G_message ( " Finding competitors for all cells in region:\n");
+
+ /* If center A is dominated by center B, then it can no longer be a competitor of center B !
+ * So we need to first identify all such situations and update the 'ruler' array accordingly.
+ */
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ east = G_easting_to_col ( GVT_get_point_x ( pmap ), &window );
+ north = G_northing_to_row ( GVT_get_point_y ( pmap ), &window );
+ G_get_c_raster_row ( fd, cell, (int) north );
+ boss = (CELL) cell[(int) east];
+ if ( boss != cats[i] ) {
+ /* this center is dominated by another center! Store that center's ID */
+ for ( j=0; j < num_centers; j ++ ) {
+ if ( cats[j] == boss ) {
+ ruler[i] = j;
+ }
+ }
+ }
+ i++;
+ }
+
+ for (row=0; row < nrows; row ++) {
+
+ /* initialize error map to be all NULL cells */
+ if ( parm.errors->answer != NULL ) {
+ G_set_d_null_value ( dcell, ncols );
+ }
+
+ if ( DISTANCE == COST ) {
+ /* if we have cost surfaces, we need to read one row from each map */
+ for ( i = 0; i < num_centers; i++ ) {
+ G_get_d_raster_row ( costf[i], costdcell[i], row );
+ }
+ }
+
+ y = G_row_to_northing ( (double) row, &window );
+
+ for (col=0; col < ncols; col ++) {
+ if ( DISTANCE == COST ) {
+ /* if we have cost surfaces, we need to read the cost value for each center in this column */
+ for ( i = 0; i < num_centers; i++ ) {
+ costs[i] = (double) costdcell[i][col];
+ }
+ }
+ x = G_col_to_easting ( (double) col, &window );
+
+ /* get second highest I */
+ if ( parm.second->answer ) {
+ if ( parm.c->answer != NULL ) {
+ I_second = get_center_xtent_second ( x, y, C_norm, a, k, &diff, costs, ruler, pmap );
+ } else {
+ I_second = get_center_xtent_second ( x, y, C, a, k, &diff, costs, ruler, pmap );
+ }
+ if ( diff > max_diff ) {
+ max_diff = diff;
+ }
+ }
+
+ if ( G_is_c_null_value ( &cell[col] ) ) {
+ /* if I map has a NULL val here, we need to copy it through */
+ /* because by definition: if I is invalid, then so is 'second I' (competitor) */
+ G_set_c_null_value ( &cell_second[col], 1 );
+ if ( parm.errors->answer != NULL ) {
+ G_set_d_null_value ( &dcell[col], 1 );
+ }
+ } else {
+ /* if not, we set the current cell to the second biggest I (=ID of competitor) */
+ cell_second[col] = (CELL) cats[I_second];
+ if ( parm.errors->answer != NULL ) {
+ /* at this point: just store the unnormalized strength of competition
+ * We will make a second pass later of the 'errors' map to normalize this measure
+ */
+ dcell[col] = diff;
+ }
+ }
+
+ /* a competitor and competitor's strength (error) can only be calculated where at least
+ * two territories overlap!
+ */
+ if ( I_second < 0 ) {
+ /* this can only happen if there is an error calculating second I or strict mode is turned on */
+ /* otherwise, the return value of I will be the ID of the second most dominant center
+ * as set above */
+ G_set_c_null_value ( &cell_second[col], 1 );
+ if ( parm.errors->answer != NULL ) {
+ G_set_d_null_value ( &dcell[col], 1 );
+ }
+ }
+
+ /* if one of the cost surfaces has a NULL value here, we will copy it thru, as well! */
+ if ( DISTANCE == COST ) {
+ for ( i = 0; i < num_centers; i++ ) {
+ if ( G_is_d_null_value (&costdcell[i][col])) {
+ G_set_c_null_value ( &cell_second[col], 1 );
+ if ( parm.errors->answer != NULL ) {
+ /* also: copy thru to error map! */
+ G_set_d_null_value ( &dcell[col], 1 );
+ }
+ }
+ }
+ }
+ }
+
+ G_put_map_row ( fd_second, cell_second );
+
+ if ( parm.errors->answer != NULL ) {
+ G_put_d_raster_row ( fd_temp, dcell );
+ }
+
+ if ( PROGRESS ) {
+ G_percent (row, (nrows-1), 2);
+ }
+ }
+
+ G_close_cell ( fd );
+ G_close_cell ( fd_second );
+ }
+
+ /* Normalize values in error output map ! */
+ if ( parm.errors->answer != NULL ) {
+ /* we need to make a second pass and normalize the error
+ * measuere over the whole error map
+ */
+
+ /* close and re-open the (temporary ) error map */
+ G_close_cell ( fd_temp );
+ fd_temp = G_open_cell_old ( temapname, G_find_cell ( temapname, "" ) );
+ if ( fd_temp < 0 ) {
+ G_fatal_error (_("Could not open temporary error map %s for reading.\n"), temapname );
+ }
+
+ /* create the new final (normalized) error map for writing */
+ G_set_fp_type(DCELL_TYPE);
+ fd_error = G_open_fp_cell_new ( parm.errors->answer );
+ if ( fd_error < 0 ) {
+ G_fatal_error (_("Could not open error map %s for writing.\n"), parm.errors->answer );
+ }
+
+ if ( max_diff == 0 ) { /* protect against extreme cases and divisio by zero */
+ G_warning (_("The maximum difference between I and second highest I in your map is 0.\n"));
+ if ( PROGRESS )
+ G_message (_( " Setting all error measures to 1.0:\n"));
+ for (row=0; row < nrows; row ++) {
+ G_get_d_raster_row ( fd_temp, dcell, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_d_null_value ( &dcell[col]) ) {
+ dcell[col] = 1.0;
+ }
+ }
+ G_put_d_raster_row ( fd_error, dcell );
+ if ( PROGRESS ) {
+ G_percent (row, (nrows-1), 2); }
+ }
+ } else {
+ if ( PROGRESS )
+ G_message ( " Normalizing error measures:\n");
+ for (row=0; row < nrows; row ++) {
+ G_get_d_raster_row ( fd_temp, dcell, row );
+ for (col=0; col < ncols; col ++) {
+ if ( !G_is_d_null_value ( &dcell[col]) ) {
+ /* this is the normalized error measure:
+ * the closer to "1", the higher the risk of error
+ */
+ dcell[col] = 1 - dcell[col] / max_diff;
+ }
+ }
+ G_put_d_raster_row ( fd_error, dcell );
+ if ( PROGRESS ) {
+ G_percent (row, (nrows-1), 2); }
+ }
+ }
+
+ /* close raster maps */
+ G_close_cell ( fd_error );
+ G_close_cell ( fd_temp );
+
+ /* make a nice gray-shade color table for output map */
+ colors = G_malloc ( sizeof (struct Colors));
+ G_init_colors ( colors );
+ v1 = (DCELL*) G_malloc (sizeof (DCELL));
+ v2 = (DCELL*) G_malloc (sizeof (DCELL));
+ *v1 = 1.0; *v2 = 0.0;
+ G_add_d_raster_color_rule ( v1,0,0,0,v2,255,255,255, colors );
+ G_write_colors ( parm.errors->answer, G_mapset (), colors );
+ }
+
+ if ( parm.labels->answer ) {
+ G_init_cats ( num_centers, "Name of dominant center", &pcats );
+ /* write map labels */
+ labels = G_malloc ( num_centers * sizeof (char* ) );
+ GVT_rewind ( pmap ); i = 0;
+ while ( GVT_next ( pmap ) ) {
+ labels[i] = GVT_get_string ( parm.labels->answer, pmap );
+ i ++;
+ }
+ /* add labels to output map categories */
+ error = 0;
+ for ( i=0; i < num_centers; i ++ ) {
+ error = G_set_cat ( cats[i], labels[i], &pcats );
+ if ( error == -1 ) {
+ G_warning (_("Error writing category label '%s': too many categories.\n"), labels[i] );
+ }
+ if ( error == 0 ) {
+ G_warning (_("Error writing category label '%s': got NULL category.\n"), labels[i] );
+ }
+ }
+ if ( error == 1 ) {
+ error = G_write_cats ( parm.output->answer, &pcats );
+ if ( error != 1 ) {
+ G_warning (_("Error writing categories for output map '%s'.\n"), parm.output->answer );
+ }
+ if ( parm.second->answer ) {
+ error = G_write_cats ( parm.second->answer, &pcats );
+ if ( error != 1 ) {
+ G_warning (_("Error writing categories for 'second highest I' map '%s'.\n"), parm.second->answer );
+ }
+ }
+ }
+ G_free_cats ( &pcats );
+ }
+
+ /* write raster map history */
+ G_short_history ( parm.output->answer, "raster", &hist );
+ sprintf ( hist.title, "XTENT: ID of center with highest local I" );
+ error = G_write_history (parm.output->answer, &hist);
+ if ( error == -1 ) {
+ G_warning (_("Could not write history file for output raster map '%s'.\n"), parm.output->answer );
+ }
+
+ if ( parm.second->answer ) {
+ G_short_history ( parm.second->answer, "raster", &hist );
+ G_command_history ( &hist );
+ sprintf ( hist.title, "XTENT: ID of center with second highest local I" );
+ error = G_write_history (parm.second->answer, &hist);
+ if ( error == -1 ) {
+ G_warning (_("Could not write history file for 'competitor ID' raster map '%s'.\n"), parm.second->answer );
+ }
+ }
+
+ if ( parm.errors->answer ) {
+ G_short_history ( parm.errors->answer, "raster", &hist );
+ G_command_history ( &hist );
+ sprintf ( hist.title, "XTENT: competition strength" );
+ error = G_write_history (parm.errors->answer, &hist);
+ if ( error == -1 ) {
+ G_warning (_("Could not write history file for 'errors' raster map '%s'.\n"), parm.errors->answer );
+ }
+ }
+
+
+ /* make colour map for output maps if RGB attribute given */
+ if ( make_colors ) {
+ if ( parm.errors->answer == NULL ) {
+ /* we don't have a color table yet, because there is no error map */
+ colors = G_malloc ( sizeof (struct Colors));
+ }
+ G_init_colors ( colors );
+ GVT_rewind ( pmap );
+ i = 0;
+ while ( GVT_next ( pmap ) ) {
+ if ( sscanf ( GVT_get_string ( parm.rgbcol->answer, pmap),"%i:%i:%i", &r, &g, &b ) != 3 ) {
+ G_warning ( "Malformed RGB attribute string in record no. %i.\n", i + 1 );
+ } else {
+ G_set_color ( (CELL) cats[i], r, g, b, colors );
+ }
+ i ++;
+ }
+ G_write_colors ( parm.output->answer, G_mapset (), colors );
+ if ( parm.second->answer ) {
+ G_write_colors ( parm.second->answer, G_mapset (), colors );
+ }
+ }
+
+ /* if user requested a report file to be written: do so now! */
+ if ( parm.report->answer ) {
+ write_report ( cats, labels, argc, argv );
+ }
+
+ /* warn, if model output is entirely empty! */
+ if ( flag.strict->answer ) {
+ if ( valid_I == 0 ) {
+ G_warning ("All cells in output map are NULL. You need to adjust your model parameters (C,a,k) or do not run in 'strict' mode.\n");
+ }
+ }
+
+ /* TODO: free colors struct */
+
+ /* free cell handles */
+ G_free ( cell );
+ if ( parm.second->answer ) {
+ G_free ( cell_second );
+ }
+ if ( parm.errors->answer ) {
+ G_free ( dcell );
+ }
+
+ /* free other arrays */
+ G_free ( cats );
+
+ /* close input vector map 'centers' */
+ GVT_close_map ( pmap );
+ GVT_free_map ( pmap );
+
+ /* THAT'S IT, FOLKS! */
+ if ( PROGRESS ) {
+ G_message ("\nJOB DONE.\n");
+ }
+
+ return (0);
+}
Added: grass-addons/raster/r.xtent/r.xtent.001.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.001.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.002.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.002.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.003.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.003.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.004.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.004.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.005.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.005.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.006.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.006.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/r.xtent.007.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.xtent/r.xtent.007.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.xtent/tools.c
===================================================================
--- grass-addons/raster/r.xtent/tools.c (rev 0)
+++ grass-addons/raster/r.xtent/tools.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,138 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdarg.h>
+
+#include <sys/types.h>
+#include <unistd.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+/* this runs an external command */
+void run_cmd ( char *cmd ) {
+ G_system ( cmd );
+}
+
+/* this runs a mapcalc instruction using r.mapcalc */
+void r_mapcalc ( char* outmap, char *cmd ) {
+ char buf[5000];
+
+ sprintf ( buf, "r.mapcalc %s=%s", outmap, cmd );
+ run_cmd (buf);
+}
+
+
+/* Function quicksort() sorts an array of double numbers
+ * by ascending size.
+ * The size of the array must be passed in 'n'
+ * Sorting is done in place, i.e. *array will be
+ * altered!
+ */
+int compare ( const void *a, const void *b ) {
+ if ( * (double *) a < * (double *) b ) {
+ return ( -1 );
+ }
+ if ( * (double *) a == * (double *) b ) {
+ return ( 0 );
+ }
+ return ( 1 );
+}
+
+void quicksort ( double *array, int n ) {
+
+ if ( array != NULL ) {
+ qsort ( array, (size_t) n, sizeof (double), compare );
+ }
+}
+
+
+/* This function uses the posix standard function
+ * and a prefix to generate a unique name for a
+ * temporary map in the current mapset.
+ * Memory for the new string is allocated and can
+ * be freed by the caller.
+*/
+char *mktmpmap ( char* prefix ) {
+ int pid;
+ static int uniq = 0;
+ char *name;
+
+ name = malloc (1024 * sizeof (char) );
+ pid = getpid();
+ sprintf ( name, "%s%i%i", prefix, pid, uniq );
+ uniq ++;
+
+ return ( name );
+}
+
+
+/* Returns 1 if map exists in current mapset and is readable,
+ * 0 if there is a problem
+ */
+int map_exists ( char *map ) {
+
+ if ( G_find_cell( map, G_mapset( )) == NULL ) {
+ return ( 0 );
+ }
+
+ return ( 1 );
+}
+
+
+/* this function tests if a raster map exists in the current
+ * mapset and can be opened for write access
+ * Returns 1 if OK, 0 if open failed
+ * On failed open, it will also generate a user-defined warning.
+ * */
+int test_map ( char *map, char *msg, ... ) {
+ char buffer[5000];
+ va_list ap;
+ int error;
+
+ va_start(ap,msg);
+ vsprintf(buffer,msg,ap);
+ va_end(ap);
+
+ if ( G_find_cell (map, "" ) == NULL ) {
+ G_warning (_("%s"), buffer );
+ }
+
+ error = G_open_cell_old ( map, G_find_cell (map, "") );
+ if ( error < 0 ) {
+ G_warning (_("%s"), buffer );
+ G_close_cell (error);
+ return ( 0 );
+ } else {
+ G_close_cell ( error );
+ }
+ return ( 1 );
+}
+
+
+/* this function tests if a raster map exists in the current
+ * mapset and can be opened for write access
+ * It aborts the program if there is a problem.
+ * */
+void test_map_f ( char *map, char *msg, ... ) {
+ char buffer[5000];
+ va_list ap;
+ int error;
+
+ va_start(ap,msg);
+ vsprintf(buffer,msg,ap);
+ va_end(ap);
+
+ if ( G_find_cell (map, "" ) == NULL ) {
+ G_fatal_error (_( "%s"), buffer );
+ }
+
+ error = G_open_cell_old ( map, G_find_cell (map, "" ) );
+ if ( error < 0 ) {
+ G_fatal_error (_( "%s"), buffer );
+ } else {
+ G_close_cell ( error );
+ }
+}
+
Added: grass-addons/raster/r.xtent/tools.h
===================================================================
--- grass-addons/raster/r.xtent/tools.h (rev 0)
+++ grass-addons/raster/r.xtent/tools.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,17 @@
+#ifndef TOOLS_H_
+#define TOOLS_H_
+
+void run_cmd ( char *cmd );
+
+void r_mapcalc ( char *outmap, char *cmd );
+
+void quicksort ( double *array, int n );
+
+char *mktmpmap ( char* prefix );
+
+int map_exists ( char *map );
+
+int test_map ( char *map, ... );
+int test_map_f ( char *map, ... );
+
+#endif /*TOOLS_H_*/
Added: grass-addons/raster/r.xtent/xtent.c
===================================================================
--- grass-addons/raster/r.xtent/xtent.c (rev 0)
+++ grass-addons/raster/r.xtent/xtent.c 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,158 @@
+/*
+ * FUNCTIONS FOR CALCULATION OF THE ORIGINAL XTENT MODEL AND A NORMALIZED VARIATION
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+
+#include "globals.h"
+#include "tools.h"
+#include "gt_vector.h"
+#include "geom.h"
+
+
+/* Returns the record number (0..n) of the vector point that is closest
+ to the location specified by x and y as given by the original Xtent model
+
+ costs is an array of doubles that stores the value of the cost surface map(s) for each center at x and y
+
+ Returns -1 on error,
+ id of most influential center (0..n), otherwise
+
+*/
+int get_center_xtent_original ( double x, double y, double *C, double a, double k, double *costs, double *reach, int *ruler, GVT_map_s *map ) {
+ int i;
+ int center;
+ double I;
+ double d;
+ double max;
+ int max_set;
+ int within_reach;
+
+ /* TODO:
+ speed-up: read point coordinates from array instead of vector map
+ better: leave attribute and coordinate caching to GT_vector lib routines
+ */
+ GVT_rewind ( map );
+ max = 0.0;
+ max_set = 0;
+ I = -1.0;
+ center = -1;
+ i = 0;
+ while ( GVT_next ( map ) ) {
+ within_reach = 0;
+ d = get_distance ( x, y , costs, map );
+ I = pow(C[i],a) - (k*d);
+ if ( reach[i] >= 0.0 ) { /* we have a center with constrained reach */
+ if ( d <= reach[i] ) {
+ /* this center is within reach! */
+ within_reach = 1;
+ }
+ } else {
+ /* reach checking was turned off */
+ within_reach = 1;
+ }
+ if ( DEBUG > 1) {
+ fprintf ( stderr, "I = %.10f^%.10f - %.10f*%.10f = %.10f\n", C[i], a, k, d, I );
+ }
+ if ( (!max_set) && (within_reach)) {
+ max = I;
+ max_set = 1;
+ center = ruler[i];
+ } else {
+ if ( (I > max) && (within_reach) ) {
+ max = I;
+ center = ruler[i];
+ }
+ }
+ i ++;
+ }
+
+ if ( flag.strict->answer ) {
+ if ( max < 0 ) {
+ return ( -1 );
+ }
+ }
+
+ return ( center );
+}
+
+
+/* Returns the record number (0..n) of the vector point that is * second * most influential
+ to the location specified by x and y as given by the Xtent model
+
+ *diff is set to the absolute difference between first and second highest I
+
+ Returns -1 on error,
+ id of second most influential center (0..n), otherwise
+*/
+int get_center_xtent_second ( double x, double y, double *C, double a, double k, double *diff, double *costs, int *ruler, GVT_map_s *map ) {
+ int i,j;
+ int second;
+ int matched;
+ double *sorted;
+ double *unsorted;
+ double I;
+ double d;
+ double val;
+
+ /* make array to store I values for all centers in input map */
+ sorted = G_malloc ( num_centers * sizeof ( double ) );
+ unsorted = G_malloc ( num_centers * sizeof ( double ) );
+
+ /* step through the input vector map and calculate I for all centers */
+ GVT_rewind ( map );
+ i = 0;
+ while ( GVT_next ( map ) ) {
+ d = get_distance ( x, y , costs, map );
+ I = pow(C[i],a) - (k*d);
+ sorted[i] = I;
+ unsorted[i] = I;
+ i ++;
+ }
+
+ /* now sort all I values ascendingly and pick the second highest one that is NOT ruled by
+ another center (if all are ruled, then there is no competitor and we return -1!) */
+ quicksort ( sorted, num_centers );
+ matched = 0;
+ second = -1;
+ val = -1.0;
+ for ( i = num_centers - 2; i >= 0; i -- ) {
+ /* get ID of the center at this position in the sorted array */
+ for ( j = 0; j < num_centers; j ++ ) {
+ if ( unsorted[j] == sorted[i] ) {
+ second = j;
+ }
+ }
+ if ( (ruler[second] == second) && (!matched) ) {
+ val = sorted [ i ];
+ matched = 1;
+ }
+ }
+ if ( !matched ) {
+ return ( -1 );
+ }
+
+ /* store absolute difference between first and second highest I */
+ *diff = fabs ( sorted[num_centers-1] - sorted[num_centers-2] );
+
+ /* get the ID of the center that corresponds to competitor's ID and return it */
+ for ( i=0; i < num_centers; i++ ) {
+ if ( unsorted[i] == val ) {
+ if ( flag.strict->answer ) {
+ /* in strict mode, second I < 0 is also invalid! */
+ if ( unsorted[i] < 0 ) {
+ return ( -1 );
+ }
+ }
+ return ( i );
+ }
+ }
+
+ return ( -1 );
+}
+
+
Added: grass-addons/raster/r.xtent/xtent.h
===================================================================
--- grass-addons/raster/r.xtent/xtent.h (rev 0)
+++ grass-addons/raster/r.xtent/xtent.h 2009-09-28 14:02:12 UTC (rev 39313)
@@ -0,0 +1,7 @@
+#ifndef XTENT_H_
+#define XTENT_H_
+
+int get_center_xtent_original ( double x, double y, double *C, double a, double k, double *costs, double *reach, int *ruler, GVT_map_s *map );
+int get_center_xtent_second ( double x, double y, double *C, double a, double k, double *diff, double *costs, int *ruler, GVT_map_s *map );
+
+#endif /*XTENT_H_*/
More information about the grass-commit
mailing list