[GRASS-SVN] r51618 - in grass/trunk/vector: . v.what.rast3
v.what.rast3/test_suite
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon May 14 04:08:54 EDT 2012
Author: huhabla
Date: 2012-05-14 01:08:53 -0700 (Mon, 14 May 2012)
New Revision: 51618
Added:
grass/trunk/vector/v.what.rast3/
grass/trunk/vector/v.what.rast3/Makefile
grass/trunk/vector/v.what.rast3/local_proto.h
grass/trunk/vector/v.what.rast3/main.c
grass/trunk/vector/v.what.rast3/search.c
grass/trunk/vector/v.what.rast3/test_suite/
grass/trunk/vector/v.what.rast3/test_suite/random_points.ref
grass/trunk/vector/v.what.rast3/test_suite/random_points.txt
grass/trunk/vector/v.what.rast3/test_suite/random_points_db.ref
grass/trunk/vector/v.what.rast3/test_suite/test.v.what.rast3.sh
grass/trunk/vector/v.what.rast3/v.what.rast3.html
Modified:
grass/trunk/vector/Makefile
Log:
Added a new module to sample 3d raster maps with vector points.
Modified: grass/trunk/vector/Makefile
===================================================================
--- grass/trunk/vector/Makefile 2012-05-14 07:57:51 UTC (rev 51617)
+++ grass/trunk/vector/Makefile 2012-05-14 08:08:53 UTC (rev 51618)
@@ -90,6 +90,7 @@
v.voronoi \
v.what \
v.what.rast \
+ v.what.rast3 \
v.vect.stats \
v.vol.rst \
v.out.ogr \
Added: grass/trunk/vector/v.what.rast3/Makefile
===================================================================
--- grass/trunk/vector/v.what.rast3/Makefile (rev 0)
+++ grass/trunk/vector/v.what.rast3/Makefile 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,14 @@
+
+MODULE_TOPDIR = ../..
+
+PGM=v.what.rast3
+
+LIBES = $(VECTORLIB) $(DBMILIB) $(RASTER3DLIB) $(GISLIB)
+DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(RASTER3DDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass/trunk/vector/v.what.rast3/local_proto.h
===================================================================
--- grass/trunk/vector/v.what.rast3/local_proto.h (rev 0)
+++ grass/trunk/vector/v.what.rast3/local_proto.h 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,16 @@
+#include <grass/gis.h>
+
+struct order
+{
+ int cat; /* point category */
+ int count; /* nuber of points with category 'cat' */
+ double x;
+ double y;
+ double z;
+ FCELL fvalue;
+ DCELL dvalue;
+};
+
+/* search.c */
+int by_cat(const void *, const void *);
+int srch_cat(const void *, const void *);
Added: grass/trunk/vector/v.what.rast3/main.c
===================================================================
--- grass/trunk/vector/v.what.rast3/main.c (rev 0)
+++ grass/trunk/vector/v.what.rast3/main.c 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,333 @@
+/* ***************************************************************
+ *
+ * MODULE: v.what.rast3
+ *
+ * AUTHOR(S): Soeren Gebbert, this code is based on a slightly modified version of
+ * v.what.rast from Radim Blazek and Michael Shapiro
+ *
+ * PURPOSE: Uploads 3d raster values at positions of vector points to the table
+ *
+ * COPYRIGHT: (C) 2001, 2011 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ * **************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <grass/raster3d.h>
+#include <grass/dbmi.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+int main(int argc, char *argv[])
+{
+ int i, j, nlines, type, field, cat;
+
+ char buf[2048];
+ struct {
+ struct Option *vect, *rast3d, *field, *col, *where;
+ } opt;
+ int Cache_size;
+ struct order *cache;
+ struct GModule *module;
+
+ RASTER3D_Region region;
+ RASTER3D_Map *map;
+
+ struct Map_info Map;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ int point;
+ int point_cnt; /* number of points in cache */
+ int outside_cnt; /* points outside region */
+ int nocat_cnt; /* points inside region but without category */
+ int dupl_cnt; /* duplicate categories */
+ int typeIntern;
+ int is_empty;
+ struct bound_box box;
+
+ int *catexst, *cex;
+ struct field_info *Fi;
+ dbString stmt;
+ dbDriver *driver;
+ int select, norec_cnt, update_cnt, upderr_cnt, col_type;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("vector"));
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("position"));
+ G_add_keyword(_("querying"));
+ G_add_keyword(_("attribute table"));
+ module->description =
+ _("Uploads 3d raster values at positions of vector points to the table.");
+
+ opt.vect = G_define_standard_option(G_OPT_V_MAP);
+ opt.vect->label =
+ _("Name of vector points map for which to edit attributes");
+
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
+
+ opt.rast3d = G_define_standard_option(G_OPT_R3_MAP);
+ opt.rast3d->key = "raster3d";
+ opt.rast3d->description = _("Name of existing 3d raster map to be queried");
+
+ opt.col = G_define_standard_option(G_OPT_DB_COLUMN);
+ opt.col->required = YES;
+ opt.col->description =
+ _("Name of attribute column to be updated with the query result");
+
+ opt.where = G_define_standard_option(G_OPT_DB_WHERE);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ db_init_string(&stmt);
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ /* Initiate the default settings */
+ Rast3d_init_defaults();
+
+ /* Figure out the current region settings */
+ Rast3d_get_window(®ion);
+
+ box.N = region.north;
+ box.S = region.south;
+ box.E = region.east;
+ box.W = region.west;
+ box.T = region.top;
+ box.B = region.bottom;
+
+ /* Open vector */
+ Vect_set_open_level(2);
+ Vect_open_old2(&Map, opt.vect->answer, "", opt.field->answer);
+
+ field = Vect_get_field_number(&Map, opt.field->answer);
+
+ Fi = Vect_get_field(&Map, field);
+ if (Fi == NULL)
+ G_fatal_error(_("Database connection not defined for layer %d"),
+ field);
+
+ /* Open driver */
+ driver = db_start_driver_open_database(Fi->driver, Fi->database);
+ if (driver == NULL) {
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+ }
+
+ map = Rast3d_open_cell_old(opt.rast3d->answer, G_find_raster3d(opt.rast3d->answer, ""), ®ion,
+ RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
+
+ if (map == NULL)
+ G_fatal_error(_("Error opening 3d raster map <%s>"), opt.rast3d->answer);
+
+ /* Check column type */
+ col_type = db_column_Ctype(driver, Fi->table, opt.col->answer);
+
+ if (col_type == -1)
+ G_fatal_error(_("Column <%s> not found"), opt.col->answer);
+
+ if (col_type != DB_C_TYPE_DOUBLE)
+ G_fatal_error(_("Column type not supported, please use a column with double type"));
+
+ /* Read vector points to cache */
+ Cache_size = Vect_get_num_primitives(&Map, GV_POINT);
+ /* Note: Some space may be wasted (outside region or no category) */
+
+ cache = (struct order *)G_calloc(Cache_size, sizeof(struct order));
+
+ point_cnt = outside_cnt = nocat_cnt = 0;
+
+ nlines = Vect_get_num_lines(&Map);
+
+ G_debug(1, "Reading %d vector features from map", nlines);
+
+ G_important_message(_("Reading features from vector map..."));
+ for (i = 1; i <= nlines; i++) {
+ type = Vect_read_line(&Map, Points, Cats, i);
+ G_debug(4, "line = %d type = %d", i, type);
+
+ G_percent(i, nlines, 2);
+
+ /* check type */
+ if (!(type & GV_POINT))
+ continue; /* Points only */
+
+ /* check region */
+ if (!Vect_point_in_box(Points->x[0], Points->y[0], Points->z[0], &box)) {
+ outside_cnt++;
+ continue;
+ }
+
+ Vect_cat_get(Cats, field, &cat);
+ if (cat < 0) { /* no category of given field */
+ nocat_cnt++;
+ continue;
+ }
+
+ G_debug(4, " cat = %d", cat);
+
+ cache[point_cnt].x = Points->x[0];
+ cache[point_cnt].y = Points->y[0];
+ cache[point_cnt].z = Points->z[0];
+ cache[point_cnt].cat = cat;
+ cache[point_cnt].count = 1;
+ point_cnt++;
+ }
+
+ Vect_set_db_updated(&Map);
+ Vect_hist_command(&Map);
+ Vect_set_db_updated(&Map);
+ Vect_close(&Map);
+
+ G_debug(1, "Read %d vector points", point_cnt);
+ /* Cache may contain duplicate categories, sort by cat, find and remove duplicates
+ * and recalc count and decrease point_cnt */
+ qsort(cache, point_cnt, sizeof(struct order), by_cat);
+
+ G_debug(1, "Points are sorted, starting duplicate removal loop");
+
+ for (i = 0, j = 1; j < point_cnt; j++)
+ if (cache[i].cat != cache[j].cat)
+ cache[++i] = cache[j];
+ else
+ cache[i].count++;
+ point_cnt = i + 1;
+
+ G_debug(1, "%d vector points left after removal of duplicates",
+ point_cnt);
+
+ /* Report number of points not used */
+ if (outside_cnt)
+ G_warning(_("%d points outside current region were skipped"),
+ outside_cnt);
+
+ if (nocat_cnt)
+ G_warning(_("%d points without category were skipped"), nocat_cnt);
+
+ /* Extract raster values from file and store in cache */
+ G_debug(1, "Extracting raster values");
+
+ typeIntern = Rast3d_tile_type_map(map);
+
+ /* Sample the 3d raster map */
+ for (point = 0; point < point_cnt; point++) {
+
+ if (cache[point].count > 1)
+ continue; /* duplicate cats */
+
+ if(typeIntern == FCELL_TYPE) {
+ Rast3d_get_window_value(map, cache[point].y, cache[point].x, cache[point].z,
+ &cache[point].fvalue, FCELL_TYPE);
+ } else {
+ Rast3d_get_window_value(map, cache[point].y, cache[point].x, cache[point].z,
+ &cache[point].dvalue, DCELL_TYPE);
+ }
+ }
+
+ Rast3d_close(map);
+
+ /* Update table from cache */
+ G_debug(1, "Updating db table");
+
+ /* select existing categories to array (array is sorted) */
+ select = db_select_int(driver, Fi->table, Fi->key, NULL, &catexst);
+
+ db_begin_transaction(driver);
+
+ norec_cnt = update_cnt = upderr_cnt = dupl_cnt = 0;
+
+ G_message("Update vector attributes...");
+ for (point = 0; point < point_cnt; point++) {
+ if (cache[point].count > 1) {
+ G_warning(_("More points (%d) of category %d, value set to 'NULL'"),
+ cache[point].count, cache[point].cat);
+ dupl_cnt++;
+ }
+
+ G_percent(point, point_cnt, 2);
+
+ /* category exist in DB ? */
+ cex =
+ (int *)bsearch((void *)&(cache[point].cat), catexst, select,
+ sizeof(int), srch_cat);
+ if (cex == NULL) { /* cat does not exist in DB */
+ norec_cnt++;
+ G_warning(_("No record for category %d in table <%s>"),
+ cache[point].cat, Fi->table);
+ continue;
+ }
+
+ G_snprintf(buf, 2048, "update %s set %s = ", Fi->table, opt.col->answer);
+
+ db_set_string(&stmt, buf);
+
+ is_empty = 0;
+
+ if (cache[point].count > 1)
+ is_empty = 1;
+ if(typeIntern == FCELL_TYPE)
+ is_empty = Rast3d_is_null_value_num(&cache[point].fvalue, FCELL_TYPE);
+ if(typeIntern == DCELL_TYPE)
+ is_empty = Rast3d_is_null_value_num(&cache[point].dvalue, DCELL_TYPE);
+
+ if (is_empty) {
+ G_snprintf(buf, 2048, "NULL");
+ }else {
+ if(typeIntern == FCELL_TYPE)
+ G_snprintf(buf, 2048, "%.10f", cache[point].fvalue);
+ if(typeIntern == DCELL_TYPE)
+ G_snprintf(buf, 2048, "%.15f", cache[point].dvalue);
+ }
+
+ db_append_string(&stmt, buf);
+
+ G_snprintf(buf, 2048, " where %s = %d", Fi->key, cache[point].cat);
+
+ db_append_string(&stmt, buf);
+ /* user provides where condition: */
+ if (opt.where->answer) {
+ G_snprintf(buf, 2048, " AND %s", opt.where->answer);
+ db_append_string(&stmt, buf);
+ }
+ G_debug(3, db_get_string(&stmt));
+
+ /* Update table */
+ if (db_execute_immediate(driver, &stmt) == DB_OK) {
+ update_cnt++;
+ }
+ else {
+ upderr_cnt++;
+ }
+ }
+ G_percent(1, 1, 1);
+
+ G_debug(1, "Committing DB transaction");
+ db_commit_transaction(driver);
+
+ G_free(catexst);
+ db_close_database_shutdown_driver(driver);
+ db_free_string(&stmt);
+
+ /* Report */
+ G_verbose_message(_("%d categories loaded from table"), select);
+ G_verbose_message(_("%d categories loaded from vector"), point_cnt);
+ G_verbose_message(_("%d categories from vector missing in table"), norec_cnt);
+ if (dupl_cnt > 0)
+ G_message(_("%d duplicate categories in vector"), dupl_cnt);
+ if (upderr_cnt > 0)
+ G_warning(_("%d update errors"), upderr_cnt);
+
+ G_done_msg(_("%d records updated."), update_cnt);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass/trunk/vector/v.what.rast3/search.c
===================================================================
--- grass/trunk/vector/v.what.rast3/search.c (rev 0)
+++ grass/trunk/vector/v.what.rast3/search.c 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,22 @@
+#include "local_proto.h"
+
+/* for qsort, order list by cat */
+int by_cat(const void *ii, const void *jj)
+{
+ const struct order *i = ii, *j = jj;
+
+ return i->cat - j->cat;
+}
+
+/* for bsearch, find cat */
+int srch_cat(const void *pa, const void *pb)
+{
+ int *p1 = (int *)pa;
+ int *p2 = (int *)pb;
+
+ if (*p1 < *p2)
+ return -1;
+ if (*p1 > *p2)
+ return 1;
+ return 0;
+}
Added: grass/trunk/vector/v.what.rast3/test_suite/random_points.ref
===================================================================
--- grass/trunk/vector/v.what.rast3/test_suite/random_points.ref (rev 0)
+++ grass/trunk/vector/v.what.rast3/test_suite/random_points.ref 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,40 @@
+ORGANIZATION:
+DIGIT DATE:
+DIGIT NAME: soeren
+MAP NAME:
+MAP DATE: Sun May 13 00:35:59 2012
+MAP SCALE: 1
+OTHER INFO:
+ZONE: 0
+MAP THRESH: 0.000000
+VERTI:
+P 1 1
+ 15.98122828 27.60680488 39.15496119
+ 1 1
+P 1 1
+ 20.15599665 63.81531506 9.87756846
+ 1 2
+P 1 1
+ 66.47772443 53.77607164 13.88873554
+ 1 3
+P 1 1
+ 44.60300442 33.41779363 31.44354624
+ 1 4
+P 1 1
+ 63.52155272 35.93806371 47.61148626
+ 1 5
+P 1 1
+ 8.3804932 44.49982096 35.86484647
+ 1 6
+P 1 1
+ 85.83974446 42.48782134 0.81502858
+ 1 7
+P 1 1
+ 75.71132294 9.60621038 40.20883771
+ 1 8
+P 1 1
+ 84.33209107 28.0661076 6.48952234
+ 1 9
+P 1 1
+ 89.1191198 69.92471626 10.91284527
+ 1 10
Added: grass/trunk/vector/v.what.rast3/test_suite/random_points.txt
===================================================================
--- grass/trunk/vector/v.what.rast3/test_suite/random_points.txt (rev 0)
+++ grass/trunk/vector/v.what.rast3/test_suite/random_points.txt 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,40 @@
+ORGANIZATION:
+DIGIT DATE:
+DIGIT NAME: soeren
+MAP NAME:
+MAP DATE: Sun May 13 00:35:59 2012
+MAP SCALE: 1
+OTHER INFO:
+ZONE: 0
+MAP THRESH: 0.000000
+VERTI:
+P 1 1
+ 15.98122828 27.60680488 39.15496119
+ 1 1
+P 1 1
+ 20.15599665 63.81531506 9.87756846
+ 1 2
+P 1 1
+ 66.47772443 53.77607164 13.88873554
+ 1 3
+P 1 1
+ 44.60300442 33.41779363 31.44354624
+ 1 4
+P 1 1
+ 63.52155272 35.93806371 47.61148626
+ 1 5
+P 1 1
+ 8.3804932 44.49982096 35.86484647
+ 1 6
+P 1 1
+ 85.83974446 42.48782134 0.81502858
+ 1 7
+P 1 1
+ 75.71132294 9.60621038 40.20883771
+ 1 8
+P 1 1
+ 84.33209107 28.0661076 6.48952234
+ 1 9
+P 1 1
+ 89.1191198 69.92471626 10.91284527
+ 1 10
Added: grass/trunk/vector/v.what.rast3/test_suite/random_points_db.ref
===================================================================
--- grass/trunk/vector/v.what.rast3/test_suite/random_points_db.ref (rev 0)
+++ grass/trunk/vector/v.what.rast3/test_suite/random_points_db.ref 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,11 @@
+cat|concentration
+1|11
+2|5
+3|11
+4|13
+5|16
+6|8
+7|13
+8|20
+9|15
+10|12
Added: grass/trunk/vector/v.what.rast3/test_suite/test.v.what.rast3.sh
===================================================================
--- grass/trunk/vector/v.what.rast3/test_suite/test.v.what.rast3.sh (rev 0)
+++ grass/trunk/vector/v.what.rast3/test_suite/test.v.what.rast3.sh 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,24 @@
+# This script tests v.what.rast3d
+
+# We specific a small region in the
+# @preprocess step
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=70 w=0 e=100 b=0 t=50 res=10 res3=10 -p3
+
+# Create the volume and the sampling vector map
+r3.mapcalc --o expr="plume = double(col() + row() + depth())"
+# This is how the input data was created
+# v.random --o -z seed=1 output=random_points n=10 zmin=0 zmax=50
+# v.out.ascii --o format=standard input=random_points output=random_points.txt
+
+v.in.ascii --o -z format=standard input=random_points.txt output=random_points
+v.db.addtable --o map=random_points column="concentration double precision"
+
+# @test the voxel sampling with vector points
+v.what.rast3 map=random_points raster3d=plume column=concentration
+
+# Some data export commands for reference data creation and visual validation
+#v.out.ascii --o format=standard input=random_points output=random_points.ref
+#v.db.select map=random_points > random_points_db.ref
+#r3.out.vtk --o input=plume output=plume.vtk null=0
+#v.out.vtk --o -n input=random_points output=random_points.vtk
Added: grass/trunk/vector/v.what.rast3/v.what.rast3.html
===================================================================
--- grass/trunk/vector/v.what.rast3/v.what.rast3.html (rev 0)
+++ grass/trunk/vector/v.what.rast3/v.what.rast3.html 2012-05-14 08:08:53 UTC (rev 51618)
@@ -0,0 +1,57 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.what.rast3</em> reads 3d raster value for each point in the vector and updates <b>col</b>
+column in vector attribute table by this value. The column should be type double. This module is based on
+<a href="v.what.rast.html">v.what.rast</a>.
+<br>
+If more points have the same category, attribute value is set to NULL.
+If 3d raster values is NULL, attribute value is set to NULL.
+
+<h2>NOTES</h2>
+
+<h2>EXAMPLES</h2>
+
+A) Reading values from 3d raster map at position of vector points, writing these values
+ into a column of the attribute table connected to the vector map:
+<br>
+<div class="code"><pre>
+v.what.rast3 map=pnts raster3d=plume column=concentration
+</pre></div>
+
+<p>
+B) In case of a vector map without attached attribute table, first add
+a new attribute table. This table is then populated with values
+queried from the raster map:
+<br>
+<div class="code"><pre>
+# create new random 3d vector points map
+v.random -z output=pnts n=100 zmin=0 zmax=50
+
+# add new table, link to map
+v.db.addtable map=pnts column="concentration double precision"
+
+# query raster map and upload values to vector table into specified column
+g.region rast3d=plume -p
+v.what.rast3 map=pnts raster3d=plume column=concentration
+
+# verify new attribute table:
+v.db.select map=pnts
+
+# verify statistics of uploaded values:
+v.univar map=pnts column=concentration type=point
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.what.rast.html">v.what.rast</a>,
+<a href="v.db.addtable.html">v.db.addtable</a>,
+<a href="v.db.select.html">v.db.select</a>,
+<a href="v.univar.html">v.univar</a>
+</em>
+
+<h2>AUTHOR</h2>
+Soeren Gebbert, heavily based on v.what.rast by Radim Blazed
+
+<p>
+<i>Last changed: $Date: 2012-02-08 22:40:06 +0100 (Mi, 08 Feb 2012) $</i>
More information about the grass-commit
mailing list