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