GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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
27#include <stdio.h>
28#include <math.h>
29#include <unistd.h>
30#include <grass/gis.h>
31#include <grass/vector.h>
32#include <grass/dbmi.h>
33
34#define POINT2D_C
35#include <grass/interpf.h>
36
37/* needed for AIX */
38#ifdef hz
39#undef hz
40#endif
41
42/*!
43 * Checks if interpolating function interp() evaluates correct z-values at
44 * given points. If smoothing is used calculate the maximum error caused
45 * by smoothing.
46 *
47 * *ertot* is a RMS deviation of the interpolated surface.
48 *
49 * \todo
50 * Alternative description:
51 * ...calculate the maximum and RMS deviation caused by smoothing.
52 */
54 struct quaddata *data, /*!< current region */
55 double *b, /*!< solution of linear equations */
56 double *ertot, /*!< total error */
57 double zmin, /*!< min z-value */
58 double dnorm,
59 struct triple skip_point
60 )
61{
62 int n_points = data->n_points; /* number of points */
63 struct triple *points = data->points; /* points for interpolation */
64 double east = data->xmax;
65 double west = data->x_orig;
66 double north = data->ymax;
67 double south = data->y_orig;
68 double rfsta2, errmax, h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
69 double skip_err;
70 int n1, mm, m, cat;
71 double fstar2;
72 int inside;
73
74 /* Site *site; */
75 char buf[1024];
76
77
78 /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
79 G_fatal_error ("Memory error for site struct"); */
80
81 fstar2 = params->fi * params->fi / 4.;
82 errmax = .0;
83 n1 = n_points + 1;
84 for (mm = 1; mm <= n_points; mm++) {
85 h = b[0];
86 for (m = 1; m <= n_points; m++) {
87 xx = points[mm - 1].x - points[m - 1].x;
88 yy = points[mm - 1].y - points[m - 1].y;
89 r2 = yy * yy + xx * xx;
90 if (r2 != 0.) {
91 rfsta2 = fstar2 * r2;
92 r = r2;
93 h = h + b[m] * params->interp(r, params->fi);
94 }
95 }
96 /* modified by helena january 1997 - normalization of z was
97 removed from segm2d.c and interp2d.c
98 hz = (h * dnorm) + zmin;
99 zz = (points[mm - 1].z * dnorm) + zmin;
100 */
101 hz = h + zmin;
102 zz = points[mm - 1].z + zmin;
103 err = hz - zz;
104 xmm = points[mm - 1].x * dnorm + params->x_orig + west;
105 ymm = points[mm - 1].y * dnorm + params->y_orig + south;
106 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
107 && (ymm >= south + params->y_orig) &&
108 (ymm <= north + params->y_orig))
109 inside = 1;
110 else
111 inside = 0;
112
113 if (params->create_devi) {
114
115 if (inside) { /* if the point is inside the region */
116 Vect_reset_line(Pnts);
117 Vect_reset_cats(Cats2);
118
119 Vect_append_point(Pnts, xmm, ymm, zz);
120 cat = count;
121 Vect_cat_set(Cats2, 1, cat);
122 Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
123
124 db_zero_string(&sql2);
125 sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
126 db_append_string(&sql2, buf);
127
128 sprintf(buf, ", %f", err);
129 db_append_string(&sql2, buf);
130 db_append_string(&sql2, ")");
131 G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
132
133 if (db_execute_immediate(driver2, &sql2) != DB_OK) {
134 db_close_database(driver2);
135 db_shutdown_driver(driver2);
136 G_fatal_error("Cannot insert new row: %s",
137 db_get_string(&sql2));
138 }
139 count++;
140
141 }
142 }
143 (*ertot) += err * err;
144 }
145
146 /* cv stuff */
147 if (params->cv) {
148 h = b[0];
149 for (m = 1; m <= n_points - 1; m++) {
150 xx = points[m - 1].x - skip_point.x;
151 yy = points[m - 1].y - skip_point.y;
152 r2 = yy * yy + xx * xx;
153 if (r2 != 0.) {
154 rfsta2 = fstar2 * r2;
155 r = r2;
156 h = h + b[m] * params->interp(r, params->fi);
157 }
158 }
159 hz = h + zmin;
160 zz = skip_point.z + zmin;
161 skip_err = hz - zz;
162 xmm = skip_point.x * dnorm + params->x_orig + west;
163 ymm = skip_point.y * dnorm + params->y_orig + south;
164
165 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
166 && (ymm >= south + params->y_orig) &&
167 (ymm <= north + params->y_orig))
168 inside = 1;
169 else
170 inside = 0;
171
172 if (inside) { /* if the point is inside the region */
173 Vect_reset_line(Pnts);
174 Vect_reset_cats(Cats2);
175
176 Vect_append_point(Pnts, xmm, ymm, zz);
177 cat = count;
178 Vect_cat_set(Cats2, 1, cat);
179 Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
180
181 db_zero_string(&sql2);
182 sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
183 db_append_string(&sql2, buf);
184
185 sprintf(buf, ", %f", skip_err);
186 db_append_string(&sql2, buf);
187 db_append_string(&sql2, ")");
188 G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
189
190 if (db_execute_immediate(driver2, &sql2) != DB_OK) {
191 db_close_database(driver2);
192 db_shutdown_driver(driver2);
193 G_fatal_error("Cannot insert new row: %s",
194 db_get_string(&sql2));
195 }
196 count++;
197 }
198 } /* cv */
199
200
201 return 1;
202}
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double b
double r
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
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:53
interp_fn * interp
Definition: interpf.h:101
double fi
Definition: interpf.h:80
double x_orig
Definition: interpf.h:87
double y_orig
Definition: interpf.h:87
bool create_devi
Definition: interpf.h:95
double ymax
Definition: dataquad.h:53
double y_orig
Definition: dataquad.h:51
double x_orig
Definition: dataquad.h:50
struct triple * points
Definition: dataquad.h:57
int n_points
Definition: dataquad.h:56
double xmax
Definition: dataquad.h:52
double z
Definition: dataquad.h:44
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220