GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
vinput2d.c
Go to the documentation of this file.
1
2/*!
3 * \file vinput2d.c
4 *
5 * \author
6 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993
7 * University of Illinois
8 * US Army Construction Engineering Research Lab
9 *
10 * \author
11 * Mitasova (University of Illinois),
12 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
13 *
14 * \author modified by McCauley in August 1995
15 * \author modified by Mitasova in August 1995
16 * \author modofied by Mitasova in Nov 1999 (dmax fix)
17 *
18 * \copyright
19 * (C) 1993-1999 by Helena Mitasova and the GRASS Development Team
20 *
21 * \copyright
22 * This program is free software under the
23 * GNU General Public License (>=v2).
24 * Read the file COPYING that comes with GRASS
25 * for details.
26 */
27
28#include <stdio.h>
29#include <stdlib.h>
30#include <math.h>
31#include <grass/bitmap.h>
32#include <grass/linkm.h>
33#include <grass/gis.h>
34#include <grass/dbmi.h>
35#include <grass/vector.h>
36#include <grass/glocale.h>
37
38#include <grass/interpf.h>
39
40/*!
41 * Insert into a quad tree
42 *
43 * Inserts input data inside the region into a quad tree. Also translates
44 * data. Returns number of segments in the quad tree.
45 *
46 * As z values may be used (in *Map*):
47 * - z coordinates in 3D file -> field = 0
48 * - categories -> field > 0, zcol = NULL
49 * - attributes -> field > 0, zcol != NULL
50 */
51int IL_vector_input_data_2d(struct interp_params *params, /*!< interpolation parameters */
52 struct Map_info *Map, /*!< input vector map */
53 int field, /*!< category field number */
54 char *zcol, /*!< name of the column containing z values */
55 char *scol, /*!< name of the column containing smooth values */
56 struct tree_info *info, /*!< quadtree info */
57 double *xmin, double *xmax,
58 double *ymin, double *ymax,
59 double *zmin, double *zmax,
60 int *n_points, /*!< number of points used for interpolation */
61 double *dmax /*!< max distance between points */
62 )
63{
64 double dmax2; /* max distance between points squared */
65 double c1, c2, c3, c4;
66 int i, line, k = 0;
67 double ns_res, ew_res;
68 int npoint, OUTRANGE;
69 int totsegm;
70 struct quaddata *data = (struct quaddata *)info->root->data;
71 double xprev, yprev, zprev, x1, y1, z1, d1, xt, yt, z, sm;
72 struct line_pnts *Points;
73 struct line_cats *Cats;
74 int times, j1, ltype, cat, zctype = 0, sctype = 0;
75 struct field_info *Fi;
76 dbDriver *driver;
77 dbHandle handle;
78 dbString stmt;
79 dbCatValArray zarray, sarray;
80
81 OUTRANGE = 0;
82 npoint = 0;
83
84 G_debug(2, "IL_vector_input_data_2d(): field = %d, zcol = %s, scol = %s",
85 field, zcol, scol);
86 ns_res = (data->ymax - data->y_orig) / data->n_rows;
87 ew_res = (data->xmax - data->x_orig) / data->n_cols;
88 dmax2 = *dmax * *dmax;
89
90 Points = Vect_new_line_struct(); /* init line_pnts struct */
91 Cats = Vect_new_cats_struct();
92
93 if (field == 0 && !Vect_is_3d(Map))
94 G_fatal_error(_("Vector map <%s> is not 3D"), Vect_get_full_name(Map));
95
96 if (field > 0 && zcol != NULL) { /* open db driver */
97 G_verbose_message(_("Loading data from attribute table ..."));
98 Fi = Vect_get_field(Map, field);
99 if (Fi == NULL)
100 G_fatal_error(_("Database connection not defined for layer %d"),
101 field);
102 G_debug(3, " driver = %s database = %s table = %s", Fi->driver,
103 Fi->database, Fi->table);
104 db_init_handle(&handle);
105 db_init_string(&stmt);
106 driver = db_start_driver(Fi->driver);
107 db_set_handle(&handle, Fi->database, NULL);
108 if (db_open_database(driver, &handle) != DB_OK)
109 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
110 Fi->database, Fi->driver);
111
112 zctype = db_column_Ctype(driver, Fi->table, zcol);
113 G_debug(3, " zcol C type = %d", zctype);
114 if (zctype == -1)
115 G_fatal_error(_("Column <%s> not found"),
116 zcol);
117 if (zctype != DB_C_TYPE_INT && zctype != DB_C_TYPE_DOUBLE)
118 G_fatal_error(_("Data type of column <%s> must be numeric"), zcol);
119
120 db_CatValArray_init(&zarray);
121 G_debug(3, "RST SQL WHERE: %s", params->wheresql);
122 db_select_CatValArray(driver, Fi->table, Fi->key, zcol,
123 params->wheresql, &zarray);
124
125 if (scol != NULL) {
126 sctype = db_column_Ctype(driver, Fi->table, scol);
127 G_debug(3, " scol C type = %d", sctype);
128 if (sctype == -1)
129 G_fatal_error(_("Column <%s> not found"), scol);
130 if (sctype != DB_C_TYPE_INT && sctype != DB_C_TYPE_DOUBLE)
131 G_fatal_error(_("Data type of column <%s> must be numeric"), scol);
132
133 db_CatValArray_init(&sarray);
134 db_select_CatValArray(driver, Fi->table, Fi->key, scol,
135 params->wheresql, &sarray);
136 }
137
138 db_close_database_shutdown_driver(driver);
139 }
140
141 /* Lines without nodes */
142 G_message(_("Reading features from vector map ..."));
143 sm = 0;
144 line = 1;
145 while ((ltype = Vect_read_next_line(Map, Points, Cats)) != -2) {
146
147 if (!(ltype & (GV_POINT | GV_LINE | GV_BOUNDARY)))
148 continue;
149
150 if (field > 0) { /* use cat or attribute */
151 Vect_cat_get(Cats, field, &cat);
152
153 /* line++; */
154 if (zcol == NULL) { /* use categories */
155 z = (double)cat;
156 }
157 else { /* read att from db */
158 int ret, intval;
159
160 if (zctype == DB_C_TYPE_INT) {
161 ret = db_CatValArray_get_value_int(&zarray, cat, &intval);
162 z = intval;
163 }
164 else { /* DB_C_TYPE_DOUBLE */
165 ret = db_CatValArray_get_value_double(&zarray, cat, &z);
166 }
167
168 if (ret != DB_OK) {
169 if (params->wheresql != NULL)
170 /* G_message(_("Database record for cat %d not used due to SQL statement")); */
171 /* do nothing in this case to not confuse user. Or implement second cat list */
172 ;
173 else
174 G_warning(_("Database record for cat %d not found"),
175 cat);
176 continue;
177 }
178
179 if (scol != NULL) {
180 if (sctype == DB_C_TYPE_INT) {
181 ret =
182 db_CatValArray_get_value_int(&sarray, cat,
183 &intval);
184 sm = intval;
185 }
186 else { /* DB_C_TYPE_DOUBLE */
187 ret =
188 db_CatValArray_get_value_double(&sarray, cat,
189 &sm);
190 }
191 if (sm < 0.0)
192 G_fatal_error(_("Negative value of smoothing detected: sm must be >= 0"));
193 }
194 G_debug(5, " z = %f sm = %f", z, sm);
195 }
196 }
197
198 /* Insert all points including nodes (end points) */
199 for (i = 0; i < Points->n_points; i++) {
200 if (field == 0)
201 z = Points->z[i];
202 process_point(Points->x[i], Points->y[i], z, sm, info,
203 params->zmult, xmin, xmax, ymin, ymax, zmin,
204 zmax, &npoint, &OUTRANGE, &k);
205
206 }
207
208 /* Check all segments */
209 xprev = Points->x[0];
210 yprev = Points->y[0];
211 zprev = Points->z[0];
212 for (i = 1; i < Points->n_points; i++) {
213 /* compare the distance between current and previous */
214 x1 = Points->x[i];
215 y1 = Points->y[i];
216 z1 = Points->z[i];
217
218 xt = x1 - xprev;
219 yt = y1 - yprev;
220 d1 = (xt * xt + yt * yt);
221 if ((d1 > dmax2) && (dmax2 != 0.)) {
222 times = (int)(d1 / dmax2 + 0.5);
223 for (j1 = 0; j1 < times; j1++) {
224 xt = x1 - j1 * ((x1 - xprev) / times);
225 yt = y1 - j1 * ((y1 - yprev) / times);
226 if (field == 0)
227 z = z1 - j1 * ((z1 - zprev) / times);
228
229 process_point(xt, yt, z, sm, info, params->zmult,
230 xmin, xmax, ymin, ymax, zmin, zmax, &npoint,
231 &OUTRANGE, &k);
232 }
233 }
234 xprev = x1;
235 yprev = y1;
236 zprev = z1;
237 }
238 }
239
240 if (field > 0 && zcol != NULL)
241 db_CatValArray_free(&zarray);
242 if (scol != NULL) {
243 db_CatValArray_free(&sarray);
244 }
245
246 c1 = *xmin - data->x_orig;
247 c2 = data->xmax - *xmax;
248 c3 = *ymin - data->y_orig;
249 c4 = data->ymax - *ymax;
250 if ((c1 > 5 * ew_res) || (c2 > 5 * ew_res) || (c3 > 5 * ns_res) ||
251 (c4 > 5 * ns_res)) {
252 static int once = 0;
253
254 if (!once) {
255 once = 1;
256 G_warning(_("Strip exists with insufficient data"));
257 }
258 }
259
260 totsegm =
261 translate_quad(info->root, data->x_orig, data->y_orig, *zmin, 4);
262 if (!totsegm)
263 return 0;
264 data->x_orig = 0;
265 data->y_orig = 0;
266
267 /* G_read_vector_timestamp(name,mapset,ts); */
268
269 if (OUTRANGE > 0)
270 G_warning(_("There are points outside specified 2D/3D region - %d points ignored"),
271 OUTRANGE);
272 if (npoint > 0)
273 G_important_message(_("Ignoring %d points (too dense)"), npoint);
274 npoint = k - npoint - OUTRANGE;
275 if (npoint < params->kmin) {
276 if (npoint != 0) {
277 G_warning(_("%d points given for interpolation (after thinning) is less than given NPMIN=%d"),
278 npoint, params->kmin);
279 params->kmin = npoint;
280 }
281 else {
282 G_warning(_("Zero points in the given region"));
283 return -1;
284 }
285 }
286 if (npoint > params->KMAX2 && params->kmin <= params->kmax) {
287 G_warning(_("Segmentation parameters set to invalid values: npmin= %d, segmax= %d "
288 "for smooth connection of segments, npmin > segmax (see manual)"),
289 params->kmin, params->kmax);
290 return -1;
291 }
292 if (npoint < params->KMAX2 && params->kmax != params->KMAX2)
293 G_warning(_("There are less than %d points for interpolation. No "
294 "segmentation is necessary, to run the program faster set "
295 "segmax=%d (see manual)"), params->KMAX2, params->KMAX2);
296
297 G_verbose_message(_("Number of points from vector map %d"), k);
298 G_verbose_message(_("Number of points outside of 2D/3D region %d"),
299 OUTRANGE);
300 G_verbose_message(_("Number of points being used %d"), npoint);
301
302 *n_points = npoint;
303 return (totsegm);
304}
305
306int process_point(double x, double y, double z, double sm, struct tree_info *info, /* quadtree info */
307 double zmult, /* multiplier for z-values */
308 double *xmin,
309 double *xmax,
310 double *ymin,
311 double *ymax,
312 double *zmin,
313 double *zmax, int *npoint, int *OUTRANGE, int *total)
314{
315 struct triple *point;
316 double c1, c2, c3, c4;
317 int a;
318 static int first_time = 1;
319 struct quaddata *data = (struct quaddata *)info->root->data;
320
321
322 (*total)++;
323
324
325 z = z * zmult;
326 c1 = x - data->x_orig;
327 c2 = data->xmax - x;
328 c3 = y - data->y_orig;
329 c4 = data->ymax - y;
330
331 if (!((c1 >= 0) && (c2 >= 0) && (c3 >= 0) && (c4 >= 0))) {
332 if (!(*OUTRANGE)) {
333 G_warning(_("Some points outside of region (ignored)"));
334 }
335 (*OUTRANGE)++;
336 }
337 else {
338 if (!(point = quad_point_new(x, y, z, sm))) {
339 G_warning(_("Unable to allocate memory"));
340 return -1;
341 }
342 a = MT_insert(point, info, info->root, 4);
343 if (a == 0) {
344 (*npoint)++;
345 }
346 if (a < 0) {
347 G_warning(_("Unable to insert %f,%f,%f a = %d"), x, y, z, a);
348 return -1;
349 }
350 free(point);
351 if (first_time) {
352 first_time = 0;
353 *xmin = x;
354 *ymin = y;
355 *zmin = z;
356 *xmax = x;
357 *ymax = y;
358 *zmax = z;
359 }
360 *xmin = amin1(*xmin, x);
361 *ymin = amin1(*ymin, y);
362 *zmin = amin1(*zmin, z);
363 *xmax = amax1(*xmax, x);
364 *ymax = amax1(*ymax, y);
365 *zmax = amax1(*zmax, z);
366 }
367 return 1;
368}
#define NULL
Definition: ccmath.h:32
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition: dataquad.c:38
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition: gis/error.c:109
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:131
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int translate_quad(struct multtree *tree, double numberx, double numbery, double numberz, int n_leafs)
Definition: input2d.c:92
double amin1(double, double)
Definition: minmax.c:67
double amax1(double, double)
Definition: minmax.c:54
int MT_insert(struct triple *point, struct tree_info *info, struct multtree *tree, int n_leafs)
Definition: qtree.c:105
Definition: driver.h:23
double zmult
Definition: interpf.h:70
const char * wheresql
Definition: interpf.h:104
struct quaddata * data
Definition: qtree.h:58
double ymax
Definition: dataquad.h:53
double y_orig
Definition: dataquad.h:51
double x_orig
Definition: dataquad.h:50
int n_points
Definition: dataquad.h:56
double xmax
Definition: dataquad.h:52
int n_cols
Definition: dataquad.h:55
int n_rows
Definition: dataquad.h:54
struct multtree * root
Definition: qtree.h:53
int process_point(double x, double y, double z, double sm, struct tree_info *info, double zmult, double *xmin, double *xmax, double *ymin, double *ymax, double *zmin, double *zmax, int *npoint, int *OUTRANGE, int *total)
Definition: vinput2d.c:306
int IL_vector_input_data_2d(struct interp_params *params, struct Map_info *Map, int field, char *zcol, char *scol, struct tree_info *info, double *xmin, double *xmax, double *ymin, double *ymax, double *zmin, double *zmax, int *n_points, double *dmax)
Definition: vinput2d.c:51
#define x