On Monday 09 July 2007 20:07, Hamish wrote:
> Dylan Beaudette wrote:
> > The v.drape module was created by suggestion of Radim, and implemented
> > by an amateur (me). This bug is one of the last remaining issues in
> > the code. Although it is mentioned in the manual page, it would be
> > good to do a simple test of
> >
> > if vector %in% current region or vector %in% raster bbox
> > then run else error and exit
> >
> > I am not sure how to implement this check, but if I get some feedback
> > I will submit a patch!
>
> I think the libgis error is appropriate, but that v.drape should not be
> trying to load a raster array i,j which is out of range. Perhaps in
> these cases the Z value should be set to NULL (nan), or if that is
> illegal either set to 0 (with warning) or G_fatal_error() out.
>
> You can check in map is contained within the current region by reading
> the map header info and the current region settings,
>
> e.g. read current region into "window":
>
> struct Cell_head window;
> G_get_window(&window);
>
> then you have doubles: window.north, .south, .east, .west, .top, .bottom
>
>
>
> To read vector bounds:
>
> struct Map_info Map;
> BOUND_BOX box;
>
> Vect_open_old_head (&Map, input->answer, mapset);
> Vect_get_map_box (&Map, &box);
>
> then you have doubles: box.N, box.S, box.E, box.W, box.T, box.B
>
>
>
> compare:
>
> if ( (box.N > window.north) || (box.S < window.south) ||
> (box.E > window.east) || (box.W < window.west) )
> G_warning("Vector exists outside of raster region");
>
> ...
>
>
> Hamish
A patch based on this suggestion is attached. I have tested and it appears to
work.
Dylan
Index: vector/v.drape/description.html
===================================================================
RCS file: /home/grass/grassrepository/grass6/vector/v.drape/description.html,v
retrieving revision 1.4
diff -u -r1.4 description.html
--- vector/v.drape/description.html 9 Jul 2007 15:15:36 -0000 1.4
+++ vector/v.drape/description.html 10 Jul 2007 06:10:26 -0000
@@ -49,16 +49,16 @@
<h2>ERROR MESSAGES</h2>
-If the following error message appears
+If the following error message appears:
<div class="code"><pre>
- WARNING: [demname in mapset] - read request for row -1 is outside region
- ERROR: problem reading raster cell file
+ERROR: Vector map [input_map] extends beyond the
+ current region.
</pre></div>
-it indicates that the vector map is spatially larger than the raster map.
+it indicates that the vector map has a bounding box larger than that of the current region.
To avoid this problem, the vector map needs to be clipped to the raster
-map extent, for example:
+map extent or to the curent region, for example:
<div class="code"><pre>
g.region rast=demname
Index: vector/v.drape/main.c
===================================================================
RCS file: /home/grass/grassrepository/grass6/vector/v.drape/main.c,v
retrieving revision 1.11
diff -u -r1.11 main.c
--- vector/v.drape/main.c 15 Apr 2007 23:16:22 -0000 1.11
+++ vector/v.drape/main.c 10 Jul 2007 06:10:26 -0000
@@ -26,9 +26,6 @@
\date 2005.09.20
\todo add support for areas
-\todo Enhanced error checking such as
- \li does the elevation raster cover the entire are of the vector map?
- \li does the current region include the entire input vector map ?
\todo Make a description.html for documentation
*/
@@ -44,111 +41,137 @@
int main(int argc, char *argv[])
{
- struct GModule *module;
- struct Option *in_opt, *out_opt, *type_opt, *rast_opt, *method_opt;
- struct Map_info In, Out;
- struct line_pnts *Points;
- struct line_cats *Cats;
- /* int layer; */
- int line, nlines, otype, ltype;
-
- /* Raster stuff from v.sample::main.c */
- char *mapset;
- int j;
- double scale, estimated_elevation;
- INTERP_TYPE method = UNKNOWN;
- int fdrast; /* file descriptor for raster map is int */
- struct Cell_head window;
- /* end raster stuff */
-
- G_gisinit(argv[0]);
-
- module = G_define_module();
- module->keywords = _("vector");
- module->description =
- _("Convert 2D vector to 3D vector by sampling of elevation raster. Default sampling by nearest neighbor");
-
- in_opt = G_define_standard_option(G_OPT_V_INPUT);
-
- type_opt = G_define_standard_option(G_OPT_V_TYPE);
- type_opt->options = "point,centroid,line,boundary,face,kernel";
- type_opt->answer = "point,centroid,line,boundary,face,kernel";
-
- /* raster sampling */
- rast_opt = G_define_option();
- rast_opt->key = "rast";
- rast_opt->type = TYPE_STRING;
- rast_opt->required = NO;
- rast_opt->description = _("Elevation raster for height extraction");
-
- method_opt = G_define_option();
- method_opt->key = "method";
- method_opt->type = TYPE_STRING;
- method_opt->required = NO;
- method_opt->multiple = NO;
- method_opt->options = "nearest,bilinear,cubic";
- method_opt->answer = "nearest";
- method_opt->descriptions = "nearest;nearest neighbor;"
+ struct GModule *module;
+ struct Option *in_opt, *out_opt, *type_opt, *rast_opt, *method_opt;
+ struct Map_info In, Out;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ /* int layer; */
+ int line, nlines, otype, ltype;
+ BOUND_BOX box;
+
+ /* Raster stuff from v.sample::main.c */
+ char *mapset;
+ int j;
+ double scale, estimated_elevation;
+ INTERP_TYPE method = UNKNOWN;
+ int fdrast; /* file descriptor for raster map is int */
+ struct Cell_head window;
+ /* end raster stuff */
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("vector");
+ module->description =
+ _("Convert 2D vector to 3D vector by sampling of elevation raster. Default sampling by nearest neighbor");
+
+ in_opt = G_define_standard_option(G_OPT_V_INPUT);
+
+ type_opt = G_define_standard_option(G_OPT_V_TYPE);
+ type_opt->options = "point,centroid,line,boundary,face,kernel";
+ type_opt->answer = "point,centroid,line,boundary,face,kernel";
+
+ /* raster sampling */
+ rast_opt = G_define_option();
+ rast_opt->key = "rast";
+ rast_opt->type = TYPE_STRING;
+ rast_opt->required = NO;
+ rast_opt->description = _("Elevation raster for height extraction");
+
+ method_opt = G_define_option();
+ method_opt->key = "method";
+ method_opt->type = TYPE_STRING;
+ method_opt->required = NO;
+ method_opt->multiple = NO;
+ method_opt->options = "nearest,bilinear,cubic";
+ method_opt->answer = "nearest";
+ method_opt->descriptions = "nearest;nearest neighbor;"
"bilinear;bilinear interpolation;"
"cubic;cubic convolution interpolation;";
- method_opt->description = _("Sampling method");
+ method_opt->description = _("Sampling method");
- out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* which interpolation method should we use */
+ if ( method_opt->answer[0] == 'b' )
+ method = BILINEAR;
+ else if ( method_opt->answer[0] == 'c' )
+ method = CUBIC;
+ else
+ {
+ /* G_message(_("Defaulting to nearest neighbor sampling.")); */
+ method = NEAREST;
+ }
- /* which interpolation method should we use */
- if ( method_opt->answer[0] == 'b' )
- method = BILINEAR;
- else if ( method_opt->answer[0] == 'c' )
- method = CUBIC;
- else
- {
- G_message(_("defaulting to nearest neighbor sampling"));
- method = NEAREST;
- }
+ /* setup the raster for sampling */
- /* setup the raster for sampling */
+ /* get the region parameters*/
+ G_get_window(&window);
- /* setup the region */
- G_get_window(&window);
-
- /* check for the elev raster, and check for error condition */
- if ((mapset = G_find_cell(rast_opt->answer, "")) == NULL) {
- G_fatal_error(_("cell file [%s] not found"), rast_opt->answer);
- }
+ /* check for the elev raster, and check for error condition */
+ if ((mapset = G_find_cell(rast_opt->answer, "")) == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), rast_opt->answer);
+ }
/* open the elev raster, and check for error condition */
if ((fdrast = G_open_cell_old(rast_opt->answer, mapset)) < 0) {
G_fatal_error(_("can't open raster map [%s]"), rast_opt->answer);
}
- /* used to scale sampled raster values: will need to add an option to modify this later */
- scale = 1;
+ /* used to scale sampled raster values: will need to add an option to modify this later */
+ scale = 1;
+
+
+ /* Check output type */
+ otype = Vect_option_to_types(type_opt);
+
+ Vect_set_open_level(2);
+
+ /* open the input vector */
+ Vect_open_old(&In, in_opt->answer, mapset);
+
+ /* get the input vector header */
+ Vect_open_old_head(&In, in_opt->answer, mapset);
+ Vect_get_map_box(&In, &box);
+
+ /*
+ * now need to check to make sure that vector is contained within the current region bounds
+ * if not, then issue a warning and optionally exit
+ */
+ if ( (box.N > window.north) || (box.S < window.south) || (box.E > window.east) || (box.W < window.west) )
+ {
+ /* clean-up */
+ /* close elevation raster: */
+ G_close_cell(fdrast);
+
+ /* close input vector */
+ Vect_close(&In);
+
+ /* issue error message */
+ G_fatal_error(_("Vector map [%s] extends beyond the current region."), in_opt->answer);
+ }
- /* Check output type */
- otype = Vect_option_to_types(type_opt);
- Vect_set_open_level(2);
- Vect_open_old(&In, in_opt->answer, "");
-
- /* setup the new vector map */
- /* remember to open the new vector map as 3D: Vect_open_new(,,1) */
- Vect_open_new(&Out, out_opt->answer, 1);
- 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);
+ /* setup the new vector map */
+ /* remember to open the new vector map as 3D: Vect_open_new(,,1) */
+ Vect_open_new(&Out, out_opt->answer, 1);
+ 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();
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
- /* line types */
- if ((otype &
+ /* line types */
+ if ((otype &
(GV_POINTS | GV_LINES | GV_BOUNDARY | GV_CENTROID | GV_FACE |
GV_KERNEL))) {
@@ -218,15 +241,15 @@
} /* end working on type=lines */
- /* close elevation raster: */
- G_close_cell(fdrast);
+ /* close elevation raster: */
+ G_close_cell(fdrast);
- /* close input vector */
- Vect_close(&In);
- /* build topology for output vector */
- Vect_build(&Out, stderr);
- /* close output vector */
- Vect_close(&Out);
+ /* close input vector */
+ Vect_close(&In);
+ /* build topology for output vector */
+ Vect_build(&Out, stderr);
+ /* close output vector */
+ Vect_close(&Out);
- exit(EXIT_SUCCESS);
+ exit(EXIT_SUCCESS);
}
_______________________________________________
grassuser mailing list
[email protected]
http://grass.itc.it/mailman/listinfo/grassuser
