GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 
7 /*------------------------------------------------------------------------------------------------*/
8 void
9 P_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 
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, &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 
183  db_commit_transaction(driver);
184 
185  return;
186 }
187 
188 
189 /*------------------------------------------------------------------------------------------------*/
190 int 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 }
db_begin_transaction
int db_begin_transaction(dbDriver *)
Begin transaction.
Definition: c_execute.c:56
bound_box
Bounding box.
Definition: dig_structs.h:65
G_fatal_error
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
if
if(!(yy_init))
Definition: sqlp.yy.c:775
bound_box::W
double W
West.
Definition: dig_structs.h:82
Cell_head
2D/3D raster map header (used also for region)
Definition: gis.h:412
db_init_string
void db_init_string(dbString *)
Initialize dbString.
Definition: string.c:25
dataInterpolateBilin
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:644
Cell_head::west
double west
Extent coordinates (west)
Definition: gis.h:464
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
Segment_get
int Segment_get(SEGMENT *, void *, off_t, off_t)
Get value from segment file.
Definition: segment/get.c:39
bound_box::E
double E
East.
Definition: dig_structs.h:78
GV_POINT
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
line_cats
Feature category info.
Definition: dig_structs.h:1702
bound_box::N
double N
North.
Definition: dig_structs.h:70
line_pnts::x
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
line_pnts
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
Cell_head::north
double north
Extent coordinates (north)
Definition: gis.h:458
db_append_string
int db_append_string(dbString *, const char *)
Append string to dbString.
Definition: string.c:205
Vect_new_line_struct
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
Cell_head::south
double south
Extent coordinates (south)
Definition: gis.h:460
db_execute_immediate
int db_execute_immediate(dbDriver *, dbString *)
Execute SQL statements.
Definition: c_execute.c:27
Segment_put
int Segment_put(SEGMENT *, const void *, off_t, off_t)
Definition: segment/put.c:45
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:178
Map_info
Vector map info.
Definition: dig_structs.h:1259
Y
#define Y
Definition: ogsf.h:138
line_pnts::y
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
_db_driver
Definition: dbmi.h:168
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
db_commit_transaction
int db_commit_transaction(dbDriver *)
Commit transaction.
Definition: c_execute.c:82
driver
Definition: driver.h:22
Rast_col_to_easting
double Rast_col_to_easting(double, const struct Cell_head *)
Column to easting.
Definition: raster/window.c:264
Cell_head::ns_res
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:452
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
bound_box::S
double S
South.
Definition: dig_structs.h:74
dataInterpolateBicubic
double dataInterpolateBicubic(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:533
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:190
line_pnts::z
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
_
#define _(str)
Definition: glocale.h:10
X
#define X
Definition: ogsf.h:137
mean
float mean(IClass_statistics *statistics, int band)
Helper function for computing mean.
Definition: iclass_statistics.c:372
_db_string
Definition: dbmi.h:147
Cell_head::ew_res
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:448
DB_OK
#define DB_OK
Definition: dbmi.h:71
Rast_row_to_northing
double Rast_row_to_northing(double, const struct Cell_head *)
Row to northing.
Definition: raster/window.c:244