GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71785
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  */
51 int 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;
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 
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 =
183  &intval);
184  sm = intval;
185  }
186  else { /* DB_C_TYPE_DOUBLE */
187  ret =
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 
306 int 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 }
int db_select_CatValArray(dbDriver *driver, const char *tab, const char *key, const char *col, const char *where, dbCatValArray *cvarr)
Select pairs key/value to array, values are sorted by key (must be integer)
CELL cat
Definition: raster3d/cats.c:82
int db_column_Ctype(dbDriver *driver, const char *tab, const char *col)
Get column ctype.
int db_CatValArray_get_value_int(dbCatValArray *arr, int key, int *val)
Find value (integer) by key.
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:130
double y_orig
Definition: dataquad.h:51
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
int db_close_database_shutdown_driver(dbDriver *driver)
Close driver/database connection.
Definition: db.c:62
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
const char * Vect_get_full_name(const struct Map_info *Map)
Get fully qualified name of vector map.
void db_CatValArray_free(dbCatValArray *arr)
Free allocated dbCatValArray.
Definition: value.c:373
double amax1(double, double)
Definition: minmax.c:54
int process_point(double, double, double, double, struct tree_info *, double, double *, double *, double *, double *, double *, double *, int *, int *, int *)
Definition: vinput2d.c:306
void free(void *)
#define NULL
Definition: ccmath.h:32
#define x
double x_orig
Definition: dataquad.h:50
int translate_quad(struct multtree *tree, double numberx, double numbery, double numberz, int n_leafs)
Definition: input2d.c:92
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition: dataquad.c:38
const char * wheresql
Definition: interpf.h:104
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
int db_CatValArray_get_value_double(dbCatValArray *arr, int key, double *val)
Find value (double) by key.
struct quaddata * data
Definition: qtree.h:58
Feature category info.
Definition: dig_structs.h:1702
#define GV_LINE
Definition: dig_defines.h:183
Layer (old: field) information.
Definition: dig_structs.h:134
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
int Vect_is_3d(const struct Map_info *Map)
Check if vector map is 3D.
int Vect_cat_get(const struct line_cats *Cats, int field, int *cat)
Get first found category of given field.
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
#define GV_BOUNDARY
Definition: dig_defines.h:184
#define DB_C_TYPE_INT
Definition: dbmi.h:108
int db_set_handle(dbHandle *handle, const char *dbName, const char *dbSchema)
Set handle (database and schema name)
Definition: handle.c:39
double ymax
Definition: dataquad.h:53
Vector map info.
Definition: dig_structs.h:1259
double amin1(double, double)
Definition: minmax.c:67
void db_CatValArray_init(dbCatValArray *arr)
Initialize dbCatValArray.
Definition: value.c:361
struct line_cats * Vect_new_cats_struct()
Creates and initializes line_cats structure.
Definition: driver.h:22
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:430
#define DB_C_TYPE_DOUBLE
Definition: dbmi.h:109
int Vect_read_next_line(const struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c)
Read next vector feature.
void db_init_handle(dbHandle *handle)
Initialize handle (i.e database/schema)
Definition: handle.c:23
#define _(str)
Definition: glocale.h:13
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:89
double zmult
Definition: interpf.h:70
int MT_insert(struct triple *point, struct tree_info *info, struct multtree *tree, int n_leafs)
Definition: qtree.c:105
double xmax
Definition: dataquad.h:52
int
Reads the categories file for map name in mapset and stores the categories in the pcats structure...
int n_rows
Definition: dataquad.h:54
int db_open_database(dbDriver *driver, dbHandle *handle)
Open database connection.
Definition: c_opendb.c:27
struct multtree * root
Definition: qtree.h:53
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:426
struct field_info * Vect_get_field(const struct Map_info *Map, int field)
Get information about link to database (by layer number)
Definition: field.c:461
int IL_vector_input_data_2d(struct interp_params *, struct Map_info *, int, char *, char *, struct tree_info *, double *, double *, double *, double *, double *, double *, int *, double *)
Definition: vinput2d.c:51
int n_cols
Definition: dataquad.h:55
dbDriver * db_start_driver(const char *name)
Initialize a new dbDriver for db transaction.
Definition: start.c:50
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
void db_init_string(dbString *x)
Initialize dbString.
Definition: string.c:25
#define DB_OK
Definition: dbmi.h:71