v.drape: major rewrite for G7

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

vector/v.drape/local_proto.h

@@ -0,0 +1,2 @@
+int sample_raster(int, const struct Cell_head *,
+		  struct line_pnts *, INTERP_TYPE, double, int, double);

vector/v.drape/main.c

@@ -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;
-
-    /* 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;
-}
+#include "local_proto.h"
 
 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;
-
-    int line, nlines, otype, ltype, ncats =
-	0, layer, i, c, *cats, field_index, id, out_num = 0, *new_cats;
+    struct cat_list *cat_list;
 
+    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 @@ int main(int argc, char *argv[])
     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.");
-
-    in_opt = G_define_standard_option(G_OPT_V_INPUT);
+	_("Converts 2D vector features to 3D by sampling of elevation raster map.");
 
-    layer_opt = G_define_standard_option(G_OPT_V_FIELD_ALL);
+    opt.input = G_define_standard_option(G_OPT_V_INPUT);
 
-    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");
+
+    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");
+
+    opt.output = G_define_standard_option(G_OPT_V_OUTPUT);
 
-    where_opt = G_define_standard_option(G_OPT_DB_WHERE);
+    opt.rast = G_define_standard_option(G_OPT_R_ELEV);
+    opt.rast->description = _("Elevation raster map for height extraction");
 
-    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.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);
-
-    /* Check output type */
-    otype = Vect_option_to_types(type_opt);
-
-    /* open the elev raster, and check for error condition */
-    fdrast = Rast_open_old(rast_opt->answer, "");
+    scale = atof(opt.scale->answer);
+
+    /* is null value defined */
+    if (opt.null->answer)
+	null_val = atof(opt.null->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, &ltype, &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, &ltype, &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);
-    
-    /* build topology for output vector */
-    if (out_num > 0) {
-	Vect_build(&Out);
+    /* 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;
         
-        /* 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"));
+        /* 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);
         }
-
-	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);
     }
+    
+    /* 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 */
+    Vect_build(&Out);
+        
+    Vect_get_map_box(&Out, &map_box);
 
-    /* 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);
     

vector/v.drape/sample.c

@@ -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;
+}

vector/v.drape/v.drape.html

@@ -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 @@ r.out.png landcover.30m out=landcover30m.png
 <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>