1
2 /****************************************************************************
3 *
4 * MODULE: v.perturb
5 * AUTHOR(S): James Darrell McCauley darrell@mccauley-usa.com
6 * http://mccauley-usa.com/
7 * PURPOSE: Random location perturbations of vector points
8 *
9 * COPYRIGHT: (C) 1994-2009 by James Darrell McCauley and the GRASS Development Team
10 *
11 * Modification History:
12 * 3/2006 added min and seed MN/SM ITC-irst
13 * 2005 updated to GRASS 6 RB ITC-irst
14 * 0.1B <18 Feb 1994> first version (jdm)
15 * 0.2B <02 Jan 1995> clean man page, added html page (jdm)
16 * 0.3B <25 Feb 1995> cleaned up 'gcc -Wall' warnings (jdm)
17 * <13 Sept 2000> released under GPL
18 *
19 * This program is free software under the GNU General
20 * Public License (>=v2). Read the file COPYING that
21 * comes with GRASS for details.
22 *
23 *****************************************************************************/
24
25 #include <stdlib.h>
26 #include <math.h>
27
28 #include <grass/gis.h>
29 #include <grass/dbmi.h>
30 #include <grass/vector.h>
31 #include <grass/glocale.h>
32
33 #include "perturb.h"
34 #include "zufall.h"
35
main(int argc,char ** argv)36 int main(int argc, char **argv)
37 {
38 double p1, p2, numbers[1000], numbers2[1000];
39 int (*rng) ();
40 int i;
41 int line, nlines, ttype, n, ret, seed, field;
42 struct field_info *Fi, *Fin;
43 double min = 0.;
44 int debuglevel = 3;
45
46 struct Map_info In, Out;
47
48 struct line_pnts *Points;
49 struct line_cats *Cats;
50
51 struct Cell_head window;
52 struct GModule *module;
53 struct
54 {
55 struct Option *in, *out, *dist, *pars, *min, *seed, *field;
56 struct Flag *no_topo;
57 } parm;
58
59 G_gisinit(argv[0]);
60
61 module = G_define_module();
62 G_add_keyword(_("vector"));
63 G_add_keyword(_("geometry"));
64 G_add_keyword(_("statistics"));
65 G_add_keyword(_("random"));
66 G_add_keyword(_("point pattern"));
67 G_add_keyword(_("level1"));
68
69 module->description =
70 _("Random location perturbations of vector points.");
71
72 parm.in = G_define_standard_option(G_OPT_V_INPUT);
73
74 parm.field = G_define_standard_option(G_OPT_V_FIELD_ALL);
75
76 parm.out = G_define_standard_option(G_OPT_V_OUTPUT);
77
78 parm.dist = G_define_option();
79 parm.dist->key = "distribution";
80 parm.dist->type = TYPE_STRING;
81 parm.dist->required = NO;
82 parm.dist->options = "uniform,normal";
83 parm.dist->answer = "uniform";
84 parm.dist->description = _("Distribution of perturbation");
85
86 parm.pars = G_define_option();
87 parm.pars->key = "parameters";
88 parm.pars->type = TYPE_DOUBLE;
89 parm.pars->required = YES;
90 parm.pars->multiple = YES;
91 parm.pars->label = _("Parameter(s) of distribution");
92 parm.pars->description = _("If the distribution "
93 "is uniform, only one parameter, the maximum, is needed. "
94 "For a normal distribution, two parameters, the mean and "
95 "standard deviation, are required.");
96
97 parm.min = G_define_option();
98 parm.min->key = "minimum";
99 parm.min->type = TYPE_DOUBLE;
100 parm.min->required = NO;
101 parm.min->answer = "0.0";
102 parm.min->description = _("Minimum deviation in map units");
103
104 parm.seed = G_define_option();
105 parm.seed->key = "seed";
106 parm.seed->type = TYPE_INTEGER;
107 parm.seed->required = NO;
108 parm.seed->answer = "0";
109 parm.seed->description = _("Seed for random number generation");
110
111 parm.no_topo = G_define_standard_flag(G_FLG_V_TOPO);
112
113 if (G_parser(argc, argv))
114 exit(EXIT_FAILURE);
115
116 min = atof(parm.min->answer);
117 seed = atoi(parm.seed->answer);
118
119 switch (parm.dist->answer[0]) {
120 case 'u':
121 rng = zufall;
122 break;
123 case 'n':
124 rng = normalen;
125 break;
126 }
127 if (rng == zufall) {
128 i = sscanf(parm.pars->answer, "%lf", &p1);
129 if (i != 1) {
130 G_fatal_error(_("Error scanning arguments"));
131 }
132 else if (p1 <= 0)
133 G_fatal_error(_("Maximum of uniform distribution must be >= zero"));
134 }
135 else {
136 if ((i = sscanf(parm.pars->answer, "%lf,%lf", &p1, &p2)) != 2) {
137 G_fatal_error(_("Error scanning arguments"));
138 }
139 if (p2 <= 0)
140 G_fatal_error(_("Standard deviation of normal distribution must be >= zero"));
141 }
142
143 G_get_window(&window);
144
145 /* Open input */
146 Vect_set_open_level(2);
147 if (Vect_open_old_head2(&In, parm.in->answer, "", parm.field->answer) < 0)
148 G_fatal_error(_("Unable to open vector map <%s>"), parm.in->answer);
149
150 field = Vect_get_field_number(&In, parm.field->answer);
151
152 /* Open output */
153 if (Vect_open_new(&Out, parm.out->answer, WITHOUT_Z) < 0) /* TODO add z support ? */
154 G_fatal_error(_("Unable to create vector map <%s>"), parm.out->answer);
155
156 Vect_hist_copy(&In, &Out);
157 Vect_hist_command(&Out);
158
159 /* Generate a bunch of random numbers */
160 zufalli(&seed);
161 myrng(numbers, 1000, rng, p1 - min, p2);
162 myrng(numbers2, 1000, rng, p1, p2);
163
164 Points = Vect_new_line_struct();
165 Cats = Vect_new_cats_struct();
166
167 nlines = Vect_get_num_lines(&In);
168
169 /* Close input, re-open on level 1 */
170 Vect_close(&In);
171 Vect_set_open_level(1);
172 if (Vect_open_old2(&In, parm.in->answer, "", parm.field->answer) < 0)
173 G_fatal_error(_("Unable to open vector map <%s>"), parm.in->answer);
174
175 i = 0;
176 line = 0;
177 while (1) {
178 int type = Vect_read_next_line(&In, Points, Cats);
179
180 /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
181 if (type == -1) {
182 G_fatal_error(_("Unable to read vector map"));
183 }
184 else if (type == -2) {
185 break;
186 }
187 G_percent(line++, nlines, 4);
188
189 if (type & GV_POINT) {
190 if (field != -1 && !Vect_cat_get(Cats, field, NULL))
191 continue;
192
193 if (i >= 800) {
194 /* Generate some more random numbers */
195 myrng(numbers, 1000, rng, p1 - min, p2);
196 myrng(numbers2, 1000, rng, p1, p2);
197 i = 0;
198 }
199
200 G_debug(debuglevel, "x: %f y: %f", Points->x[0],
201 Points->y[0]);
202
203 /* perturb */
204 /* TODO: tends to concentrate in box corners when min is used */
205 if (numbers2[i] >= 0) {
206 if (numbers[i] >= 0) {
207 G_debug(debuglevel, "deltax: %f", numbers[i] + min);
208 Points->x[0] += numbers[i++] + min;
209 }
210 else {
211 G_debug(debuglevel, "deltax: %f", numbers[i] - min);
212 Points->x[0] += numbers[i++] - min;
213 }
214 Points->y[0] += numbers2[i++];
215 }
216 else {
217 if (numbers[i] >= 0) {
218 G_debug(debuglevel, "deltay: %f", numbers[i] + min);
219 Points->y[0] += numbers[i++] + min;
220 }
221 else {
222 G_debug(debuglevel, "deltay: %f", numbers[i] - min);
223 Points->y[0] += numbers[i++] - min;
224 }
225 Points->x[0] += numbers2[i++];
226 }
227
228 G_debug(debuglevel, "x_pert: %f y_pert: %f", Points->x[0],
229 Points->y[0]);
230 }
231
232 Vect_write_line(&Out, type, Points, Cats);
233 }
234 G_percent(1, 1, 1);
235
236 /* Copy tables */
237 n = Vect_get_num_dblinks(&In);
238 ttype = GV_1TABLE;
239 if (n > 1)
240 ttype = GV_MTABLE;
241
242 for (i = 0; i < n; i++) {
243 Fi = Vect_get_dblink(&In, i);
244 if (Fi == NULL) {
245 G_fatal_error(_("Cannot get db link info"));
246 }
247 Fin = Vect_default_field_info(&Out, Fi->number, Fi->name, ttype);
248 Vect_map_add_dblink(&Out, Fi->number, Fi->name, Fin->table, Fi->key,
249 Fin->database, Fin->driver);
250
251 ret = db_copy_table(Fi->driver, Fi->database, Fi->table,
252 Fin->driver, Vect_subst_var(Fin->database, &Out),
253 Fin->table);
254 if (ret == DB_FAILED) {
255 G_warning("Cannot copy table");
256 }
257 }
258
259 Vect_close(&In);
260
261 if (!parm.no_topo->answer)
262 Vect_build(&Out);
263 Vect_close(&Out);
264
265 return (EXIT_SUCCESS);
266 }
267