GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-3f124c5c78
point2d.c
Go to the documentation of this file.
1 /*!
2  * \file point2d.c
3  *
4  * \author
5  * Lubos Mitas (original program and various modifications)
6  *
7  * \author
8  * H. Mitasova,
9  * I. Kosinovsky,
10  * D. Gerdes,
11  * D. McCauley
12  * (GRASS4.1 version of the program and GRASS4.2 modifications)
13  *
14  * \author modified by McCauley in August 1995
15  * \author modified by Mitasova in August 1995, Nov. 1996
16  *
17  * \copyright
18  * (C) 1993-2006 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 for details.
24  */
25 
26 #include <stdio.h>
27 #include <math.h>
28 #include <unistd.h>
29 #include <grass/gis.h>
30 #include <grass/vector.h>
31 #include <grass/dbmi.h>
32 
33 #define POINT2D_C
34 #include <grass/interpf.h>
35 
36 /* needed for AIX */
37 #ifdef hz
38 #undef hz
39 #endif
40 
41 /*!
42  * Checks if interpolating function interp() evaluates correct z-values at
43  * given points. If smoothing is used calculate the maximum error caused
44  * by smoothing.
45  *
46  * *ertot* is a RMS deviation of the interpolated surface.
47  *
48  * \todo
49  * Alternative description:
50  * ...calculate the maximum and RMS deviation caused by smoothing.
51  */
53  struct quaddata *data, /*!< current region */
54  double *b, /*!< solution of linear equations */
55  double *ertot, /*!< total error */
56  double zmin, /*!< min z-value */
57  double dnorm, struct triple skip_point)
58 {
59  int n_points = data->n_points; /* number of points */
60  struct triple *points = data->points; /* points for interpolation */
61  double east = data->xmax;
62  double west = data->x_orig;
63  double north = data->ymax;
64  double south = data->y_orig;
65  double /* rfsta2, errmax, */ h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
66  double skip_err;
67  int /* n1, */ mm, m, cat;
68 
69  /* double fstar2; */
70  int inside;
71 
72  /* Site *site; */
73  char buf[1024];
74 
75  /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
76  G_fatal_error ("Memory error for site struct"); */
77 
78  /* fstar2 = params->fi * params->fi / 4.; */
79  /* errmax = .0; */
80  /* n1 = n_points + 1; */
81  for (mm = 1; mm <= n_points; mm++) {
82  h = b[0];
83  for (m = 1; m <= n_points; m++) {
84  xx = points[mm - 1].x - points[m - 1].x;
85  yy = points[mm - 1].y - points[m - 1].y;
86  r2 = yy * yy + xx * xx;
87  if (r2 != 0.) {
88  /* rfsta2 = fstar2 * r2; */
89  r = r2;
90  h = h + b[m] * params->interp(r, params->fi);
91  }
92  }
93  /* modified by helena january 1997 - normalization of z was
94  removed from segm2d.c and interp2d.c
95  hz = (h * dnorm) + zmin;
96  zz = (points[mm - 1].z * dnorm) + zmin;
97  */
98  hz = h + zmin;
99  zz = points[mm - 1].z + zmin;
100  err = hz - zz;
101  xmm = points[mm - 1].x * dnorm + params->x_orig + west;
102  ymm = points[mm - 1].y * dnorm + params->y_orig + south;
103  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
104  (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
105  inside = 1;
106  else
107  inside = 0;
108 
109  if (params->create_devi) {
110 
111  if (inside) { /* if the point is inside the region */
114 
115  Vect_append_point(Pnts, xmm, ymm, zz);
116  cat = count;
117  Vect_cat_set(Cats2, 1, cat);
119 
121  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
122  db_append_string(&sql2, buf);
123 
124  sprintf(buf, ", %f", err);
125  db_append_string(&sql2, buf);
126  db_append_string(&sql2, ")");
127  G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
128 
132  G_fatal_error("Cannot insert new row: %s",
133  db_get_string(&sql2));
134  }
135  count++;
136  }
137  }
138  (*ertot) += err * err;
139  }
140 
141  /* cv stuff */
142  if (params->cv) {
143  h = b[0];
144  for (m = 1; m <= n_points - 1; m++) {
145  xx = points[m - 1].x - skip_point.x;
146  yy = points[m - 1].y - skip_point.y;
147  r2 = yy * yy + xx * xx;
148  if (r2 != 0.) {
149  /* rfsta2 = fstar2 * r2; */
150  r = r2;
151  h = h + b[m] * params->interp(r, params->fi);
152  }
153  }
154  hz = h + zmin;
155  zz = skip_point.z + zmin;
156  skip_err = hz - zz;
157  xmm = skip_point.x * dnorm + params->x_orig + west;
158  ymm = skip_point.y * dnorm + params->y_orig + south;
159 
160  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
161  (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
162  inside = 1;
163  else
164  inside = 0;
165 
166  if (inside) { /* if the point is inside the region */
169 
170  Vect_append_point(Pnts, xmm, ymm, zz);
171  cat = count;
172  Vect_cat_set(Cats2, 1, cat);
174 
176  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
177  db_append_string(&sql2, buf);
178 
179  sprintf(buf, ", %f", skip_err);
180  db_append_string(&sql2, buf);
181  db_append_string(&sql2, ")");
182  G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
183 
187  G_fatal_error("Cannot insert new row: %s",
188  db_get_string(&sql2));
189  }
190  count++;
191  }
192  } /* cv */
193 
194  return 1;
195 }
#define DB_OK
Definition: dbmi.h:71
int db_shutdown_driver(dbDriver *)
Closedown the driver, and free the driver structure.
Definition: shutdown.c:36
void db_zero_string(dbString *)
Zero string.
Definition: string.c:79
int db_execute_immediate(dbDriver *, dbString *)
Execute SQL statements.
Definition: c_execute.c:27
int db_close_database(dbDriver *)
Close database connection.
Definition: c_closedb.c:26
int db_append_string(dbString *, const char *)
Append string to dbString.
Definition: string.c:205
char * db_get_string(const dbString *)
Get string.
Definition: string.c:140
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_reset_cats(struct line_cats *)
Reset category structure to make sure cats structure is clean to be re-used.
int Vect_cat_set(struct line_cats *, int, int)
Add new field/cat to category structure if doesn't exist yet.
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:183
struct Map_info Map2
struct field_info * ff
dbString sql2
dbDriver * driver2
struct line_cats * Cats2
int count
struct line_pnts * Pnts
int IL_check_at_points_2d(struct interp_params *params, struct quaddata *data, double *b, double *ertot, double zmin, double dnorm, struct triple skip_point)
Definition: point2d.c:52
double b
Definition: r_raster.c:39
double r
Definition: r_raster.c:39
char * table
Name of DB table.
Definition: dig_structs.h:151
interp_fn * interp
Definition: interpf.h:124
double fi
Definition: interpf.h:88
double x_orig
Definition: interpf.h:101
double y_orig
Definition: interpf.h:101
bool create_devi
Definition: interpf.h:114
double ymax
Definition: dataquad.h:49
double y_orig
Definition: dataquad.h:47
double x_orig
Definition: dataquad.h:46
struct triple * points
Definition: dataquad.h:53
int n_points
Definition: dataquad.h:52
double xmax
Definition: dataquad.h:48
double z
Definition: dataquad.h:41
double x
Definition: dataquad.h:39
double y
Definition: dataquad.h:40
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:216