GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2653d7b383
Go to the documentation of this file.
1 /*****************************************************************************
2  *
3  * MODULE: Grass PDE Numerical Library
4  * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
5  * soerengebbert <at> gmx <dot> de
6  *
7  * PURPOSE: gradient management functions
8  * part of the gpde library
9  *
10  * COPYRIGHT: (C) 2000 by the GRASS Development Team
11  *
12  * This program is free software under the GNU General Public
13  * License (>=v2). Read the file COPYING that comes with GRASS
14  * for details.
15  *
16  *****************************************************************************/
17
18 #include <grass/N_pde.h>
19
20 /*! \brief Calculate basic statistics of a gradient field
21  *
22  * The statistic is stored in the gradient field struct
23  *
24  * \param field N_gradient_2d_field *
25  * \return void
26  *
27  * */
29 {
30  double minx, miny;
31  double maxx, maxy;
32  double sumx, sumy;
33  int nonullx, nonully;
34
36
37  N_calc_array_2d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
38  N_calc_array_2d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
39
40  if (minx < miny)
41  field->min = minx;
42  else
43  field->min = miny;
44
45  if (maxx > maxy)
46  field->max = maxx;
47  else
48  field->max = maxy;
49
50  field->sum = sumx + sumy;
51  field->nonull = nonullx + nonully;
52  field->mean = field->sum / (double)field->nonull;
53
54  return;
55 }
56
57 /*!
58  * \brief This function computes the gradient based on the input N_array_2d pot
59  * (potential), a weighting factor N_array_2d named weight and the distance
60  between two cells
61  * saved in the N_geom_data struct.
62  *
63  * The gradient is calculated between cells for each cell and direction.
64  * An existing gradient field can be filled with new data or, if a NULL pointer
65  is
66  * given, a new gradient field will be allocated with the appropriate size.
67  *
68  *
69  \verbatim
70  ______________
71  | | | |
72  | | | |
73  |----|-NC-|----|
74  | | | |
75  | WC EC |
76  | | | |
77  |----|-SC-|----|
78  | | | |
79  |____|____|____|
80
81
82  x - direction:
83
84  r = 2 * weight[row][col]*weight[row][col + 1] /
85  (weight[row][col]*weight[row][col + 1]) EC = r * (pot[row][col] - pot[row][col
86  + 1])/dx
87
88  y - direction:
89
90  r = 2 * weight[row][col]*weight[row + 1][col] / (weight[row][col]*weight[row +
91  1][col]) SC = r * (pot[row][col] - pot[row + 1][col])/dy
92
93  the values SC and EC are the values of the next row/col
94
95
96  \endverbatim
97  * \param pot N_array_2d * - the potential N_array_2d
98  * \param weight_x N_array_2d * - the weighting factor N_array_2d used to modify
100  * \param weight_y N_array_2d * - the weighting factor N_array_2d used to modify
102  * \param geom N_geom_data * - geometry data structure
104  size, if a NULL pointer is provided this gradient field will be new allocated
105  * \return N_gradient_field_2d * - the pointer to the computed gradient field
106
107  *
108  * */
110  N_array_2d *weight_x,
111  N_array_2d *weight_y,
112  N_geom_data *geom,
114 {
115  int i, j;
116  int rows, cols;
117  double dx, dy, p1, p2, r1, r2, mean, grad, res;
119
120  if (pot->cols != weight_x->cols || pot->cols != weight_y->cols)
122  "N_compute_gradient_field_2d: the arrays are not of equal size");
123
124  if (pot->rows != weight_x->rows || pot->rows != weight_y->rows)
126  "N_compute_gradient_field_2d: the arrays are not of equal size");
127
128  if (pot->cols != geom->cols || pot->rows != geom->rows)
129  G_fatal_error("N_compute_gradient_field_2d: array sizes and geometry "
130  "data are different");
131
133
134  rows = pot->rows;
135  cols = pot->cols;
136  dx = geom->dx;
137  dy = geom->dy;
138
139  if (field == NULL) {
141  }
142  else {
143  if (field->cols != geom->cols || field->rows != geom->rows)
145  "and geometry data are different");
146  }
147
148  for (j = 0; j < rows; j++)
149  for (i = 0; i < cols - 1; i++) {
151  mean = 0;
152
153  /* Only compute if the arrays are not null */
154  if (!N_is_array_2d_value_null(pot, i, j) &&
155  !N_is_array_2d_value_null(pot, i + 1, j)) {
156  p1 = N_get_array_2d_d_value(pot, i, j);
157  p2 = N_get_array_2d_d_value(pot, i + 1, j);
159  }
160  if (!N_is_array_2d_value_null(weight_x, i, j) &&
161  !N_is_array_2d_value_null(weight_x, i + 1, j)) {
162  r1 = N_get_array_2d_d_value(weight_x, i, j);
163  r2 = N_get_array_2d_d_value(weight_x, i + 1, j);
164  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
165  }
166
167  res = mean * grad;
168
169  N_put_array_2d_d_value(field->x_array, i + 1, j, res);
170  }
171
172  for (j = 0; j < rows - 1; j++)
173  for (i = 0; i < cols; i++) {
175  mean = 0;
176
177  /* Only compute if the arrays are not null */
178  if (!N_is_array_2d_value_null(pot, i, j) &&
179  !N_is_array_2d_value_null(pot, i, j + 1)) {
180  p1 = N_get_array_2d_d_value(pot, i, j);
181  p2 = N_get_array_2d_d_value(pot, i, j + 1);
183  }
184  if (!N_is_array_2d_value_null(weight_y, i, j) &&
185  !N_is_array_2d_value_null(weight_y, i, j + 1)) {
186  r1 = N_get_array_2d_d_value(weight_y, i, j);
187  r2 = N_get_array_2d_d_value(weight_y, i, j + 1);
188  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
189  }
190
191  res = -1 * mean * grad;
192
193  N_put_array_2d_d_value(field->y_array, i, j + 1, res);
194  }
195
196  /*Compute gradient field statistics */
198
199  return field;
200 }
201
202 /*!
203  * \brief Calculate the x and y vector components from a gradient field for each
204  * cell and stores them in the provided N_array_2d structures
205  *
206  * The arrays must have the same size as the gradient field.
207
208  \verbatim
209
210  Based on this storages scheme the gradient vector for each cell is
211  calculated and stored in the provided N_array_2d structures
212
213  ______________
214  | | | |
215  | | | |
216  |----|-NC-|----|
217  | | | |
218  | WC EC |
219  | | | |
220  |----|-SC-|----|
221  | | | |
222  |____|____|____|
223
224  x vector component:
225
226  x = (WC + EC) / 2
227
228  y vector component:
229
230  y = (NC + SC) / 2
231
232  \endverbatim
233  *
234  * \param field N_gradient_field_2d *
235  * \param x_comp N_array_2d * - the array in which the x component will be
236  written
237  * \param y_comp N_array_2d * - the array in which the y component will be
238  written
239  *
240  * \return void
241  * */
243  N_array_2d *x_comp,
244  N_array_2d *y_comp)
245 {
246  int i, j;
247
248  int rows, cols;
249
250  double vx, vy;
251
252  N_array_2d *x = x_comp;
253
254  N_array_2d *y = y_comp;
255
257
258  if (!x)
259  G_fatal_error("N_compute_gradient_components_2d: x array is empty");
260  if (!y)
261  G_fatal_error("N_compute_gradient_components_2d: y array is empty");
262
263  cols = field->x_array->cols;
264  rows = field->x_array->rows;
265
266  /*Check the array sizes */
267  if (x->cols != cols || x->rows != rows)
268  G_fatal_error("N_compute_gradient_components_2d: the size of the x "
269  "array doesn't fit the gradient field size");
270  if (y->cols != cols || y->rows != rows)
271  G_fatal_error("N_compute_gradient_components_2d: the size of the y "
272  "array doesn't fit the gradient field size");
273
274  for (j = 0; j < rows; j++)
275  for (i = 0; i < cols; i++) {
277
278  /* in case a gradient is zero, we expect a no flow boundary */
281  else
285  else
287
288  N_put_array_2d_d_value(x, i, j, vx);
289  N_put_array_2d_d_value(y, i, j, vy);
290  }
291
292  return;
293 }
294
295 /*! \brief Calculate basic statistics of a gradient field
296  *
297  * The statistic is stored in the gradient field struct
298  *
299  * \param field N_gradient_3d_field *
300  * \return void
301  *
302  * */
304 {
305  double minx, miny, minz;
306
307  double maxx, maxy, maxz;
308
309  double sumx, sumy, sumz;
310
311  int nonullx, nonully, nonullz;
312
314
315  N_calc_array_3d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0);
316  N_calc_array_3d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0);
317  N_calc_array_3d_stats(field->z_array, &minz, &maxz, &sumz, &nonullz, 0);
318
319  if (minx <= minz && minx <= miny)
320  field->min = minx;
321  if (miny <= minz && miny <= minx)
322  field->min = miny;
323  if (minz <= minx && minz <= miny)
324  field->min = minz;
325
326  if (maxx >= maxz && maxx >= maxy)
327  field->max = maxx;
328  if (maxy >= maxz && maxy >= maxx)
329  field->max = maxy;
330  if (maxz >= maxx && maxz >= maxy)
331  field->max = maxz;
332
333  field->sum = sumx + sumy + sumz;
334  field->nonull = nonullx + nonully + nonullz;
335  field->mean = field->sum / (double)field->nonull;
336
337  return;
338 }
339
340 /*!
341  * \brief This function computes the gradient based on the input N_array_3d pot
342  * (that means potential), a weighting factor N_array_3d named weight and the
343  distance between two cells
344  * saved in the N_geom_data struct.
345  *
346  * The gradient is calculated between cells for each cell and direction.
347  * An existing gradient field can be filled with new data or, if a NULL pointer
348  is
349  * given, a new gradient field will be allocated with the appropriate size.
350  *
351  *
352  *
353  *
354  \verbatim
355
356  | /
357  TC NC
358  |/
359  --WC-----EC--
360  /|
361  SC BC
362  / |
363
364  x - direction:
365
366  r = 2 * weight_x[depth][row][col]*weight_x[depth][row][col + 1] /
367  (weight_X[depth][row][col]*weight_x[depth][row][col + 1]) EC = r *
368  (pot[depth][row][col] - pot[depth][row][col + 1])/dx
369
370  y - direction:
371
372  r = 2 * weight_y[depth][row][col]*weight_y[depth][row + 1][col] /
373  (weight_y[depth][row][col]*weight_y[depth][row + 1][col]) SC = r *
374  (pot[depth][row][col] - pot[depth][row + 1][col])/dy
375
376  z - direction:
377
378  r = 2 * weight_z[depth][row][col]*weight_z[depth + 1][row][col] /
379  (weight_z[depth][row][col]*weight_z[depth + 1][row][col]) TC = r *
380  (pot[depth][row][col] - pot[depth + 1][row][col])/dy
381
382  the values BC, NC, WC are the values of the next depth/row/col
383
384
385  \endverbatim
386  * \param pot N_array_3d * - the potential N_array_2d
387  * \param weight_x N_array_3d * - the weighting factor N_array_3d used to modify
389  * \param weight_y N_array_3d * - the weighting factor N_array_3d used to modify
391  * \param weight_z N_array_3d * - the weighting factor N_array_3d used to modify
393  * \param geom N_geom_data * - geometry data structure
395  size, if a NULL pointer is provided this gradient field will be new allocated
396  * \return N_gradient_field_3d * - the pointer to the computed gradient field
397  *
398  * */
401  N_array_3d *weight_y, N_array_3d *weight_z,
403 {
404  int i, j, k;
405
406  int cols, rows, depths;
407
408  double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
409
411
412  if (pot->cols != weight_x->cols || pot->cols != weight_y->cols ||
413  pot->cols != weight_z->cols)
415  "N_compute_gradient_field_3d: the arrays are not of equal size");
416
417  if (pot->rows != weight_x->rows || pot->rows != weight_y->rows ||
418  pot->rows != weight_z->rows)
420  "N_compute_gradient_field_3d: the arrays are not of equal size");
421
422  if (pot->depths != weight_x->depths || pot->depths != weight_y->depths ||
423  pot->depths != weight_z->depths)
425  "N_compute_gradient_field_3d: the arrays are not of equal size");
426
427  if (pot->cols != geom->cols || pot->rows != geom->rows ||
428  pot->depths != geom->depths)
429  G_fatal_error("N_compute_gradient_field_3d: array sizes and geometry "
430  "data are different");
431
433
434  cols = geom->cols;
435  rows = geom->rows;
436  depths = geom->depths;
437  dx = geom->dx;
438  dy = geom->dy;
439  dz = geom->dz;
440
441  if (gradfield == NULL) {
442  field = N_alloc_gradient_field_3d(cols, rows, depths);
443  }
444  else {
445  if (field->cols != geom->cols || field->rows != geom->rows ||
446  field->depths != geom->depths)
448  "and geometry data are different");
449  }
450
451  for (k = 0; k < depths; k++)
452  for (j = 0; j < rows; j++)
453  for (i = 0; i < cols - 1; i++) {
455  mean = 0;
456
457  /*Only compute if the arrays are not null */
458  if (!N_is_array_3d_value_null(pot, i, j, k) &&
459  !N_is_array_3d_value_null(pot, i + 1, j, k)) {
460  p1 = N_get_array_3d_d_value(pot, i, j, k);
461  p2 = N_get_array_3d_d_value(pot, i + 1, j, k);
463  }
464  if (!N_is_array_3d_value_null(weight_x, i, j, k) &&
465  !N_is_array_3d_value_null(weight_x, i + 1, j, k)) {
466  r1 = N_get_array_3d_d_value(weight_x, i, j, k);
467  r2 = N_get_array_3d_d_value(weight_x, i + 1, j, k);
468  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
469  }
470
471  res = mean * grad;
472
473  G_debug(6,
474  "N_compute_gradient_field_3d: X-direction insert value "
475  "%6.5g at %i %i %i ",
476  res, k, j, i + 1);
477
478  N_put_array_3d_d_value(field->x_array, i + 1, j, k, res);
479  }
480
481  for (k = 0; k < depths; k++)
482  for (j = 0; j < rows - 1; j++)
483  for (i = 0; i < cols; i++) {
485  mean = 0;
486
487  /* Only compute if the arrays are not null */
488  if (!N_is_array_3d_value_null(pot, i, j, k) &&
489  !N_is_array_3d_value_null(pot, i, j + 1, k)) {
490  p1 = N_get_array_3d_d_value(pot, i, j, k);
491  p2 = N_get_array_3d_d_value(pot, i, j + 1, k);
493  }
494  if (!N_is_array_3d_value_null(weight_y, i, j, k) &&
495  !N_is_array_3d_value_null(weight_y, i, j + 1, k)) {
496  r1 = N_get_array_3d_d_value(weight_y, i, j, k);
497  r2 = N_get_array_3d_d_value(weight_y, i, j + 1, k);
498  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
499  }
500
501  res = -1 * mean * grad; /*invert the direction, because we count
502  * from north to south, but the gradient
503  * is defined in y direction */
504
505  G_debug(6,
506  "N_compute_gradient_field_3d: Y-direction insert value "
507  "%6.5g at %i %i %i ",
508  res, k, j + 1, i);
509
510  N_put_array_3d_d_value(field->y_array, i, j + 1, k, res);
511  }
512
513  for (k = 0; k < depths - 1; k++)
514  for (j = 0; j < rows; j++)
515  for (i = 0; i < cols; i++) {
517  mean = 0;
518
519  /* Only compute if the arrays are not null */
520  if (!N_is_array_3d_value_null(pot, i, j, k) &&
521  !N_is_array_3d_value_null(pot, i, j, k + 1)) {
522  p1 = N_get_array_3d_d_value(pot, i, j, k);
523  p2 = N_get_array_3d_d_value(pot, i, j, k + 1);
525  }
526  if (!N_is_array_3d_value_null(weight_z, i, j, k) &&
527  !N_is_array_3d_value_null(weight_z, i, j, k + 1)) {
528  r1 = N_get_array_3d_d_value(weight_z, i, j, k);
529  r2 = N_get_array_3d_d_value(weight_z, i, j, k + 1);
530  mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */
531  }
532
533  res = mean * grad;
534
535  G_debug(6,
536  "N_compute_gradient_field_3d: Z-direction insert value "
537  "%6.5g at %i %i %i ",
538  res, k + 1, j, i);
539
540  N_put_array_3d_d_value(field->z_array, i, j, k + 1, res);
541  }
542
543  /*Compute gradient field statistics */
545
546  return field;
547 }
548
549 /*!
550  * \brief Calculate the x, y and z vector components from a gradient field for
551  each cell
552  * and store them in the provided N_array_3d structures
553  *
554  * The arrays must have the same size as the gradient field.
555  *
556  \verbatim
557
558  Based on this storages scheme the gradient vector for each cell is
559  calculated and stored in the provided N_array_3d structures
560
561
562  | /
563  TC NC
564  |/
565  --WC-----EC--
566  /|
567  SC BC
568  / |
569
570
571  x vector component:
572
573  x = (WC + EC) / 2
574
575  y vector component:
576
577  y = (NC + SC) / 2
578
579  z vector component:
580
581  z = (TC + BC) / 2
582
583  \endverbatim
584
585  * \param field N_gradient_field_3d *
586  * \param x_comp N_array_3d * - the array in which the x component will be
587  written
588  * \param y_comp N_array_3d * - the array in which the y component will be
589  written
590  * \param z_comp N_array_3d * - the array in which the z component will be
591  written
592  *
593  * \return void
594  * */
596  N_array_3d *x_comp,
597  N_array_3d *y_comp,
598  N_array_3d *z_comp)
599 {
600  int i, j, k;
601
602  int rows, cols, depths;
603
604  double vx, vy, vz;
605
606  N_array_3d *x = x_comp;
607
608  N_array_3d *y = y_comp;
609
610  N_array_3d *z = z_comp;
611
613
614  if (!x)
615  G_fatal_error("N_compute_gradient_components_3d: x array is empty");
616  if (!y)
617  G_fatal_error("N_compute_gradient_components_3d: y array is empty");
618  if (!z)
619  G_fatal_error("N_compute_gradient_components_3d: z array is empty");
620
621  cols = field->x_array->cols;
622  rows = field->x_array->rows;
623  depths = field->x_array->depths;
624
625  /*Check the array sizes */
626  if (x->cols != cols || x->rows != rows || x->depths != depths)
627  G_fatal_error("N_compute_gradient_components_3d: the size of the x "
628  "array doesn't fit the gradient field size");
629  if (y->cols != cols || y->rows != rows || y->depths != depths)
630  G_fatal_error("N_compute_gradient_components_3d: the size of the y "
631  "array doesn't fit the gradient field size");
632  if (z->cols != cols || z->rows != rows || z->depths != depths)
633  G_fatal_error("N_compute_gradient_components_3d: the size of the z "
634  "array doesn't fit the gradient field size");
635
636  for (k = 0; k < depths; k++)
637  for (j = 0; j < rows; j++)
638  for (i = 0; i < cols; i++) {
640  /* in case a gradient is zero, we expect a no flow boundary */
643  else
647  else
651  else
653
654  N_put_array_3d_d_value(x, i, j, k, vx);
655  N_put_array_3d_d_value(y, i, j, k, vy);
656  N_put_array_3d_d_value(z, i, j, k, vz);
657  }
658
659  return;
660 }
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: n_tools.c:115
#define NULL
Definition: ccmath.h:32
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
float mean(IClass_statistics *statistics, int band)
Helper function for computing mean.
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: n_arrays.c:380
int N_is_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function returns 1 if value of N_array_3d data at position col, row, depth is of type null,...
Definition: n_arrays.c:873
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at position col, row is of type null, otherwise 0.
Definition: n_arrays.c:231
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: n_arrays.c:1148
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:979
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:576
void N_calc_array_3d_stats(N_array_3d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_3d struct.
void N_calc_array_2d_stats(N_array_2d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_2d struct.
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col].
Calculate the x, y and z vector components from a gradient field for each cell and store them in the ...
Calculate the x and y vector components from a gradient field for each cell and stores them in the pr...
Calculate basic statistics of a gradient field.
This function computes the gradient based on the input N_array_3d pot (that means potential),...
This function computes the gradient based on the input N_array_2d pot (potential),...
Calculate basic statistics of a gradient field.
int cols
Definition: N_pde.h:134
int rows
Definition: N_pde.h:134
int rows
Definition: N_pde.h:177
int depths
Definition: N_pde.h:177
int cols
Definition: N_pde.h:177
Geometric information about the structured grid.
Definition: N_pde.h:101
double dx
Definition: N_pde.h:108
int depths
Definition: N_pde.h:114
double dy
Definition: N_pde.h:109
double dz
Definition: N_pde.h:110
int rows
Definition: N_pde.h:115
int cols
Definition: N_pde.h:116
Gradient between the cells in X and Y direction.
Definition: N_pde.h:457
double EC
Definition: N_pde.h:459
double SC
Definition: N_pde.h:459
double WC
Definition: N_pde.h:459
double NC
Definition: N_pde.h:459
Gradient between the cells in X, Y and Z direction.
Definition: N_pde.h:464
double WC
Definition: N_pde.h:466
double NC
Definition: N_pde.h:466
double EC
Definition: N_pde.h:466
double SC
Definition: N_pde.h:466
double BC
Definition: N_pde.h:466
double TC
Definition: N_pde.h:466
N_array_2d * y_array
Definition: N_pde.h:562
N_array_2d * x_array
Definition: N_pde.h:561
N_array_3d * y_array
Definition: N_pde.h:573
N_array_3d * x_array
Definition: N_pde.h:572
N_array_3d * z_array
Definition: N_pde.h:574
#define x