1 /* ***************************************************************
2 *
3 * MODULE: v.what.rast3
4 *
5 * AUTHOR(S): Soeren Gebbert, this code is based on a slightly modified version of
6 * v.what.rast from Radim Blazek and Michael Shapiro
7 *
8 * PURPOSE: Uploads 3d raster values at positions of vector points to the table
9 *
10 * COPYRIGHT: (C) 2001, 2011 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General
13 * Public License (>=v2). Read the file COPYING that
14 * comes with GRASS for details.
15 *
16 * **************************************************************/
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <string.h>
20
21 #include <grass/raster3d.h>
22 #include <grass/dbmi.h>
23 #include <grass/vector.h>
24 #include <grass/glocale.h>
25
26 #include "local_proto.h"
27
main(int argc,char * argv[])28 int main(int argc, char *argv[])
29 {
30 int i, j, nlines, type, field, cat;
31
32 char buf[2048];
33 struct {
34 struct Option *vect, *rast3d, *field, *col, *where;
35 } opt;
36 int Cache_size;
37 struct order *cache;
38 struct GModule *module;
39
40 RASTER3D_Region region;
41 RASTER3D_Map *map;
42
43 struct Map_info Map;
44 struct line_pnts *Points;
45 struct line_cats *Cats;
46 int point;
47 int point_cnt; /* number of points in cache */
48 int outside_cnt; /* points outside region */
49 int nocat_cnt; /* points inside region but without category */
50 int dupl_cnt; /* duplicate categories */
51 int typeIntern;
52 int is_empty;
53 struct bound_box box;
54
55 int *catexst, *cex;
56 struct field_info *Fi;
57 dbString stmt;
58 dbDriver *driver;
59 int select, norec_cnt, update_cnt, upderr_cnt, col_type;
60
61 G_gisinit(argv[0]);
62
63 module = G_define_module();
64 G_add_keyword(_("vector"));
65 G_add_keyword(_("sampling"));
66 G_add_keyword(_("raster"));
67 G_add_keyword(_("position"));
68 G_add_keyword(_("querying"));
69 G_add_keyword(_("attribute table"));
70 G_add_keyword(_("surface information"));
71 module->description =
72 _("Uploads 3D raster values at positions of vector points to the table.");
73
74 opt.vect = G_define_standard_option(G_OPT_V_MAP);
75 opt.vect->label =
76 _("Name of vector points map for which to edit attributes");
77
78 opt.field = G_define_standard_option(G_OPT_V_FIELD);
79
80 opt.rast3d = G_define_standard_option(G_OPT_R3_MAP);
81 opt.rast3d->key = "raster_3d";
82 opt.rast3d->description = _("Name of existing 3D raster map to be queried");
83
84 opt.col = G_define_standard_option(G_OPT_DB_COLUMN);
85 opt.col->required = YES;
86 opt.col->description =
87 _("Name of attribute column to be updated with the query result");
88
89 opt.where = G_define_standard_option(G_OPT_DB_WHERE);
90
91 if (G_parser(argc, argv))
92 exit(EXIT_FAILURE);
93
94 db_init_string(&stmt);
95 Points = Vect_new_line_struct();
96 Cats = Vect_new_cats_struct();
97
98 /* Initiate the default settings */
99 Rast3d_init_defaults();
100
101 /* Figure out the current region settings */
102 Rast3d_get_window(®ion);
103
104 box.N = region.north;
105 box.S = region.south;
106 box.E = region.east;
107 box.W = region.west;
108 box.T = region.top;
109 box.B = region.bottom;
110
111 /* Open vector */
112 Vect_set_open_level(2);
113 if (Vect_open_old2(&Map, opt.vect->answer, "", opt.field->answer) < 0)
114 G_fatal_error(_("Unable to open vector map <%s>"), opt.vect->answer);
115
116 field = Vect_get_field_number(&Map, opt.field->answer);
117
118 Fi = Vect_get_field(&Map, field);
119 if (Fi == NULL)
120 G_fatal_error(_("Database connection not defined for layer %d"),
121 field);
122
123 /* Open driver */
124 driver = db_start_driver_open_database(Fi->driver, Fi->database);
125 if (driver == NULL) {
126 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
127 Fi->database, Fi->driver);
128 }
129 db_set_error_handler_driver(driver);
130
131 map = Rast3d_open_cell_old(opt.rast3d->answer, G_find_raster3d(opt.rast3d->answer, ""), ®ion,
132 RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
133
134 if (map == NULL)
135 G_fatal_error(_("Error opening 3D raster map <%s>"), opt.rast3d->answer);
136
137 /* Check column type */
138 col_type = db_column_Ctype(driver, Fi->table, opt.col->answer);
139
140 if (col_type == -1)
141 G_fatal_error(_("Column <%s> not found"), opt.col->answer);
142
143 if (col_type != DB_C_TYPE_DOUBLE)
144 G_fatal_error(_("Column type not supported, please use a column with double type"));
145
146 /* Read vector points to cache */
147 Cache_size = Vect_get_num_primitives(&Map, GV_POINT);
148 /* Note: Some space may be wasted (outside region or no category) */
149
150 cache = (struct order *)G_calloc(Cache_size, sizeof(struct order));
151
152 point_cnt = outside_cnt = nocat_cnt = 0;
153
154 nlines = Vect_get_num_lines(&Map);
155
156 G_debug(1, "Reading %d vector features from map", nlines);
157
158 G_important_message(_("Reading features from vector map..."));
159 for (i = 1; i <= nlines; i++) {
160 type = Vect_read_line(&Map, Points, Cats, i);
161 G_debug(4, "line = %d type = %d", i, type);
162
163 G_percent(i, nlines, 2);
164
165 /* check type */
166 if (!(type & GV_POINT))
167 continue; /* Points only */
168
169 /* check region */
170 if (!Vect_point_in_box(Points->x[0], Points->y[0], Points->z[0], &box)) {
171 outside_cnt++;
172 continue;
173 }
174
175 Vect_cat_get(Cats, field, &cat);
176 if (cat < 0) { /* no category of given field */
177 nocat_cnt++;
178 continue;
179 }
180
181 G_debug(4, " cat = %d", cat);
182
183 cache[point_cnt].x = Points->x[0];
184 cache[point_cnt].y = Points->y[0];
185 cache[point_cnt].z = Points->z[0];
186 cache[point_cnt].cat = cat;
187 cache[point_cnt].count = 1;
188 point_cnt++;
189 }
190
191 Vect_set_db_updated(&Map);
192 Vect_hist_command(&Map);
193 Vect_set_db_updated(&Map);
194 Vect_close(&Map);
195
196 G_debug(1, "Read %d vector points", point_cnt);
197 /* Cache may contain duplicate categories, sort by cat, find and remove duplicates
198 * and recalc count and decrease point_cnt */
199 qsort(cache, point_cnt, sizeof(struct order), by_cat);
200
201 G_debug(1, "Points are sorted, starting duplicate removal loop");
202
203 for (i = 0, j = 1; j < point_cnt; j++)
204 if (cache[i].cat != cache[j].cat)
205 cache[++i] = cache[j];
206 else
207 cache[i].count++;
208 point_cnt = i + 1;
209
210 G_debug(1, "%d vector points left after removal of duplicates",
211 point_cnt);
212
213 /* Report number of points not used */
214 if (outside_cnt)
215 G_warning(n_("%d point outside current region was skipped",
216 "%d points outside current region were skipped",
217 outside_cnt), outside_cnt);
218
219 if (nocat_cnt)
220 G_warning(n_("%d point without category was skipped",
221 "%d points without category were skipped",
222 nocat_cnt), nocat_cnt);
223
224 /* Extract raster values from file and store in cache */
225 G_debug(1, "Extracting raster values");
226
227 typeIntern = Rast3d_tile_type_map(map);
228
229 /* Sample the 3d raster map */
230 for (point = 0; point < point_cnt; point++) {
231
232 if (cache[point].count > 1)
233 continue; /* duplicate cats */
234
235 if(typeIntern == FCELL_TYPE) {
236 Rast3d_get_window_value(map, cache[point].y, cache[point].x, cache[point].z,
237 &cache[point].fvalue, FCELL_TYPE);
238 } else {
239 Rast3d_get_window_value(map, cache[point].y, cache[point].x, cache[point].z,
240 &cache[point].dvalue, DCELL_TYPE);
241 }
242 }
243
244 Rast3d_close(map);
245
246 /* Update table from cache */
247 G_debug(1, "Updating db table");
248
249 /* select existing categories to array (array is sorted) */
250 select = db_select_int(driver, Fi->table, Fi->key, NULL, &catexst);
251
252 db_begin_transaction(driver);
253
254 norec_cnt = update_cnt = upderr_cnt = dupl_cnt = 0;
255
256 G_message("Update vector attributes...");
257 for (point = 0; point < point_cnt; point++) {
258 if (cache[point].count > 1) {
259 G_warning(_("More points (%d) of category %d, value set to 'NULL'"),
260 cache[point].count, cache[point].cat);
261 dupl_cnt++;
262 }
263
264 G_percent(point, point_cnt, 2);
265
266 /* category exist in DB ? */
267 cex =
268 (int *)bsearch((void *)&(cache[point].cat), catexst, select,
269 sizeof(int), srch_cat);
270 if (cex == NULL) { /* cat does not exist in DB */
271 norec_cnt++;
272 G_warning(_("No record for category %d in table <%s>"),
273 cache[point].cat, Fi->table);
274 continue;
275 }
276
277 G_snprintf(buf, 2048, "update %s set %s = ", Fi->table, opt.col->answer);
278
279 db_set_string(&stmt, buf);
280
281 is_empty = 0;
282
283 if (cache[point].count > 1)
284 is_empty = 1;
285 if(typeIntern == FCELL_TYPE)
286 is_empty = Rast3d_is_null_value_num(&cache[point].fvalue, FCELL_TYPE);
287 if(typeIntern == DCELL_TYPE)
288 is_empty = Rast3d_is_null_value_num(&cache[point].dvalue, DCELL_TYPE);
289
290 if (is_empty) {
291 G_snprintf(buf, 2048, "NULL");
292 }else {
293 if(typeIntern == FCELL_TYPE)
294 G_snprintf(buf, 2048, "%.10f", cache[point].fvalue);
295 if(typeIntern == DCELL_TYPE)
296 G_snprintf(buf, 2048, "%.15f", cache[point].dvalue);
297 }
298
299 db_append_string(&stmt, buf);
300
301 G_snprintf(buf, 2048, " where %s = %d", Fi->key, cache[point].cat);
302
303 db_append_string(&stmt, buf);
304 /* user provides where condition: */
305 if (opt.where->answer) {
306 G_snprintf(buf, 2048, " AND %s", opt.where->answer);
307 db_append_string(&stmt, buf);
308 }
309 G_debug(3, "%s", db_get_string(&stmt));
310
311 /* Update table */
312 if (db_execute_immediate(driver, &stmt) == DB_OK) {
313 update_cnt++;
314 }
315 else {
316 upderr_cnt++;
317 }
318 }
319 G_percent(1, 1, 1);
320
321 G_debug(1, "Committing DB transaction");
322 db_commit_transaction(driver);
323
324 G_free(catexst);
325 db_close_database_shutdown_driver(driver);
326 db_free_string(&stmt);
327
328 /* Report */
329 G_verbose_message(n_("%d category loaded from table",
330 "%d categories loaded from table",
331 select), select);
332 G_verbose_message(n_("%d category loaded from vector",
333 "%d categories loaded from vector",
334 point_cnt), point_cnt);
335 G_verbose_message(n_("%d category from vector missing in table",
336 "%d categories from vector missing in table",
337 norec_cnt), norec_cnt);
338 if (dupl_cnt > 0)
339 G_message(n_("%d duplicate category in vector",
340 "%d duplicate categories in vector",
341 dupl_cnt), dupl_cnt);
342 if (upderr_cnt > 0)
343 G_warning(n_("%d update error",
344 "%d update errors",
345 upderr_cnt), upderr_cnt);
346
347 G_done_msg(n_("%d record updated.",
348 "%d records updated.",
349 update_cnt), update_cnt);
350
351 exit(EXIT_SUCCESS);
352 }
353