GRASS GIS 7 Programmer's Manual  7.0.2(2015)-r00000
n_gradient_calc.c
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 
30 {
31  double minx, miny;
32  double maxx, maxy;
33  double sumx, sumy;
34  int nonullx, nonully;
35 
36  G_debug(3,
37  "N_calc_gradient_field_2d_stats: compute gradient field stats");
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 
106  N_array_2d * weight_x,
107  N_array_2d * weight_y,
108  N_geom_data * geom,
110  gradfield)
111 {
112  int i, j;
113  int rows, cols;
114  double dx, dy, p1, p2, r1, r2, mean, grad, res;
115  N_gradient_field_2d *field = gradfield;
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 
131  G_debug(3, "N_compute_gradient_field_2d: compute gradient field");
132 
133  rows = pot->rows;
134  cols = pot->cols;
135  dx = geom->dx;
136  dy = geom->dy;
137 
138  if (field == NULL) {
139  field = N_alloc_gradient_field_2d(cols, rows);
140  }
141  else {
142  if (field->cols != geom->cols || field->rows != geom->rows)
144  ("N_compute_gradient_field_2d: gradient field sizes and geometry data are different");
145  }
146 
147 
148  for (j = 0; j < rows; j++)
149  for (i = 0; i < cols - 1; i++) {
150  grad = 0;
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);
158  grad = (p1 - p2) / dx; /* gradient */
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++) {
175  grad = 0;
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);
183  grad = (p1 - p2) / dy; /* gradient */
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 
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 
257  N_gradient_2d grad;
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++) {
278  N_get_gradient_2d(field, &grad, i, j);
279 
280  /* in case a gradient is zero, we expect a no flow boundary */
281  if (grad.WC == 0.0 || grad.EC == 0.0)
282  vx = (grad.WC + grad.EC);
283  else
284  vx = (grad.WC + grad.EC) / 2;
285  if (grad.NC == 0.0 || grad.SC == 0.0)
286  vy = (grad.NC + grad.SC);
287  else
288  vy = (grad.NC + grad.SC) / 2;
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 
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,
316  "N_calc_gradient_field_3d_stats: compute gradient field stats");
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 
395  N_array_3d * weight_x,
396  N_array_3d * weight_y,
397  N_array_3d * weight_z,
398  N_geom_data * geom,
400  gradfield)
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 
408  N_gradient_field_3d *field = gradfield;
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 
431  G_debug(3, "N_compute_gradient_field_3d: compute gradient field");
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)
447  ("N_compute_gradient_field_3d: gradient field sizes and geometry data are different");
448  }
449 
450  for (k = 0; k < depths; k++)
451  for (j = 0; j < rows; j++)
452  for (i = 0; i < cols - 1; i++) {
453  grad = 0;
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);
461  grad = (p1 - p2) / dx; /* gradient */
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++) {
483  grad = 0;
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);
491  grad = (p1 - p2) / dy; /* gradient */
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++) {
514  grad = 0;
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);
522  grad = (p1 - p2) / dz; /* gradient */
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 
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 
607  N_gradient_3d grad;
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++) {
635  N_get_gradient_3d(field, &grad, i, j, k);
636  /* in case a gradient is zero, we expect a no flow boundary */
637  if (grad.WC == 0.0 || grad.EC == 0.0)
638  vx = (grad.WC + grad.EC);
639  else
640  vx = (grad.WC + grad.EC) / 2;
641  if (grad.NC == 0.0 || grad.SC == 0.0)
642  vy = (grad.NC + grad.SC);
643  else
644  vy = (grad.NC + grad.SC) / 2;
645  if (grad.TC == 0.0 || grad.BC == 0.0)
646  vz = (grad.TC + grad.BC);
647  else
648  vz = (grad.TC + grad.BC) / 2;
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 }
N_gradient_2d * N_get_gradient_2d(N_gradient_field_2d *field, N_gradient_2d *gradient, int col, int row)
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col]...
Definition: n_gradient.c:115
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
N_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
Definition: n_gradient.c:1018
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
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:159
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
void N_calc_gradient_field_2d_stats(N_gradient_field_2d *field)
Calculate basic statistics of a gradient field.
int rows
Definition: N_pde.h:116
double TC
Definition: N_pde.h:419
N_gradient_field_3d * N_compute_gradient_field_3d(N_array_3d *pot, N_array_3d *weight_x, N_array_3d *weight_y, N_array_3d *weight_z, N_geom_data *geom, N_gradient_field_3d *gradfield)
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 postion col, row is of type null, otherwise 0...
Definition: n_arrays.c:231
int cols
Definition: N_pde.h:117
Gradient between the cells in X, Y and Z direction.
Definition: N_pde.h:416
void N_calc_gradient_field_3d_stats(N_gradient_field_3d *field)
Calculate basic statistics of a gradient field.
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition: n_gradient.c:920
void N_compute_gradient_field_components_2d(N_gradient_field_2d *field, N_array_2d *x_comp, N_array_2d *y_comp)
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
N_gradient_field_2d * N_compute_gradient_field_2d(N_array_2d *pot, N_array_2d *weight_x, N_array_2d *weight_y, N_geom_data *geom, N_gradient_field_2d *gradfield)
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
void N_compute_gradient_field_components_3d(N_gradient_field_3d *field, N_array_3d *x_comp, N_array_3d *y_comp, N_array_3d *z_comp)
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
double NC
Definition: N_pde.h:411
N_array_2d * x_array
Definition: N_pde.h:523
int rows
Definition: N_pde.h:168
N_gradient_3d * N_get_gradient_3d(N_gradient_field_3d *field, N_gradient_3d *gradient, int col, int row, int depth)
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
Definition: n_gradient.c:248
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:881