[GRASS-SVN] r59981 - grass/trunk/vector/v.random
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Apr 30 07:09:56 PDT 2014
Author: mmetz
Date: 2014-04-30 07:09:56 -0700 (Wed, 30 Apr 2014)
New Revision: 59981
Modified:
grass/trunk/vector/v.random/main.c
grass/trunk/vector/v.random/v.random.html
Log:
v.random: add option to restrict points to vector areas
Modified: grass/trunk/vector/v.random/main.c
===================================================================
--- grass/trunk/vector/v.random/main.c 2014-04-30 13:58:17 UTC (rev 59980)
+++ grass/trunk/vector/v.random/main.c 2014-04-30 14:09:56 UTC (rev 59981)
@@ -6,6 +6,7 @@
* AUTHOR(S): James Darrell McCauley darrell at 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 @@
#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 @@
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 @@
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 @@
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 @@
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 @@
}
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 @@
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 (!Vect_get_area_centroid(&In, area))
+ continue;
- 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);
+ if (field > 0) {
+ if (Vect_get_area_cats(&In, area, Cats))
+ continue;
+
+ if (!Vect_cats_in_constraint(Cats, field, cat_list)) {
+ continue;
+ }
+ }
+
+ 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
- sprintf(buf, "%.0f )", z);
- db_append_string(&sql, buf);
+ Vect_append_point(Points, x, y, 0.0);
- 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));
+ 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);
Modified: grass/trunk/vector/v.random/v.random.html
===================================================================
--- grass/trunk/vector/v.random/v.random.html 2014-04-30 13:58:17 UTC (rev 59980)
+++ grass/trunk/vector/v.random/v.random.html 2014-04-30 14:09:56 UTC (rev 59981)
@@ -6,15 +6,30 @@
<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 @@
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 at PERMANENT layer=1 column=sample
+v.what.rast vector=random_samples raster=elevation at 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 at 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>.
More information about the grass-commit
mailing list