GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
vinput2d.c
Go to the documentation of this file.
1 
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 
52  struct Map_info *Map,
53  int field,
54  char *zcol,
55  char *scol,
56  struct tree_info *info,
57  double *xmin, double *xmax,
58  double *ymin, double *ymax,
59  double *zmin, double *zmax,
60  int *n_points,
61  double *dmax
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 
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 }
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
double amax1(double, double)
Definition: minmax.c:54
#define NULL
Definition: ccmath.h:32
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
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
struct quaddata * data
Definition: qtree.h:58
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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
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
double ymax
Definition: dataquad.h:53
double amin1(double, double)
Definition: minmax.c:67
Definition: driver.h:22
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 n_rows
Definition: dataquad.h:54
struct multtree * root
Definition: qtree.h:53
int n_cols
Definition: dataquad.h:55
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203