grass  7.8.6
About: GRASS (Geographic Resources Analysis Support System) is a raster- and vector-based GIS, image processing system, graphics production system and spatial modeling system.
  Fossies Dox: grass-7.8.6.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

raster.c
Go to the documentation of this file.
1#include <grass/config.h>
2#include <stdio.h>
3#include <stdlib.h>
4#include <string.h>
5#include <grass/lidar.h>
6
7/*------------------------------------------------------------------------------------------------*/
8void
9P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
10 struct bound_box General, struct bound_box Overlap, double **obs,
11 double *param, int *line_num, double pe, double pn,
12 double overlap, int nsplx, int nsply, int num_points,
13 int bilin, struct line_cats *categories, dbDriver * driver,
14 double mean, char *tab_name)
15{
16 int i;
17 char buf[1024];
18 dbString sql;
19
20 double interpolation, csi, eta, weight;
21 struct line_pnts *point;
22
24
26
27 for (i = 0; i < num_points; i++) {
28
29 if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) { /*Here mean is just for asking if obs point is in box */
30
31 if (bilin)
32 interpolation =
33 dataInterpolateBilin(obs[i][0], obs[i][1], pe, pn, nsplx,
34 nsply, Elaboration->west,
35 Elaboration->south, param);
36 else
37 interpolation =
38 dataInterpolateBicubic(obs[i][0], obs[i][1], pe, pn,
39 nsplx, nsply, Elaboration->west,
40 Elaboration->south, param);
41
42 interpolation += mean;
43 Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1],
44 &interpolation, 1);
45
46 if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation, &Overlap)) { /*(5) */
47 Vect_write_line(Out, GV_POINT, point, categories);
48 }
49 else {
50 db_init_string(&sql);
51
52 sprintf(buf, "INSERT INTO %s (ID, X, Y, Interp)", tab_name);
53 db_append_string(&sql, buf);
54
55 sprintf(buf, " VALUES (");
56 db_append_string(&sql, buf);
57 sprintf(buf, "%d, %f, %f, ", line_num[i], obs[i][0],
58 obs[i][1]);
59 db_append_string(&sql, buf);
60
61 if ((*point->x > Overlap.E) && (*point->x < General.E)) {
62 if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
63 csi = (General.E - *point->x) / overlap;
64 eta = (General.N - *point->y) / overlap;
65 weight = csi * eta;
66 *point->z = weight * interpolation;
67
68 sprintf(buf, "%lf", *point->z);
69 db_append_string(&sql, buf);
70 sprintf(buf, ")");
71 db_append_string(&sql, buf);
72
73 if (db_execute_immediate(driver, &sql) != DB_OK)
74 G_fatal_error(_("Unable to access table <%s>"),
75 tab_name);
76 }
77 else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
78 csi = (General.E - *point->x) / overlap;
79 eta = (*point->y - General.S) / overlap;
80 weight = csi * eta;
81 *point->z = weight * interpolation;
82
83 sprintf(buf, "%lf", *point->z);
84 db_append_string(&sql, buf);
85 sprintf(buf, ")");
86 db_append_string(&sql, buf);
87
88 if (db_execute_immediate(driver, &sql) != DB_OK)
89 G_fatal_error(_("Unable to access table <%s>"),
90 tab_name);
91 }
92 else if ((*point->y <= Overlap.N) && (*point->y >= Overlap.S)) { /*(1) */
93 weight = (General.E - *point->x) / overlap;
94 *point->z = weight * interpolation;
95
96 sprintf(buf, "%lf", *point->z);
97 db_append_string(&sql, buf);
98 sprintf(buf, ")");
99 db_append_string(&sql, buf);
100
101 if (db_execute_immediate(driver, &sql) != DB_OK)
102 G_fatal_error(_("Unable to access table <%s>"),
103 tab_name);
104 }
105 }
106 else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
107 if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(4) */
108 csi = (*point->x - General.W) / overlap;
109 eta = (General.N - *point->y) / overlap;
110 weight = eta * csi;
111 *point->z = weight * interpolation;
112
113 sprintf(buf, "%lf", *point->z);
114 db_append_string(&sql, buf);
115 sprintf(buf, ")");
116 db_append_string(&sql, buf);
117
118 if (db_execute_immediate(driver, &sql) != DB_OK)
119 G_fatal_error(_("Unable to access table <%s>"),
120 tab_name);
121 }
122 else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(2) */
123 csi = (*point->x - General.W) / overlap;
124 eta = (*point->y - General.S) / overlap;
125 weight = csi * eta;
126 *point->z = weight * interpolation;
127
128 sprintf(buf, "%lf", *point->z);
129 db_append_string(&sql, buf);
130 sprintf(buf, ")");
131 db_append_string(&sql, buf);
132
133 if (db_execute_immediate(driver, &sql) != DB_OK)
134 G_fatal_error(_("Unable to access table <%s>"),
135 tab_name);
136 }
137 else if ((*point->y >= Overlap.S) && (*point->y <= Overlap.N)) { /*(2) */
138 weight = (*point->x - General.W) / overlap;
139 *point->z = weight * interpolation;
140
141 sprintf(buf, "%lf", *point->z);
142 db_append_string(&sql, buf);
143 sprintf(buf, ")");
144 db_append_string(&sql, buf);
145
146 if (db_execute_immediate(driver, &sql) != DB_OK)
147 G_fatal_error(_("Unable to access table <%s>"),
148 tab_name);
149 }
150 }
151 else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)){
152 if ((*point->y > Overlap.N) && (*point->y < General.N)) { /*(3) */
153 weight = (General.N - *point->y) / overlap;
154 *point->z = weight * interpolation;
155
156 sprintf(buf, "%lf", *point->z);
157 db_append_string(&sql, buf);
158 sprintf(buf, ")");
159 db_append_string(&sql, buf);
160
161 if (db_execute_immediate(driver, &sql) != DB_OK)
162 G_fatal_error(_("Unable to access table <%s>"),
163 tab_name);
164 }
165 else if ((*point->y < Overlap.S) && (*point->y > General.S)) { /*(1) */
166 weight = (*point->y - General.S) / overlap;
167 *point->z = (1 - weight) * interpolation;
168
169 sprintf(buf, "%lf", *point->z);
170 db_append_string(&sql, buf);
171 sprintf(buf, ")");
172 db_append_string(&sql, buf);
173
174 if (db_execute_immediate(driver, &sql) != DB_OK)
175 G_fatal_error(_("Unable to access table <%s>"),
176 tab_name);
177 }
178 }
179 }
180 } /*IF*/
181 } /*FOR*/
182
184
185 return;
186}
187
188
189/*------------------------------------------------------------------------------------------------*/
190int P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original,
191 struct bound_box General, struct bound_box Overlap,
192 SEGMENT *out_seg, double *param,
193 double passoN, double passoE, double overlap,
194 double mean, int nsplx, int nsply,
195 int nrows, int ncols, int bilin)
196{
197
198 int col, row, startcol, endcol, startrow, endrow;
199 double X, Y, interpolation, weight, csi, eta, dval;
200
201 /* G_get_window(&Original); */
202 if (Original->north > General.N)
203 startrow = (Original->north - General.N) / Original->ns_res - 1;
204 else
205 startrow = 0;
206 if (Original->north > General.S) {
207 endrow = (Original->north - General.S) / Original->ns_res + 1;
208 if (endrow > nrows)
209 endrow = nrows;
210 }
211 else
212 endrow = nrows;
213 if (General.W > Original->west)
214 startcol = (General.W - Original->west) / Original->ew_res - 1;
215 else
216 startcol = 0;
217 if (General.E > Original->west) {
218 endcol = (General.E - Original->west) / Original->ew_res + 1;
219 if (endcol > ncols)
220 endcol = ncols;
221 }
222 else
223 endcol = ncols;
224
225 for (row = startrow; row < endrow; row++) {
226 for (col = startcol; col < endcol; col++) {
227
228 X = Rast_col_to_easting((double)(col) + 0.5, Original);
229 Y = Rast_row_to_northing((double)(row) + 0.5, Original);
230
231 if (Vect_point_in_box(X, Y, mean, &General)) { /* Here, mean is just for asking if obs point is in box */
232
233 if (bilin)
234 interpolation =
235 dataInterpolateBilin(X, Y, passoE, passoN, nsplx,
236 nsply, Elaboration->west,
237 Elaboration->south, param);
238 else
239 interpolation =
240 dataInterpolateBicubic(X, Y, passoE, passoN, nsplx,
241 nsply, Elaboration->west,
242 Elaboration->south, param);
243
244 interpolation += mean;
245
246 if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
247 dval = interpolation;
248 }
249 else {
250 Segment_get(out_seg, &dval, row, col);
251 if ((X > Overlap.E) && (X < General.E)) {
252 if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
253 csi = (General.E - X) / overlap;
254 eta = (General.N - Y) / overlap;
255 weight = csi * eta;
256 interpolation *= weight;
257 dval += interpolation;
258 }
259 else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
260 csi = (General.E - X) / overlap;
261 eta = (Y - General.S) / overlap;
262 weight = csi * eta;
263 interpolation *= weight;
264 dval = interpolation;
265 }
266 else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (1) */
267 weight = (General.E - X ) / overlap;
268 interpolation *= weight;
269 dval = interpolation;
270 }
271 }
272 else if ((X < Overlap.W) && (X > General.W)) {
273 if ((Y > Overlap.N) && (Y < General.N)) { /* (4) */
274 csi = (X - General.W) / overlap;
275 eta = (General.N - Y) / overlap;
276 weight = eta * csi;
277 interpolation *= weight;
278 dval += interpolation;
279 }
280 else if ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
281 csi = (X - General.W) / overlap;
282 eta = (Y - General.S) / overlap;
283 weight = csi * eta;
284 interpolation *= weight;
285 dval += interpolation;
286 }
287 else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (2) */
288 weight = (X - General.W) / overlap;
289 interpolation *= weight;
290 dval += interpolation;
291 }
292 }
293 else if ((X >= Overlap.W) && (X <= Overlap.E)) {
294 if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
295 weight = (General.N - Y) / overlap;
296 interpolation *= weight;
297 dval += interpolation;
298 }
299 else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
300 weight = (Y - General.S) / overlap;
301 interpolation *= weight;
302 dval = interpolation;
303 }
304 }
305 }
306 Segment_put(out_seg, &dval, row, col);
307 }
308 } /* END COL */
309 } /* END ROW */
310 return 1;
311}
double dataInterpolateBicubic(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:533
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:644
int db_begin_transaction(dbDriver *driver)
Begin transaction.
Definition: c_execute.c:56
int db_execute_immediate(dbDriver *driver, dbString *SQLstatement)
Execute SQL statements.
Definition: c_execute.c:27
int db_commit_transaction(dbDriver *driver)
Commit transaction.
Definition: c_execute.c:82
if(!DBFLoadRecord(psDBF, hEntity)) return NULL
#define DB_OK
Definition: dbmi.h:71
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: error.c:160
#define _(str)
Definition: glocale.h:13
float mean(IClass_statistics *statistics, int band)
Helper function for computing mean.
void P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration, struct bound_box General, struct bound_box Overlap, double **obs, double *param, int *line_num, double pe, double pn, double overlap, int nsplx, int nsply, int num_points, int bilin, struct line_cats *categories, dbDriver *driver, double mean, char *tab_name)
Definition: raster.c:9
int P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original, struct bound_box General, struct bound_box Overlap, SEGMENT *out_seg, double *param, double passoN, double passoE, double overlap, double mean, int nsplx, int nsply, int nrows, int ncols, int bilin)
Definition: raster.c:190
int Vect_copy_xyz_to_pnts(struct line_pnts *Points, const double *x, const double *y, const double *z, int n)
Copy points from array to line_pnts structure.
Definition: line.c:99
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
#define X
Definition: ogsf.h:137
#define Y
Definition: ogsf.h:138
static Rast3d_paramType * param
Definition: param.c:21
static int ncols
Definition: raster.c:29
double Rast_row_to_northing(double row, const struct Cell_head *window)
Row to northing.
Definition: window.c:244
double Rast_col_to_easting(double col, const struct Cell_head *window)
Column to easting.
Definition: window.c:264
int Segment_get(SEGMENT *SEG, void *buf, off_t row, off_t col)
Get value from segment file.
Definition: get.c:39
int Segment_put(SEGMENT *SEG, const void *buf, off_t row, off_t col)
Definition: put.c:45
void db_init_string(dbString *x)
Initialize dbString.
Definition: string.c:25
int db_append_string(dbString *x, const char *s)
Append string to dbString.
Definition: string.c:205
2D/3D raster map header (used also for region)
Definition: gis.h:414
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:449
double north
Extent coordinates (north)
Definition: gis.h:459
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:453
double south
Extent coordinates (south)
Definition: gis.h:461
double west
Extent coordinates (west)
Definition: gis.h:465
Vector map info.
Definition: dig_structs.h:1260
Bounding box.
Definition: dig_structs.h:66
double W
West.
Definition: dig_structs.h:82
double S
South.
Definition: dig_structs.h:74
double N
North.
Definition: dig_structs.h:70
double E
East.
Definition: dig_structs.h:78
Definition: driver.h:23
Feature category info.
Definition: dig_structs.h:1703
Feature geometry info - coordinates.
Definition: dig_structs.h:1676
Definition: plot.c:44
int y
Definition: plot.c:46
double x
Definition: plot.c:45
int Vect_point_in_box(double x, double y, double z, const struct bound_box *Box)
Tests if point is in 3D box.
Definition: box.c:48
off_t Vect_write_line(struct Map_info *Map, int type, const struct line_pnts *points, const struct line_cats *cats)
Writes a new feature.
Definition: write.c:178