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