v.random: add option to restrict points to vector areas

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

vector/v.random/main.c

@@ -6,6 +6,7 @@
  * AUTHOR(S):    James Darrell McCauley darrell@mccauley-usa.com
  * 	         http://mccauley-usa.com/
  *               OGR support by Martin Landa <landa.martin gmail.com>
+ *               Area support by Markus Metz
  *
  * PURPOSE:      Randomly generate a 2D/3D GRASS vector points map.
  *
@@ -48,31 +49,57 @@ double myrand(void);
 #if defined(__CYGWIN__) || defined(__APPLE__) || defined(__MINGW32__)
 double drand48()
 {
-    return (rand() / 32767.0);
+    return (rand() / (1.0 * RAND_MAX));
 }
 
 #define srand48(sv) (srand((unsigned)(sv)))
 #endif
 
+/* for qsort */
+
+typedef struct {
+    int i;
+    double size;
+    struct bound_box box;
+} BOX_SIZE;
+
+static int sort_by_size(const void *a, const void *b)
+{
+    BOX_SIZE *as = (BOX_SIZE *)a;
+    BOX_SIZE *bs = (BOX_SIZE *)b;
+    
+    if (as->size < bs->size)
+	return -1;
+
+    return (as->size > bs->size);
+}
+
 int main(int argc, char *argv[])
 {
     char *output, buf[DB_SQL_MAX];
     double (*rng) ();
     double max, zmin, zmax;
     int seed;
-    int i, n, b, type, usefloat;
-    struct Map_info Out;
+    int i, j, k, n, b, type, usefloat;
+    int area, nareas, field;
+    struct boxlist *List = NULL;
+    BOX_SIZE *size_list = NULL;
+    int alloc_size_list = 0;
+    struct Map_info In, Out;
     struct line_pnts *Points;
     struct line_cats *Cats;
+    struct cat_list *cat_list;
+    struct bound_box box;
     struct Cell_head window;
     struct GModule *module;
     struct
     {
-	struct Option *output, *nsites, *zmin, *zmax, *zcol, *ztype, *seed;
+	struct Option *input, *field, *cats, *where, *output, *nsites,
+		      *zmin, *zmax, *zcol, *ztype, *seed;
     } parm;
     struct
     {
-	struct Flag *rand, *drand48, *z, *notopo;
+	struct Flag *rand, *drand48, *z, *notopo, *a;
     } flag;
     struct field_info *Fi;
     dbDriver *driver;
@@ -86,7 +113,7 @@ int main(int argc, char *argv[])
     G_add_keyword(_("sampling"));
     G_add_keyword(_("statistics"));
     G_add_keyword(_("random"));
-    module->description = _("Generates randomly 2D/3D vector points map.");
+    module->description = _("Generates random 2D/3D vector points.");
 
     parm.output = G_define_standard_option(G_OPT_V_OUTPUT);
 
@@ -96,6 +123,20 @@ int main(int argc, char *argv[])
     parm.nsites->required = YES;
     parm.nsites->description = _("Number of points to be created");
 
+    parm.input = G_define_standard_option(G_OPT_V_INPUT);
+    parm.input->required = NO;
+    parm.input->description = _("Restrict points to areas in input vector");
+    parm.input->guisection = _("Selection");
+
+    parm.field = G_define_standard_option(G_OPT_V_FIELD_ALL);
+    parm.field->guisection = _("Selection");
+
+    parm.cats = G_define_standard_option(G_OPT_V_CATS);
+    parm.cats->guisection = _("Selection");
+    
+    parm.where = G_define_standard_option(G_OPT_DB_WHERE);
+    parm.where->guisection = _("Selection");
+
     parm.zmin = G_define_option();
     parm.zmin->key = "zmin";
     parm.zmin->type = TYPE_DOUBLE;
@@ -142,6 +183,10 @@ int main(int argc, char *argv[])
     flag.z->description = _("Create 3D output");
     flag.z->guisection = _("3D output");
 
+    flag.a = G_define_flag();
+    flag.a->key = 'a';
+    flag.a->description = _("Generate n points for each individual area");
+
     flag.drand48 = G_define_flag();
     flag.drand48->key = 'd';
     flag.drand48->description = _("Use drand48() function instead of rand()");
@@ -162,12 +207,45 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Number of points must be > 0 (%d given)"), n);
     }
 
+    nareas = 0;
+    cat_list = NULL;
+    field = -1;
+    if (parm.input->answer) {
+	Vect_set_open_level(2); /* topology required */
+	if (2 > Vect_open_old2(&In, parm.input->answer, "", parm.field->answer))
+	    G_fatal_error(_("Unable to open vector map <%s>"),
+			  parm.input->answer);
+
+	if (parm.field->answer)
+	    field = Vect_get_field_number(&In, parm.field->answer);
+
+	if ((parm.cats->answer || parm.where->answer) && field == -1) {
+	    G_warning(_("Invalid layer number (%d). Parameter '%s' or '%s' specified, assuming layer '1'."),
+		      field, parm.cats->key, parm.where->key);
+	    field = 1;
+	}
+	if (field > 0)
+	    cat_list = Vect_cats_set_constraint(&In, field, parm.where->answer,
+						parm.cats->answer);
+	nareas = Vect_get_num_areas(&In);
+	if (nareas == 0) {
+	    Vect_close(&In);
+	    G_fatal_error(_("No areas in vector map <%s>"), parm.input->answer);
+	}
+    }
+    else {
+	if (flag.a->answer)
+	    G_fatal_error(_("The <-%c> flag requires an input vector with areas"),
+	                  flag.a->key);
+    }
+
     /* create new vector map */
     if (-1 == Vect_open_new(&Out, output, flag.z->answer ? WITH_Z : WITHOUT_Z))
         G_fatal_error(_("Unable to create vector map <%s>"), output);
     Vect_set_error_handler_io(NULL, &Out);
 
     /* Do we need to write random values into attribute table? */
+    usefloat = -1;
     if (parm.zcol->answer) {
 	Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
 	driver =
@@ -214,7 +292,6 @@ int main(int argc, char *argv[])
 	}
 
 	type = db_get_column_sqltype(db_get_table_column(table, 1));
-	usefloat = -1;
 	if (type == DB_SQL_TYPE_SMALLINT || type == DB_SQL_TYPE_INTEGER)
 	    usefloat = 0;
 	if (type == DB_SQL_TYPE_REAL || type == DB_SQL_TYPE_DOUBLE_PRECISION)
@@ -248,53 +325,303 @@ int main(int argc, char *argv[])
 
     G_get_window(&window);
 
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    if (nareas > 0) {
+	int first = 1, count;
+	struct bound_box abox, bbox;
+
+	box.W = window.west;
+	box.E = window.east;
+	box.S = window.south;
+	box.N = window.north;
+	box.B = -PORT_DOUBLE_MAX;
+	box.T = PORT_DOUBLE_MAX;
+
+	count = 0;
+
+	for (i = 1; i <= nareas; i++) {
+	    
+	    if (!Vect_get_area_centroid(&In, i))
+		continue;
+
+	    if (field > 0) {
+		if (Vect_get_area_cats(&In, i, Cats))
+		    continue;
+
+		if (!Vect_cats_in_constraint(Cats, field, cat_list))
+		    continue;
+	    }
+
+	    Vect_get_area_box(&In, i, &abox);
+	    if (!Vect_box_overlap(&abox, &box))
+		continue;
+
+	    if (first) {
+		Vect_box_copy(&bbox, &abox);
+		first = 0;
+	    }
+	    else
+		Vect_box_extend(&bbox, &abox);
+	    count++;
+	}
+	if (count == 0) {
+	    Vect_close(&In);
+	    Vect_close(&Out);
+	    Vect_delete(output);
+	    G_fatal_error(_("Selected areas in input vector <%s> do not overlap with the current region"),
+			  parm.input->answer);
+	}
+	Vect_box_copy(&box, &bbox);
+
+	/* does the vector overlap with the current region ? */
+	if (box.W >= window.east || box.E <= window.west ||
+	    box.S >= window.north || box.N <= window.south) {
+
+	    Vect_close(&In);
+	    Vect_close(&Out);
+	    Vect_delete(output);
+	    G_fatal_error(_("Input vector <%s> does not overlap with the current region"),
+	                  parm.input->answer);
+	}
+
+	/* try to reduce the current region */
+	if (window.east > box.E)
+	    window.east = box.E;
+	if (window.west < box.W)
+	    window.west = box.W;
+	if (window.north > box.N)
+	    window.north = box.N;
+	if (window.south < box.S)
+	    window.south = box.S;
+
+	List = Vect_new_boxlist(1);
+	alloc_size_list = 10;
+	size_list = G_malloc(alloc_size_list * sizeof(BOX_SIZE));
+    }
+
+    zmin = zmax = 0;
     if (flag.z->answer || parm.zcol->answer) {
 	zmax = atof(parm.zmax->answer);
 	zmin = atof(parm.zmin->answer);
     }
 
-    Points = Vect_new_line_struct();
-    Cats = Vect_new_cats_struct();
-
     G_message(_("Generating points..."));
-    for (i = 0; i < n; ++i) {
-	double x, y, z;
+    if (flag.a->answer && nareas > 0) {
+	struct bound_box abox, bbox;
+	int cat = 1;
 
-	G_percent(i, n, 5);
+	/* n points for each area */
+	nareas = Vect_get_num_areas(&In);
+	
+	G_percent(0, nareas, 1);
+	for (area = 1; area <= nareas; area++) {
 
-	Vect_reset_line(Points);
-	Vect_reset_cats(Cats);
+	    G_percent(area, nareas, 1);
 
-	x = rng() / max * (window.west - window.east) + window.east;
-	y = rng() / max * (window.north - window.south) + window.south;
-	z = rng() / max * (zmax - zmin) + zmin;
-        
-	if (flag.z->answer)
-	    Vect_append_point(Points, x, y, z);
-	else
-	    Vect_append_point(Points, x, y, 0.0);
-
-	if (parm.zcol->answer) {
-	    sprintf(buf, "insert into %s values ( %d, ", Fi->table, i + 1);
-	    db_set_string(&sql, buf);
-	    /* Round random value if column is integer type */
-	    if (usefloat)
-		sprintf(buf, "%f )", z);
-	    else
-		sprintf(buf, "%.0f )", z);
-	    db_append_string(&sql, buf);
+	    if (!Vect_get_area_centroid(&In, area))
+		continue;
+
+	    if (field > 0) {
+		if (Vect_get_area_cats(&In, area, Cats))
+		    continue;
+
+		if (!Vect_cats_in_constraint(Cats, field, cat_list)) {
+		    continue;
+		}
+	    }
 
-	    G_debug(3, db_get_string(&sql));
-	    if (db_execute_immediate(driver, &sql) != DB_OK) {
-		G_fatal_error(_("Cannot insert new row: %s"),
-			      db_get_string(&sql));
+	    box.W = window.west;
+	    box.E = window.east;
+	    box.S = window.south;
+	    box.N = window.north;
+	    box.B = -PORT_DOUBLE_MAX;
+	    box.T = PORT_DOUBLE_MAX;
+	    
+	    Vect_get_area_box(&In, area, &abox);
+	    if (!Vect_box_overlap(&box, &abox))
+		continue;
+		
+	    bbox = abox;
+	    if (bbox.W < box.W)
+		bbox.W = box.W;
+	    if (bbox.E > box.E)
+		bbox.E = box.E;
+	    if (bbox.S < box.S)
+		bbox.S = box.S;
+	    if (bbox.N > box.N)
+		bbox.N = box.N;
+
+	    for (i = 0; i < n; ++i) {
+		double x, y, z;
+		int outside = 1;
+		int ret;
+
+		Vect_reset_line(Points);
+		Vect_reset_cats(Cats);
+
+		while (outside) {
+		    x = rng() / max * (bbox.W - bbox.E) + bbox.E;
+		    y = rng() / max * (bbox.N - bbox.S) + bbox.S;
+		    z = rng() / max * (zmax - zmin) + zmin;
+
+		    ret = Vect_point_in_area(x, y, &In, area, &abox);
+
+		    G_debug(3, "    area = %d Vect_point_in_area() = %d", area, ret);
+
+		    if (ret >= 1) {
+			outside = 0;
+		    }
+		}
+
+		if (flag.z->answer)
+		    Vect_append_point(Points, x, y, z);
+		else
+		    Vect_append_point(Points, x, y, 0.0);
+
+		if (parm.zcol->answer) {
+		    sprintf(buf, "insert into %s values ( %d, ", Fi->table, i + 1);
+		    db_set_string(&sql, buf);
+		    /* Round random value if column is integer type */
+		    if (usefloat)
+			sprintf(buf, "%f )", z);
+		    else
+			sprintf(buf, "%.0f )", z);
+		    db_append_string(&sql, buf);
+
+		    G_debug(3, db_get_string(&sql));
+		    if (db_execute_immediate(driver, &sql) != DB_OK) {
+			G_fatal_error(_("Cannot insert new row: %s"),
+				      db_get_string(&sql));
+		    }
+		}
+
+		Vect_cat_set(Cats, 1, cat++);
+		Vect_write_line(&Out, GV_POINT, Points, Cats);
 	    }
 	}
+    }
+    else {
+	/* n points in total */
+	for (i = 0; i < n; ++i) {
+	    double x, y, z;
+
+	    G_percent(i, n, 4);
+
+	    Vect_reset_line(Points);
+	    Vect_reset_cats(Cats);
+
+	    x = rng() / max * (window.west - window.east) + window.east;
+	    y = rng() / max * (window.north - window.south) + window.south;
+	    z = rng() / max * (zmax - zmin) + zmin;
+	    
+	    if (nareas) {
+		int outside = 1;
+
+		do {
+		    /* select areas by box */
+		    box.E = x;
+		    box.W = x;
+		    box.N = y;
+		    box.S = y;
+		    box.T = PORT_DOUBLE_MAX;
+		    box.B = -PORT_DOUBLE_MAX;
+		    Vect_select_areas_by_box(&In, &box, List);
+		    G_debug(3, "  %d areas selected by box", List->n_values);
+
+		    /* sort areas by size, the smallest is likely to be the nearest */
+		    if (alloc_size_list < List->n_values) {
+			alloc_size_list = List->n_values;
+			size_list = G_realloc(size_list, alloc_size_list * sizeof(BOX_SIZE));
+		    }
+
+		    k = 0;
+		    for (j = 0; j < List->n_values; j++) {
+			area = List->id[j];
+
+			if (!Vect_get_area_centroid(&In, area))
+			    continue;
+
+			if (field > 0) {
+			    if (Vect_get_area_cats(&In, area, Cats))
+				continue;
+
+			    if (!Vect_cats_in_constraint(Cats, field, cat_list)) {
+				continue;
+			    }
+			}
+
+			List->id[k] = List->id[j];
+			List->box[k] = List->box[j];
+			size_list[k].i = List->id[k];
+			box = List->box[k];
+			size_list[k].box = List->box[k];
+			size_list[k].size = (box.N - box.S) * (box.E - box.W);
+			k++;
+		    }
+		    List->n_values = k;
+		    
+		    if (List->n_values == 2) {
+			/* simple swap */
+			if (size_list[1].size < size_list[0].size) {
+			    size_list[0].i = List->id[1];
+			    size_list[1].i = List->id[0];
+			    size_list[0].box = List->box[1];
+			    size_list[1].box = List->box[0];
+			}
+		    }
+		    else if (List->n_values > 2)
+			qsort(size_list, List->n_values, sizeof(BOX_SIZE), sort_by_size);
+
+		    for (j = 0; j < List->n_values; j++) {
+			int ret;
+
+			area = size_list[j].i;
+			ret = Vect_point_in_area(x, y, &In, area, &size_list[j].box);
+
+			G_debug(3, "    area = %d Vect_point_in_area() = %d", area, ret);
+
+			if (ret >= 1) {
+			    outside = 0;
+			    break;
+			}
+		    }
+		    if (outside) {
+			x = rng() / max * (window.west - window.east) + window.east;
+			y = rng() / max * (window.north - window.south) + window.south;
+			z = rng() / max * (zmax - zmin) + zmin;
+		    }
+		} while (outside);
+	    }
+
+	    if (flag.z->answer)
+		Vect_append_point(Points, x, y, z);
+	    else
+		Vect_append_point(Points, x, y, 0.0);
+
+	    if (parm.zcol->answer) {
+		sprintf(buf, "insert into %s values ( %d, ", Fi->table, i + 1);
+		db_set_string(&sql, buf);
+		/* Round random value if column is integer type */
+		if (usefloat)
+		    sprintf(buf, "%f )", z);
+		else
+		    sprintf(buf, "%.0f )", z);
+		db_append_string(&sql, buf);
+
+		G_debug(3, db_get_string(&sql));
+		if (db_execute_immediate(driver, &sql) != DB_OK) {
+		    G_fatal_error(_("Cannot insert new row: %s"),
+				  db_get_string(&sql));
+		}
+	    }
 
-	Vect_cat_set(Cats, 1, i + 1);
-	Vect_write_line(&Out, GV_POINT, Points, Cats);
+	    Vect_cat_set(Cats, 1, i + 1);
+	    Vect_write_line(&Out, GV_POINT, Points, Cats);
+	}
+	G_percent(1, 1, 1);
     }
-    G_percent(1, 1, 1);
     
     if (parm.zcol->answer) {
 	db_commit_transaction(driver);

vector/v.random/v.random.html

@@ -6,15 +6,30 @@ current region using the selected random number generator.
 <p><em>v.random</em> can generate also 3D vector points or
 write random value to attribute table. Point height range or
 attribute value range is controlled by specifying zmin and zmax values.
-Both z values are included in range (<em>zmin >= z <= zmax</em>).
+Both z values are included in range (<em>zmin <= z <= zmax</em>).
 Generated random attribute value type can be controlled by column
 data type. Use <b>INTEGER</b> column type for integers and 
 <b>DOUBLE PRECISION</b> for floating point numbers. Integer values are
 calculated by rounding random floating point number.
 
-<p>To produce repeatable results a random seed can be set using the option <em>seed</em>>.
-    
+<p>To produce repeatable results a random seed can be set using the
+option <em>seed</em>.
+
+<h3>Restriction to vector areas</h3>
+<p>
+If an <em>input</em> vector map with areas is specified, the location of 
+random points is restricted to the selected areas. By default, the 
+requested number of points are distributed across all areas.
+
+<p>
+If the <em>-a</em> flag is given, the requested number of points is 
+generated for each individual area. For example, if 20 points should be 
+generated and the input map has 100 individual areas, 2000 points will 
+be generated in total.
+
 <h2>EXAMPLES</h2>
+All examples use the North Carolina sample dataset.
+<p>
 Generate 20 random points with binary attribute (only 0 or 1):
 <div class="code"><pre>
 v.random output=binary_random n=20 zmin=0 zmax=1 column='binary INTEGER'
@@ -27,18 +42,35 @@ v.random seed=501 output=height_random n=20 zmin=0 zmax=1 column=height
 
 Get 20 random samples from raster map:
 <div class="code"><pre>
+g.region -p rast=elevation
 v.random output=random_samples n=20
 v.db.addtable map=random_samples layer=1 columns='cat INTEGER, sample DOUBLE PRECISION'
-v.what.rast vector=random_samples raster=elevation.10m@PERMANENT layer=1 column=sample 
+v.what.rast vector=random_samples raster=elevation@PERMANENT layer=1 column=sample 
 </pre></div>
 
 Generate 20 random points and sample attribute data from geology (vector) map:
 <div class="code"><pre>
+g.region -p vect=geology
 v.random output=random_samples n=20
 v.db.addtable map=random_samples layer=1 columns='cat integer, geology varchar(100)'
 v.what.vect vector=random_samples layer=1 column=geology qvector=geology@PERMANENT qlayer=1 qcolumn=label 
 </pre></div>
 
+Generate 20 random points in forested areas
+<div class="code"><pre>
+g.region -p rast=landclass96
+r.to.vect -v input=landclass96 output=landclass96 type=area
+v.random input=landclass96 output=random_samples n=20 where="label = 'forest'"
+</pre></div>
+
+Generate 20 random points in each forest patch
+<div class="code"><pre>
+g.region -p rast=landclass96
+r.to.vect -v input=landclass96 output=landclass96 type=area
+v.random input=landclass96 output=random_samples n=20 where="label = 'forest'" -a
+</pre></div>
+
+
 <h2>SEE ALSO</h2>
 
 UNIX man pages for <em>rand(3)</em> and <em>drand48(3)</em>.