[GRASS-SVN] r34528 - grass/trunk/vector/v.drape
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Nov 27 06:16:49 EST 2008
Author: marisn
Date: 2008-11-27 06:16:49 -0500 (Thu, 27 Nov 2008)
New Revision: 34528
Modified:
grass/trunk/vector/v.drape/main.c
grass/trunk/vector/v.drape/v.drape.html
Log:
v.drape WHERE support and ignore vectors without height information (NULL in raster) (merge from develbranch_6 r34150)
Modified: grass/trunk/vector/v.drape/main.c
===================================================================
--- grass/trunk/vector/v.drape/main.c 2008-11-27 10:18:44 UTC (rev 34527)
+++ grass/trunk/vector/v.drape/main.c 2008-11-27 11:16:49 UTC (rev 34528)
@@ -4,7 +4,7 @@
* MODULE: v.drape
*
* AUTHOR(S): Radim Blazek, Dylan Beaudette
- *
+ * Maris Nartiss - WHERE support and raster NULL support
* PURPOSE: Convert 2D vector to 3D vector by sampling of elevation raster.
*
* COPYRIGHT: (C) 2005 by the GRASS Development Team
@@ -29,32 +29,137 @@
*/
#include <stdlib.h>
-#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
+#include <grass/dbmi.h>
#include <grass/glocale.h>
+/* Samples raster map */
+int sample_raster(const int ltype, int fdrast, struct Cell_head window,
+ struct line_pnts *Points, const INTERP_TYPE method,
+ const double scale, const struct Option *null_opt,
+ const double null_val)
+{
+ double estimated_elevation;
+ int j;
+
+ /* adjust flow based on specific type of line */
+ switch (ltype) {
+ /* points (at least 1 vertex) */
+ case GV_POINT:
+ case GV_CENTROID:
+ case GV_KERNEL:
+ /* sample raster at this point, and update the z-coordinate
+ * (note that input vector should not be 3D!)
+ */
+ estimated_elevation = scale * G_get_raster_sample(fdrast,
+ &window,
+ NULL,
+ Points->
+ y[0],
+ Points->
+ x[0], 0, method);
+ /* Elevation value has to be meaningfull */
+ if (G_is_d_null_value(&estimated_elevation)) {
+ if (null_opt->answer) {
+ estimated_elevation = null_val;
+ }
+ else {
+ return 0;
+ }
+ }
+ /* update the elevation value for each data point */
+ Points->z[0] = estimated_elevation;
+ break;
+ /* standard lines (at least 2 vertexes) */
+ case GV_LINE:
+ case GV_BOUNDARY:
+ if (Points->n_points < 2)
+ break; /* At least 2 points */
+
+ /* loop through each point in a line */
+ for (j = 0; j < Points->n_points; j++) {
+ /* sample raster at this point, and update the z-coordinate (note that input vector should not be 3D!) */
+ estimated_elevation = scale * G_get_raster_sample(fdrast,
+ &window,
+ NULL,
+ Points->
+ y[j],
+ Points->
+ x[j], 0,
+ method);
+
+ if (G_is_d_null_value(&estimated_elevation)) {
+ if (null_opt->answer) {
+ estimated_elevation = null_val;
+ }
+ else {
+ return 0;
+ }
+ }
+ /* update the elevation value for each data point */
+ Points->z[j] = estimated_elevation;
+ } /* end looping through point in a line */
+ break;
+
+ /* lines with at least 3 vertexes */
+ case GV_FACE:
+ if (Points->n_points < 3)
+ break; /* At least 3 points */
+
+ /* loop through each point in a line */
+ for (j = 0; j < Points->n_points; j++) {
+ /* sample raster at this point, and update the z-coordinate (note that input vector should not be 3D!) */
+ estimated_elevation = scale * G_get_raster_sample(fdrast,
+ &window,
+ NULL,
+ Points->
+ y[j],
+ Points->
+ x[j], 0,
+ method);
+
+ if (G_is_d_null_value(&estimated_elevation)) {
+ if (null_opt->answer) {
+ estimated_elevation = null_val;
+ }
+ else {
+ return 0;
+ }
+ }
+ /* update the elevation value for each data point */
+ Points->z[j] = estimated_elevation;
+ }
+ break;
+ } /* end line type switch */
+ return 1;
+}
+
int main(int argc, char *argv[])
{
struct GModule *module;
struct Option *in_opt, *out_opt, *type_opt, *rast_opt, *method_opt,
- *scale_opt;
+ *scale_opt, *where_opt, *layer_opt, *null_opt;
struct Map_info In, Out;
struct line_pnts *Points;
struct line_cats *Cats;
+ struct field_info *Fi;
- /* int layer; */
- int line, nlines, otype, ltype;
- BOUND_BOX in_bbox, region_bbox, rast_bbox;
+ int line, nlines, otype, ltype, ncats =
+ 0, layer, i, c, *cats, field_index, id, out_num = 0, *new_cats;
char *mapset;
- int j;
- double scale, estimated_elevation;
+ double scale, null_val;
INTERP_TYPE method = UNKNOWN;
int fdrast; /* file descriptor for raster map is int */
- struct Cell_head window, rast_window;
+ struct Cell_head window;
+ dbDriver *driver;
+ dbHandle handle;
+ dbTable *table;
+ dbString table_name, dbsql, valstr;
+
G_gisinit(argv[0]);
module = G_define_module();
@@ -73,13 +178,9 @@
rast_opt->key = "rast";
rast_opt->required = NO;
rast_opt->description = _("Elevation raster map for height extraction");
+
+ out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
- scale_opt = G_define_option();
- scale_opt->key = "scale";
- scale_opt->type = TYPE_DOUBLE;
- scale_opt->description = _("Scale sampled raster values");
- scale_opt->answer = "1.0";
-
method_opt = G_define_option();
method_opt->key = "method";
method_opt->type = TYPE_STRING;
@@ -91,12 +192,36 @@
"bilinear;bilinear interpolation;"
"cubic;cubic convolution interpolation;";
method_opt->description = _("Sampling method");
+
+ scale_opt = G_define_option();
+ scale_opt->key = "scale";
+ scale_opt->type = TYPE_DOUBLE;
+ scale_opt->description = _("Scale sampled raster values");
+ scale_opt->answer = "1.0";
- out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ where_opt = G_define_standard_option(G_OPT_WHERE);
+ layer_opt = G_define_standard_option(G_OPT_V_FIELD);
+ layer_opt->answer = "1";
+ layer_opt->description = _("Layer is only used for WHERE SQL statement");
+
+ null_opt = G_define_option();
+ null_opt->key = "null_value";
+ null_opt->type = TYPE_DOUBLE;
+ null_opt->label = _("Vector Z value for unknown height");
+ null_opt->description =
+ _("Will set Z to this value, if value from raster map can not be read");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+ layer = atoi(layer_opt->answer);
+ if (layer == 0)
+ G_fatal_error(_("Layer 0 not supported"));
+
+ if (null_opt->answer)
+ null_val = atof(null_opt->answer);
+
/* which interpolation method should we use */
if (method_opt->answer[0] == 'b')
method = BILINEAR;
@@ -106,8 +231,6 @@
method = NEAREST;
}
- /* setup the raster for sampling */
-
/* setup the region */
G_get_window(&window);
@@ -127,9 +250,6 @@
G_fatal_error(_("Unable to open raster map <%s>"), rast_opt->answer);
}
- /* read raster header */
- G_get_cellhd(rast_opt->answer, mapset, &rast_window);
-
Vect_set_open_level(2);
/* check input/output vector maps */
@@ -138,37 +258,49 @@
Vect_open_old(&In, in_opt->answer, "");
- /* checks
- does the elevation raster cover the entire are of the vector map?
- does the current region include the entire input vector map ?
- */
- Vect_get_map_box(&In, &in_bbox);
- Vect_region_box(&window, ®ion_bbox);
- Vect_region_box(&rast_window, &rast_bbox);
- if (in_bbox.W < region_bbox.W ||
- in_bbox.E > region_bbox.E ||
- in_bbox.S < region_bbox.S || in_bbox.N > region_bbox.N) {
- G_warning(_("Current region does not include the entire input vector map <%s>"),
- in_opt->answer);
+ if (where_opt->answer) {
+ /* Let's get vector layers db connections information */
+ Fi = Vect_get_field(&In, layer);
+ if (!Fi) {
+ Vect_close(&In);
+ G_fatal_error(_("Database connection not defined for layer %d"),
+ layer);
+ }
+ G_debug(1,
+ "Field number:%d; Name:<%s>; Driver:<%s>; Database:<%s>; Table:<%s>; Key:<%s>",
+ Fi->number, Fi->name, Fi->driver, Fi->database, Fi->table,
+ Fi->key);
+ /* Prepeare strings for use in db_* calls */
+ db_init_string(&dbsql);
+ db_init_string(&valstr);
+ db_init_string(&table_name);
+ db_init_handle(&handle);
+
+ /* Prepearing database for use */
+ driver = db_start_driver(Fi->driver);
+ if (driver == NULL) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to start driver <%s>"), Fi->driver);
+ }
+ db_set_handle(&handle, Fi->database, NULL);
+ if (db_open_database(driver, &handle) != DB_OK) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, driver);
+ }
+ db_set_string(&table_name, Fi->table);
+ if (db_describe_table(driver, &table_name, &table) != DB_OK) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to describe table <%s>"), Fi->table);
+ }
+ ncats =
+ db_select_int(driver, Fi->table, Fi->key, where_opt->answer,
+ &cats);
+ if (ncats < 1)
+ G_fatal_error(_("No features match Your query"));
+ G_debug(3, "Number of features matching query: %d", ncats);
}
- if (in_bbox.W < rast_bbox.W ||
- in_bbox.E > rast_bbox.E ||
- in_bbox.S < rast_bbox.S || in_bbox.N > rast_bbox.N) {
- G_warning(_("Elevation raster map <%s> does not cover the entire area "
- "of the input vector map <%s>. "), rast_opt->answer,
- in_opt->answer);
- }
- /* setup the new vector map */
- /* remember to open the new vector map as 3D */
- Vect_open_new(&Out, out_opt->answer, WITH_Z);
- Vect_copy_head_data(&In, &Out);
- Vect_hist_copy(&In, &Out);
- Vect_hist_command(&Out);
- /* copy the input vector's attribute table to the new vector */
- /* This works for both level 1 and 2 */
- Vect_copy_tables(&In, &Out, 0);
-
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
@@ -177,99 +309,119 @@
(GV_POINTS | GV_LINES | GV_BOUNDARY | GV_CENTROID | GV_FACE |
GV_KERNEL))) {
- /* loop through each line in the dataset */
- nlines = Vect_get_num_lines(&In);
+ /* process all lines matching WHERE statement */
+ if (where_opt->answer) {
+ field_index = Vect_cidx_get_field_index(&In, layer);
+ for (i = 0; i < ncats; i++) {
+ G_percent(i, ncats, 2);
+
+ c = Vect_cidx_find_next(&In, field_index, cats[i],
+ otype, 0, <ype, &id);
+ while (c >= 0) {
+ c++;
+ if (ltype & otype) {
+ if (Vect_read_line(&In, Points, Cats, id) == ltype)
+ if (sample_raster
+ (ltype, fdrast, window, Points, method,
+ scale, null_opt, null_val)) {
+
+ /* Open new file only if there is something to write in */
+ if (out_num < 1) {
+ if (0 >
+ Vect_open_new(&Out, out_opt->answer,
+ WITH_Z)) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to create vector map <%s>"),
+ out_opt->answer);
+ }
+ Vect_copy_head_data(&In, &Out);
+ Vect_hist_copy(&In, &Out);
+ Vect_hist_command(&Out);
+ }
+
+ Vect_write_line(&Out, ltype, Points, Cats);
+ out_num++;
+ }
+ }
+ else
+ G_fatal_error
+ ("Error in Vect_cidx_find_next function! Report a bug.");
+ c = Vect_cidx_find_next(&In, field_index, cats[i],
+ otype, c, <ype, &id);
+ }
+ }
+ }
+ else {
+ /* loop through each line in the dataset */
+ nlines = Vect_get_num_lines(&In);
- for (line = 1; line <= nlines; line++) {
+ for (line = 1; line <= nlines; line++) {
- /* progress feedback */
- G_percent(line, nlines, 2);
+ /* progress feedback */
+ G_percent(line, nlines, 2);
- /* get the line type */
- ltype = Vect_read_line(&In, Points, Cats, line);
+ /* get the line type */
+ ltype = Vect_read_line(&In, Points, Cats, line);
- /* adjust flow based on specific type of line */
- switch (ltype) {
- /* points (at least 1 vertex) */
- case GV_POINT:
- case GV_CENTROID:
- case GV_KERNEL:
- /* sample raster at this point, and update the z-coordinate
- * (note that input vector should not be 3D!)
- */
- estimated_elevation = scale * G_get_raster_sample(fdrast,
- &window,
- NULL,
- Points->
- y[0],
- Points->
- x[0], 0,
- method);
-
- /* update the elevation value for each data point */
- Points->z[0] = estimated_elevation;
- break;
- /* standard lines (at least 2 vertexes) */
- case GV_LINE:
- case GV_BOUNDARY:
- if (Points->n_points < 2)
- break; /* At least 2 points */
-
- /* loop through each point in a line */
- for (j = 0; j < Points->n_points; j++) {
- /* sample raster at this point, and update the z-coordinate (note that input vector should not be 3D!) */
- estimated_elevation = scale * G_get_raster_sample(fdrast,
- &window,
- NULL,
- Points->
- y[j],
- Points->
- x[j], 0,
- method);
-
- /* update the elevation value for each data point */
- Points->z[j] = estimated_elevation;
- } /* end looping through point in a line */
- break;
-
- /* lines with at least 3 vertexes */
- case GV_FACE:
- if (Points->n_points < 3)
- break; /* At least 3 points */
-
- /* loop through each point in a line */
- for (j = 0; j < Points->n_points; j++) {
- /* sample raster at this point, and update the z-coordinate (note that input vector should not be 3D!) */
- estimated_elevation = scale * G_get_raster_sample(fdrast,
- &window,
- NULL,
- Points->
- y[j],
- Points->
- x[j], 0,
- method);
-
- /* update the elevation value for each data point */
- Points->z[j] = estimated_elevation;
+ /* write the new line file, with the updated Points struct */
+ if (sample_raster
+ (ltype, fdrast, window, Points, method, scale, null_opt,
+ null_val)) {
+
+ /* Open new file only if there is something to write in */
+ if (out_num < 1) {
+ if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to create vector map <%s>"),
+ out_opt->answer);
+ }
+ Vect_copy_head_data(&In, &Out);
+ Vect_hist_copy(&In, &Out);
+ Vect_hist_command(&Out);
+ }
+
+ Vect_write_line(&Out, ltype, Points, Cats);
+ out_num++;
}
- break;
- } /* end line type switch */
+ } /* end looping thru lines */
+ }
- /* write the new line file, with the updated Points struct */
- Vect_write_line(&Out, ltype, Points, Cats);
- } /* end looping thru lines */
-
} /* end working on type=lines */
/* close elevation raster: */
G_close_cell(fdrast);
-
+
+ /* build topology for output vector */
+ if (out_num > 0) {
+ Vect_build(&Out);
+
+ /* Now let's move attribute data from all old map layers to new map */
+ for (i = 0; i < Out.plus.n_cidx; i++) {
+ new_cats = (int *)G_malloc(Out.plus.cidx[i].n_cats * sizeof(int));
+ if (!new_cats) {
+ G_warning(_("Due to error attribute data to new map are not transferred"));
+ break;
+ }
+ /* Vect_copy_table_by_cats does not accept Map.plus.cidx[].cat array as input.
+ Thus we construct new array of cats. */
+ for (c = 0; c < Out.plus.cidx[i].n_cats; c++) {
+ new_cats[c] = Out.plus.cidx[i].cat[c][0];
+ }
+ if (0 > Vect_copy_table_by_cats(&In, &Out, In.plus.cidx[i].field, Out.plus.cidx[i].field,
+ NULL, otype, new_cats, Out.plus.cidx[i].n_cats))
+ G_warning(_("Due to error attribute data to new map are not transferred"));
+ }
+ /* close output vector */
+ Vect_close(&Out);
+ }
+ else {
+ /* close input vector */
+ Vect_close(&In);
+ G_warning(_("No features drapped. Check Your computational region and input raster map."));
+ exit(EXIT_FAILURE);
+ }
/* close input vector */
Vect_close(&In);
- /* build topology for output vector */
- Vect_build(&Out);
- /* close output vector */
- Vect_close(&Out);
-
+
exit(EXIT_SUCCESS);
}
Modified: grass/trunk/vector/v.drape/v.drape.html
===================================================================
--- grass/trunk/vector/v.drape/v.drape.html 2008-11-27 10:18:44 UTC (rev 34527)
+++ grass/trunk/vector/v.drape/v.drape.html 2008-11-27 11:16:49 UTC (rev 34528)
@@ -1,46 +1,46 @@
<h2>DESCRIPTION</h2>
-<em>v.drape</em> converts 2D/3D vector data into 3D vector format via
+<p><em>v.drape</em> converts 2D/3D vector data into 3D vector format via
sampling of an elevation surface. Three sampling algorithms adapted
from <a href="v.sample.html">v.sample</a> were incorporated into this
-module: nearest neighbor, bilinear, and cubic convultion.
+module: nearest neighbor, bilinear, and cubic convultion.</p>
-<h2>NOTES</h2>
+<p>v.drape will skip vector features outside of current computational region or
+where raster map has NULL value. It's possible to include all vector features
+by specifying height value that will be assigned to verticles whose values
+can not be determined from raster map.</p>
-Please run beforehand
-
-<div class="code"><pre>
-g.region vect=2D_vector
-</pre></div>
-
-and make sure that the extent of the elevation raster is at least as
-big as the vector to convert.
-
-<P>
+<h2>NOTES</h2>
+<p>
Additional vertices can be added to the input 2D vector map
-with <a href="v.split.html">v.split</a>.
+with <a href="v.split.html">v.split</a>.</p>
-<P>
+<p>
The module can be used in conjunction
with <a href="v.out.pov.html">v.out.pov</a> and
<a href="r.out.pov.html">r.out.pov</a> to export a complete set of
-vector and raster data for display in POVRAY.
+vector and raster data for display in POVRAY.</p>
-<h2>EXAMPLE</h2>
+<h2>EXAMPLES</h2>
-Spearfish example:
-
+<p>Spearfish example:
<div class="code"><pre>
-g.region -p vect=roads align=elevation.10m
v.drape in=roads rast=elevation.10m method=bilinear out=roads3d
v.info roads3d
</pre></div>
+</p>
+<p>Create 3D vector roads map containing only "unimproved" roads.
+Set road height to 1000 m for all parts without height information.
+<div class="code"><pre>
+v.drape input=roads type=line rast=elevation.dem output=roads_3d \
+method=nearest scale=1.0 where='cat=5' layer=1 null_value=1000
+</pre></div>
+</p>
+
<h2>POVRAY EXAMPLE</h2>
<div class="code"><pre>
-#setup the region
-g.region vect=roads align=elevation.10m -p
#export the vector data
v.drape in=roads out=roads3d rast=elevation.10m
v.out.pov roads3d out=roads3d.pov
@@ -51,54 +51,17 @@
# now write a complete povray-script and launch povray
</pre></div>
-<h2>ERROR MESSAGES</h2>
-
-If the following error message appears
-
-<div class="code"><pre>
-WARNING: Current region does not include the entire input vector map
-WARNING: Reading raster map request for row 466 is outside region
-ERROR: Problem reading raster map
-</pre></div>
-
-it indicates that the vector map is spatially larger than the current
-region settings. To avoid this problem, set region from the input
-vector map and re-run <em>v.drape</em>.
-
-<div class="code"><pre>
-g.region vect=vectmap
-</pre></div>
-
-If the following error message appears
-
-<div class="code"><pre>
-WARNING: Elevation raster map does not cover the entire area of the input vector map.
-</pre></div>
-
-it indicates that the vector map is spatially larger than the raster map.
-To avoid this problem, the vector map needs to be clipped to the raster
-map extent, for example:
-
-<div class="code"><pre>
-g.region rast=demname
-v.in.region clipbox
-v.overlay ain=clipbox bin=vectmap out=vectmap_clipped op=and
-v.drape vectmap_clipped out=vectdrape rast=demname
-</pre></div>
-
-Then <em>v.drape</em> should perform the draping.
-
<h2>SEE ALSO</h2>
-<EM>
-<A HREF="r.out.pov.html">r.out.pov</A>,
-<A HREF="v.in.region.html">v.in.region</A>,
-<A HREF="v.out.pov.html">v.out.pov</A>,
-<A HREF="v.overlay.html">v.overlay</A>,
-<A HREF="v.split.html">v.split</A>,
-<A HREF="v.what.rast.html">v.what.rast</A>,
-<A HREF="v.extrude.html">v.extrude</A>
-</EM>
+<em>
+<a href="r.out.pov.html">r.out.pov</a>,
+<a href="v.in.region.html">v.in.region</a>,
+<a href="v.out.pov.html">v.out.pov</a>,
+<a href="v.overlay.html">v.overlay</a>,
+<a href="v.split.html">v.split</a>,
+<a href="v.what.rast.html">v.what.rast</a>,
+<a href="v.extrude.html">v.extrude</a>
+</em>
<h2>AUTHOR</h2>
More information about the grass-commit
mailing list