[GRASS-SVN] r56262 - grass/trunk/vector/v.drape
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 15 06:29:39 PDT 2013
Author: martinl
Date: 2013-05-15 06:29:38 -0700 (Wed, 15 May 2013)
New Revision: 56262
Added:
grass/trunk/vector/v.drape/local_proto.h
grass/trunk/vector/v.drape/sample.c
Modified:
grass/trunk/vector/v.drape/main.c
grass/trunk/vector/v.drape/v.drape.html
Log:
v.drape: major rewrite for G7
Added: grass/trunk/vector/v.drape/local_proto.h
===================================================================
--- grass/trunk/vector/v.drape/local_proto.h (rev 0)
+++ grass/trunk/vector/v.drape/local_proto.h 2013-05-15 13:29:38 UTC (rev 56262)
@@ -0,0 +1,2 @@
+int sample_raster(int, const struct Cell_head *,
+ struct line_pnts *, INTERP_TYPE, double, int, double);
Property changes on: grass/trunk/vector/v.drape/local_proto.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass/trunk/vector/v.drape/main.c
===================================================================
--- grass/trunk/vector/v.drape/main.c 2013-05-15 12:45:35 UTC (rev 56261)
+++ grass/trunk/vector/v.drape/main.c 2013-05-15 13:29:38 UTC (rev 56262)
@@ -5,11 +5,12 @@
*
* AUTHOR(S): Radim Blazek, Dylan Beaudette
* Maris Nartiss - WHERE support and raster NULL support
- * OGR support by Martin Landa <landa.martin gmail.com>
+ * OGR support & major rewrite for GRASS 7 by Martin
+ * Landa <landa.martin gmail.com>
*
* PURPOSE: Convert 2D vector to 3D vector by sampling of elevation raster.
*
- * COPYRIGHT: (C) 2005-2009 by the GRASS Development Team
+ * COPYRIGHT: (C) 2005-2009, 2013 by the GRASS Development Team
*
* This program is free software under the GNU General
* Public License (>=v2). Read the file COPYING that
@@ -25,115 +26,30 @@
#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;
+#include "local_proto.h"
- /* 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 * Rast_get_sample(
- fdrast, &window, NULL, Points->y[0], Points->x[0], 0, method);
- /* Elevation value has to be meaningfull */
- if (Rast_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 * Rast_get_sample(
- fdrast, &window, NULL, Points->y[j], Points->x[j], 0, method);
-
- if (Rast_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 * Rast_get_sample(
- fdrast, &window, NULL, Points->y[j], Points->x[j], 0, method);
-
- if (Rast_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, *where_opt, *layer_opt, *null_opt;
+ struct {
+ struct Option *input, *output, *type, *rast, *method,
+ *scale, *where, *layer, *null, *cats;
+ } opt;
struct Map_info In, Out;
struct line_pnts *Points;
struct line_cats *Cats;
- struct field_info *Fi;
+ struct cat_list *cat_list;
- int line, nlines, otype, ltype, ncats =
- 0, layer, i, c, *cats, field_index, id, out_num = 0, *new_cats;
-
+ int otype, layer;
double scale, null_val;
- INTERP_TYPE method = UNKNOWN;
+ INTERP_TYPE method;
int fdrast; /* file descriptor for raster map is int */
+ int nlines, line, ltype;
+
struct Cell_head window;
struct bound_box map_box;
- dbDriver *driver;
- dbHandle handle;
- dbTable *table;
- dbString table_name, dbsql, valstr;
-
G_gisinit(argv[0]);
module = G_define_module();
@@ -141,238 +57,162 @@
G_add_keyword(_("geometry"));
G_add_keyword(_("sampling"));
G_add_keyword(_("3D"));
-
module->description =
- _("Converts vector map to 3D by sampling of elevation raster map.");
+ _("Converts 2D vector features to 3D by sampling of elevation raster map.");
- in_opt = G_define_standard_option(G_OPT_V_INPUT);
+ opt.input = G_define_standard_option(G_OPT_V_INPUT);
- layer_opt = G_define_standard_option(G_OPT_V_FIELD_ALL);
-
- type_opt = G_define_standard_option(G_OPT_V3_TYPE);
+ opt.layer = G_define_standard_option(G_OPT_V_FIELD_ALL);
+ opt.layer->guisection = _("Selection");
- /* raster sampling */
- rast_opt = G_define_standard_option(G_OPT_R_MAP);
- rast_opt->key = "rast";
- rast_opt->description = _("Elevation raster map for height extraction");
-
- out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.cats = G_define_standard_option(G_OPT_V_CATS);
+ opt.cats->guisection = _("Selection");
- method_opt = G_define_standard_option(G_OPT_R_INTERP_TYPE);
- method_opt->answer = "nearest";
- 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";
+ opt.where = G_define_standard_option(G_OPT_DB_WHERE);
+ opt.where->guisection = _("Selection");
- where_opt = G_define_standard_option(G_OPT_DB_WHERE);
+ opt.type = G_define_standard_option(G_OPT_V_TYPE);
+ opt.type->options = "point,line,boundary,centroid";
+ opt.type->answer = "point,line,boundary,centroid";
+ opt.type->guisection = _("Selection");
- 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");
+ opt.output = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.rast = G_define_standard_option(G_OPT_R_ELEV);
+ opt.rast->description = _("Elevation raster map for height extraction");
+
+ opt.method = G_define_standard_option(G_OPT_R_INTERP_TYPE);
+ opt.method->answer = "nearest";
+ opt.method->guisection = _("Elevation");
+
+ opt.scale = G_define_option();
+ opt.scale->key = "scale";
+ opt.scale->type = TYPE_DOUBLE;
+ opt.scale->description = _("Scale factor sampled raster values");
+ opt.scale->answer = "1.0";
+ opt.scale->guisection = _("Elevation");
+
+ opt.null = G_define_option();
+ opt.null->key = "null_value";
+ opt.null->type = TYPE_DOUBLE;
+ opt.null->description =
+ _("Height for sampled raster NULL values");
+ opt.null->guisection = _("Elevation");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- if (null_opt->answer)
- null_val = atof(null_opt->answer);
-
/* which interpolation method should we use */
- method = Rast_option_to_interp_type(method_opt);
+ method = Rast_option_to_interp_type(opt.method);
/* setup the region */
G_get_window(&window);
/* used to scale sampled raster values */
- scale = atof(scale_opt->answer);
+ scale = atof(opt.scale->answer);
- /* Check output type */
- otype = Vect_option_to_types(type_opt);
+ /* is null value defined */
+ if (opt.null->answer)
+ null_val = atof(opt.null->answer);
- /* open the elev raster, and check for error condition */
- fdrast = Rast_open_old(rast_opt->answer, "");
+ /* check output type */
+ otype = Vect_option_to_types(opt.type);
+ if (otype & GV_AREA) {
+ if (otype & GV_BOUNDARY) {
+ /* process area -> skip boundaries */
+ otype &= ~GV_BOUNDARY;
+ }
+ if (otype & GV_CENTROID) {
+ /* process area -> skip centroids */
+ otype &= ~GV_CENTROID;
+ }
+ }
+ /* open the elev raster map */
+ fdrast = Rast_open_old(opt.rast->answer, "");
+
/* check input/output vector maps */
- Vect_check_input_output_name(in_opt->answer, out_opt->answer,
+ Vect_check_input_output_name(opt.input->answer, opt.output->answer,
G_FATAL_EXIT);
- Vect_set_open_level(2);
- Vect_open_old2(&In, in_opt->answer, "", layer_opt->answer);
- layer = Vect_get_field_number(&In, layer_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);
+ /* open input vector map */
+ Vect_set_open_level(2); /* topology required ? */
+ Vect_open_old2(&In, opt.input->answer, "", opt.layer->answer);
+ Vect_set_error_handler_io(&In, &Out);
+
+ /* get layer number */
+ layer = Vect_get_field_number(&In, opt.layer->answer);
+ if ((opt.cats->answer || opt.where->answer) && layer == -1) {
+ G_warning(_("Invalid layer number (%d). "
+ "Parameter '%s' or '%s' specified, assuming layer '1'."),
+ layer, opt.cats->key, opt.where->key);
+ layer = 1;
}
+ /* create output */
+ Vect_open_new(&Out, opt.output->answer, WITH_Z);
+
+ Vect_copy_head_data(&In, &Out);
+ Vect_hist_copy(&In, &Out);
+ Vect_hist_command(&Out);
+
+ /* set constraint for cats or where options */
+ cat_list = NULL;
+ if (layer > 0)
+ cat_list = Vect_cats_set_constraint(&In, layer, opt.where->answer,
+ opt.cats->answer);
+
+ /* allocate space for points and cats */
Points = Vect_new_line_struct();
Cats = Vect_new_cats_struct();
- /* line types */
- if ((otype &
- (GV_POINTS | GV_LINES | GV_BOUNDARY | GV_CENTROID | GV_FACE |
- GV_KERNEL))) {
-
- /* 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++;
- }
- }
- 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++) {
-
- /* progress feedback */
- G_percent(line, nlines, 2);
-
- /* get the line type */
- ltype = Vect_read_line(&In, Points, Cats, line);
- if (layer != -1 && !Vect_cat_get(Cats, layer, NULL))
- continue;
-
- /* 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++;
- }
- } /* end looping thru lines */
- }
-
- } /* end working on type=lines */
-
- /* close elevation raster: */
- Rast_close(fdrast);
+ /* loop through each line */
+ nlines = Vect_get_num_lines(&In);
+ G_important_message(_("Processing features..."));
+ for (line = 1; line <= nlines; line++) {
+ /* progress feedback */
+ G_percent(line, nlines, 2);
+
+ if (!Vect_line_alive(&In, line))
+ continue;
+
+ /* get the line type */
+ ltype = Vect_read_line(&In, Points, Cats, line);
+ if (!(ltype & otype))
+ continue;
+ if (layer > 0 && !Vect_cats_in_constraint(Cats, layer, cat_list))
+ continue;
+
+ /* write the new line file, with the updated Points struct */
+ if (sample_raster(fdrast, &window, Points,
+ method, scale,
+ opt.null->answer ? TRUE : FALSE, null_val)) {
+ Vect_write_line(&Out, ltype, Points, Cats);
+ }
+ else {
+ G_warning(_("Undefined height for feature %d. Skipping."), line);
+ }
+ }
+ /* copy attribute data */
+ G_important_message(_("Copying attribute tables..."));
+ if (layer < 0)
+ Vect_copy_tables(&In, &Out, 0);
+ else
+ Vect_copy_table_by_cat_list(&In, &Out, layer, layer, NULL,
+ GV_1TABLE, cat_list);
+
/* build topology for output vector */
- if (out_num > 0) {
- Vect_build(&Out);
+ 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"));
- }
+ Vect_get_map_box(&Out, &map_box);
- Vect_get_map_box(&Out, &map_box);
-
- /* close output vector */
- Vect_close(&Out);
- }
- else {
- /* close input vector */
- Vect_close(&In);
- G_warning(_("No features drapped. Check your computational region and input vector map."));
- exit(EXIT_SUCCESS);
- }
-
- /* close input vector */
+ /* close elevation raster map */
+ Rast_close(fdrast);
+ /* close input vector map */
Vect_close(&In);
+ /* close output vector map */
+ Vect_close(&Out);
G_done_msg("T: %f B: %f.", map_box.T, map_box.B);
Added: grass/trunk/vector/v.drape/sample.c
===================================================================
--- grass/trunk/vector/v.drape/sample.c (rev 0)
+++ grass/trunk/vector/v.drape/sample.c 2013-05-15 13:29:38 UTC (rev 56262)
@@ -0,0 +1,35 @@
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+/* Samples raster map */
+int sample_raster(int fdrast, const struct Cell_head *window,
+ struct line_pnts *Points,
+ INTERP_TYPE method, double scale,
+ int null_defined, double null_val)
+{
+ double estimated_elevation;
+ int j;
+
+ /* 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 * Rast_get_sample(fdrast, window, NULL,
+ Points->y[j], Points->x[j],
+ 0, method);
+
+ if (Rast_is_d_null_value(&estimated_elevation)) {
+ if (null_defined)
+ estimated_elevation = null_val;
+ else
+ return 0;
+ }
+
+ /* update the elevation value for each data point */
+ Points->z[j] = estimated_elevation;
+ }
+
+ return 1;
+}
Property changes on: grass/trunk/vector/v.drape/sample.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass/trunk/vector/v.drape/v.drape.html
===================================================================
--- grass/trunk/vector/v.drape/v.drape.html 2013-05-15 12:45:35 UTC (rev 56261)
+++ grass/trunk/vector/v.drape/v.drape.html 2013-05-15 13:29:38 UTC (rev 56262)
@@ -1,40 +1,48 @@
<h2>DESCRIPTION</h2>
-<p><em>v.drape</em> converts 2D/3D vector data into 3D vector format via
+<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.
+from <em><a href="v.sample.html">v.sample</a></em> were incorporated
+into this module: nearest neighbor, bilinear, and cubic convultion.
-<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>
+<em>v.drape</em> 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.
<h2>NOTES</h2>
-<p>Additional vertices can be added to the input 2D vector map
-with <a href="v.split.html">v.split</a>.
-<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.
+Additional vertices can be added to the input 2D vector map
+with <em><a href="v.split.html">v.split</a></em>.
+<p>
+The module can be used in conjunction
+with <em><a href="v.out.pov.html">v.out.pov</a></em> and
+<em><a href="r.out.pov.html">r.out.pov</a></em> to export a complete
+set of vector and raster data for display
+in <a href="http://www.povray.org/">POVRAY</a>.
+
<h2>EXAMPLES</h2>
-<p>Spearfish example:
+Spearfish example:
+
<div class="code"><pre>
v.drape in=roads rast=elevation.10m method=bilinear out=roads3d
v.info roads3d
</pre></div>
-<p>Create 3D vector roads map containing only "unimproved" roads.
-Set road height to 1000 m for all parts without height information.
+<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
+ method=nearest scale=1.0 where='cat=5' layer=1 null_value=1000
</pre></div>
-<h2>POVRAY EXAMPLE</h2>
+<h3>POVRAY example</h3>
<div class="code"><pre>
#export the vector data
@@ -50,17 +58,21 @@
<h2>SEE ALSO</h2>
<em>
+<a href="v.extrude.html">v.extrude</a>,
+<a href="v.to.3d.html">v.to.3d</a>,
<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>
+<a href="v.what.rast.html">v.what.rast</a>
</em>
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
-Dylan Beaudette, University of California at Davis.
+Dylan Beaudette, University of California at Davis.<br>
+Updated for GRASS 7 by Martin Landa, Czech Technical University in
+Prague, Czech Republic
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>
More information about the grass-commit
mailing list