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