GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-498676c526
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;
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 
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 
306 int 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
#define DB_C_TYPE_INT
Definition: dbmi.h:108
#define DB_C_TYPE_DOUBLE
Definition: dbmi.h:109
#define DB_OK
Definition: dbmi.h:71
void db_CatValArray_free(dbCatValArray *)
Free allocated dbCatValArray.
Definition: value.c:373
int db_CatValArray_get_value_int(dbCatValArray *, int, int *)
Find value (integer) by key.
void db_CatValArray_init(dbCatValArray *)
Initialize dbCatValArray.
Definition: value.c:361
int db_column_Ctype(dbDriver *, const char *, const char *)
Get column ctype.
dbDriver * db_start_driver(const char *)
Initialize a new dbDriver for db transaction.
Definition: start.c:51
int db_open_database(dbDriver *, dbHandle *)
Open database connection.
Definition: c_opendb.c:27
int db_close_database_shutdown_driver(dbDriver *)
Close driver/database connection.
Definition: db.c:61
int db_select_CatValArray(dbDriver *, const char *, const char *, const char *, const char *, dbCatValArray *)
Select pairs key/value to array, values are sorted by key (must be integer)
int db_set_handle(dbHandle *, const char *, const char *)
Set handle (database and schema name)
Definition: handle.c:39
void db_init_handle(dbHandle *)
Initialize handle (i.e database/schema)
Definition: handle.c:23
void db_init_string(dbString *)
Initialize dbString.
Definition: string.c:25
int db_CatValArray_get_value_double(dbCatValArray *, int, double *)
Find value (double) by key.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
void void G_verbose_message(const char *,...) __attribute__((format(printf
void void void G_important_message(const char *,...) __attribute__((format(printf
void G_message(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
const char * Vect_get_full_name(struct Map_info *)
Get fully qualified name of vector map.
int Vect_cat_get(const struct line_cats *, int, int *)
Get first found category of given field.
struct field_info * Vect_get_field(struct Map_info *, int)
Get information about link to database (by layer number)
Definition: field.c:515
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_read_next_line(struct Map_info *, struct line_pnts *, struct line_cats *)
Read next vector feature.
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
int Vect_is_3d(struct Map_info *)
Check if vector map is 3D.
#define GV_LINE
Definition: dig_defines.h:184
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:183
#define GV_BOUNDARY
Definition: dig_defines.h:185
#define _(str)
Definition: glocale.h:10
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
void free(void *)
Vector map info.
Definition: dig_structs.h:1243
Definition: driver.h:21
Layer (old: field) information.
Definition: dig_structs.h:131
double zmult
Definition: interpf.h:67
const char * wheresql
Definition: interpf.h:137
Feature category info.
Definition: dig_structs.h:1677
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
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
int n_points
Definition: dataquad.h:52
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