GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r73111
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
n_arrays_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: Higher level array 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 #include <math.h>
19 
20 #include <grass/N_pde.h>
21 #include <grass/raster.h>
22 #include <grass/glocale.h>
23 
24 
25 /* ******************** 2D ARRAY FUNCTIONS *********************** */
26 
27 /*!
28  * \brief Copy the source N_array_2d struct to the target N_array_2d struct
29  *
30  * The arrays must have the same size and the same offset.
31  *
32  * The array types can be mixed, the values are automatically casted
33  * and the null values are set accordingly.
34  * <br><br>
35  * If you copy a cell array into a dcell array, the values are casted to dcell and
36  * the null values are converted from cell-null to dcell-null
37  * <br><br>
38  * This function can be called in a parallel region defined with OpenMP.
39  * The copy loop is parallelize with a openmp for pragma.
40  *
41  * \param source N_array_2d *
42  * \param target N_array_2d *
43  * \return void
44  * */
46 {
47  int i;
48  int null = 0;
49 
50 #pragma omp single
51  {
52  if (source->cols_intern != target->cols_intern)
54  ("N_copy_array_2d: the arrays are not of equal size");
55 
56  if (source->rows_intern != target->rows_intern)
58  ("N_copy_array_2d: the arrays are not of equal size");
59 
60  G_debug(3,
61  "N_copy_array_2d: copy source array to target array size %i",
62  source->cols_intern * source->rows_intern);
63  }
64 
65 #pragma omp for
66  for (i = 0; i < source->cols_intern * source->rows_intern; i++) {
67  null = 0;
68  if (source->type == CELL_TYPE) {
69  if (Rast_is_c_null_value((void *)&source->cell_array[i]))
70  null = 1;
71 
72  if (target->type == CELL_TYPE) {
73  target->cell_array[i] = source->cell_array[i];
74  }
75  if (target->type == FCELL_TYPE) {
76  if (null)
77  Rast_set_f_null_value((void *)&(target->fcell_array[i]), 1);
78  else
79  target->fcell_array[i] = (FCELL) source->cell_array[i];
80  }
81  if (target->type == DCELL_TYPE) {
82  if (null)
83  Rast_set_d_null_value((void *)&(target->dcell_array[i]), 1);
84  else
85  target->dcell_array[i] = (DCELL) source->cell_array[i];
86  }
87 
88  }
89  if (source->type == FCELL_TYPE) {
90  if (Rast_is_f_null_value((void *)&source->fcell_array[i]))
91  null = 1;
92 
93  if (target->type == CELL_TYPE) {
94  if (null)
95  Rast_set_c_null_value((void *)&(target->cell_array[i]), 1);
96  else
97  target->cell_array[i] = (CELL) source->fcell_array[i];
98  }
99  if (target->type == FCELL_TYPE) {
100  target->fcell_array[i] = source->fcell_array[i];
101  }
102  if (target->type == DCELL_TYPE) {
103  if (null)
104  Rast_set_d_null_value((void *)&(target->dcell_array[i]), 1);
105  else
106  target->dcell_array[i] = (DCELL) source->fcell_array[i];
107  }
108  }
109  if (source->type == DCELL_TYPE) {
110  if (Rast_is_d_null_value((void *)&source->dcell_array[i]))
111  null = 1;
112 
113  if (target->type == CELL_TYPE) {
114  if (null)
115  Rast_set_c_null_value((void *)&(target->cell_array[i]), 1);
116  else
117  target->cell_array[i] = (CELL) source->dcell_array[i];
118  }
119  if (target->type == FCELL_TYPE) {
120  if (null)
121  Rast_set_f_null_value((void *)&(target->fcell_array[i]), 1);
122  else
123  target->fcell_array[i] = (FCELL) source->dcell_array[i];
124  }
125  if (target->type == DCELL_TYPE) {
126  target->dcell_array[i] = source->dcell_array[i];
127  }
128  }
129  }
130 
131  return;
132 }
133 
134 /*!
135  * \brief Calculate the norm of the two input arrays
136  *
137  * The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
138  * All arrays must have equal sizes and offsets.
139  * The complete data array inclusively offsets is used for norm calucaltion.
140  * Only non-null values are used to calculate the norm.
141  *
142 
143  * \param a N_array_2d *
144  * \param b N_array_2d *
145  * \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
146  * \return double the calculated norm
147  * */
148 double N_norm_array_2d(N_array_2d * a, N_array_2d * b, int type)
149 {
150  int i = 0;
151  double norm = 0.0, tmp = 0.0;
152  double v1 = 0.0, v2 = 0.0;
153 
154  if (a->cols_intern != b->cols_intern)
155  G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
156 
157  if (a->rows_intern != b->rows_intern)
158  G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
159 
160  G_debug(3, "N_norm_array_2d: norm of a and b size %i",
161  a->cols_intern * a->rows_intern);
162 
163  for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
164  v1 = 0.0;
165  v2 = 0.0;
166 
167  if (a->type == CELL_TYPE) {
168  if (!Rast_is_f_null_value((void *)&(a->cell_array[i])))
169  v1 = (double)a->cell_array[i];
170  }
171  if (a->type == FCELL_TYPE) {
172  if (!Rast_is_f_null_value((void *)&(a->fcell_array[i])))
173  v1 = (double)a->fcell_array[i];
174  }
175  if (a->type == DCELL_TYPE) {
176  if (!Rast_is_f_null_value((void *)&(a->dcell_array[i])))
177  v1 = (double)a->dcell_array[i];
178  }
179  if (b->type == CELL_TYPE) {
180  if (!Rast_is_f_null_value((void *)&(b->cell_array[i])))
181  v2 = (double)b->cell_array[i];
182  }
183  if (b->type == FCELL_TYPE) {
184  if (!Rast_is_f_null_value((void *)&(b->fcell_array[i])))
185  v2 = (double)b->fcell_array[i];
186  }
187  if (b->type == DCELL_TYPE) {
188  if (!Rast_is_f_null_value((void *)&(b->dcell_array[i])))
189  v2 = (double)b->dcell_array[i];
190  }
191 
192  if (type == N_MAXIMUM_NORM) {
193  tmp = fabs(v2 - v1);
194  if ((tmp > norm))
195  norm = tmp;
196  }
197  if (type == N_EUKLID_NORM) {
198  norm += fabs(v2 - v1);
199  }
200  }
201 
202  return norm;
203 }
204 
205 /*!
206  * \brief Calculate basic statistics of the N_array_2d struct
207  *
208  * Calculates the minimum, maximum, sum and the number of
209  * non null values. The array offset can be included in the calculation.
210  *
211  * \param a N_array_2d * - input array
212  * \param min double* - variable to store the computed minimum
213  * \param max double* - variable to store the computed maximum
214  * \param sum double* - variable to store the computed sum
215  * \param nonull int* - variable to store the number of non null values
216  * \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
217  * \return void
218  * */
219 void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
220  double *sum, int *nonull, int withoffset)
221 {
222  int i, j;
223  double val;
224 
225  *sum = 0.0;
226  *nonull = 0;
227 
228  if (withoffset == 1) {
229 
230  *min =
231  (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
232  *max =
233  (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
234 
235  for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
236  for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
237  if (!N_is_array_2d_value_null(a, i, j)) {
238  val = (double)N_get_array_2d_d_value(a, i, j);
239  if (*min > val)
240  *min = val;
241  if (*max < val)
242  *max = val;
243  *sum += val;
244  (*nonull)++;
245  }
246  }
247  }
248  }
249  else {
250 
251  *min = (double)N_get_array_2d_d_value(a, 0, 0);
252  *max = (double)N_get_array_2d_d_value(a, 0, 0);
253 
254 
255  for (j = 0; j < a->rows; j++) {
256  for (i = 0; i < a->cols; i++) {
257  if (!N_is_array_2d_value_null(a, i, j)) {
258  val = (double)N_get_array_2d_d_value(a, i, j);
259  if (*min > val)
260  *min = val;
261  if (*max < val)
262  *max = val;
263  *sum += val;
264  (*nonull)++;
265  }
266  }
267  }
268  }
269 
270  G_debug(3,
271  "N_calc_array_2d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
272  *min, *max, *sum, *nonull);
273  return;
274 }
275 
276 
277 /*!
278  * \brief Perform calculations with two input arrays,
279  * the result is written to a third array.
280  *
281  * All arrays must have equal sizes and offsets.
282  * The complete data array inclusively offsets is used for calucaltions.
283  * Only non-null values are computed. If one array value is null,
284  * the result array value will be null too.
285  * <br><br>
286  * If a division with zero is detected, the resulting arrays
287  * value will set to null and not to NaN.
288  * <br><br>
289  * The result array is optional, if the result arrays points to NULL,
290  * a new array will be allocated with the largest arrays data type
291  * (CELL, FCELL or DCELL) used by the input arrays.
292  * <br><br>
293  * the array computations can be of the following forms:
294  *
295  * <ul>
296  * <li>result = a + b -> N_ARRAY_SUM</li>
297  * <li>result = a - b -> N_ARRAY_DIF</li>
298  * <li>result = a * b -> N_ARRAY_MUL</li>
299  * <li>result = a / b -> N_ARRAY_DIV</li>
300  * </ul>
301  *
302  * \param a N_array_2d * - first input array
303  * \param b N_array_2d * - second input array
304  * \param result N_array_2d * - the optional result array
305  * \param type - the type of calculation
306  * \return N_array_2d * - the pointer to the result array
307  * */
309  N_array_2d * result, int type)
310 {
311  N_array_2d *c;
312  int i, j, setnull = 0;
313  double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
314 
315  /*Set the pointer */
316  c = result;
317 
318 #pragma omp single
319  {
320  /*Check the array sizes */
321  if (a->cols_intern != b->cols_intern)
323  ("N_math_array_2d: the arrays are not of equal size");
324  if (a->rows_intern != b->rows_intern)
326  ("N_math_array_2d: the arrays are not of equal size");
327  if (a->offset != b->offset)
329  ("N_math_array_2d: the arrays have different offsets");
330 
331  G_debug(3, "N_math_array_2d: mathematical calculations, size: %i",
332  a->cols_intern * a->rows_intern);
333 
334  /*if the result array is null, allocate a new one, use the
335  * largest data type of the input arrays*/
336  if (c == NULL) {
337  if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
338  c = N_alloc_array_2d(a->cols, a->rows, a->offset, DCELL_TYPE);
339  G_debug(3,
340  "N_math_array_2d: array of type DCELL_TYPE created");
341  }
342  else if (a->type == FCELL_TYPE || b->type == FCELL_TYPE) {
343  c = N_alloc_array_2d(a->cols, a->rows, a->offset, FCELL_TYPE);
344  G_debug(3,
345  "N_math_array_2d: array of type FCELL_TYPE created");
346  }
347  else {
348  c = N_alloc_array_2d(a->cols, a->rows, a->offset, CELL_TYPE);
349  G_debug(3,
350  "N_math_array_2d: array of type CELL_TYPE created");
351  }
352  }
353  else {
354  /*Check the array sizes */
355  if (a->cols_intern != c->cols_intern)
357  ("N_math_array_2d: the arrays are not of equal size");
358  if (a->rows_intern != c->rows_intern)
360  ("N_math_array_2d: the arrays are not of equal size");
361  if (a->offset != c->offset)
363  ("N_math_array_2d: the arrays have different offsets");
364  }
365  }
366 
367 #pragma omp for private(va, vb, vc, setnull)
368  for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
369  for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
370  if (!N_is_array_2d_value_null(a, i, j) &&
371  !N_is_array_2d_value_null(b, i, j)) {
372  /*we always calculate internally with double values */
373  va = (double)N_get_array_2d_d_value(a, i, j);
374  vb = (double)N_get_array_2d_d_value(b, i, j);
375  vc = 0;
376  setnull = 0;
377 
378  switch (type) {
379  case N_ARRAY_SUM:
380  vc = va + vb;
381  break;
382  case N_ARRAY_DIF:
383  vc = va - vb;
384  break;
385  case N_ARRAY_MUL:
386  vc = va * vb;
387  break;
388  case N_ARRAY_DIV:
389  if (vb != 0)
390  vc = va / vb;
391  else
392  setnull = 1;
393  break;
394  }
395 
396  if (c->type == CELL_TYPE) {
397  if (setnull)
398  N_put_array_2d_value_null(c, i, j);
399  else
400  N_put_array_2d_c_value(c, i, j, (CELL) vc);
401  }
402  if (c->type == FCELL_TYPE) {
403  if (setnull)
404  N_put_array_2d_value_null(c, i, j);
405  else
406  N_put_array_2d_f_value(c, i, j, (FCELL) vc);
407  }
408  if (c->type == DCELL_TYPE) {
409  if (setnull)
410  N_put_array_2d_value_null(c, i, j);
411  else
412  N_put_array_2d_d_value(c, i, j, (DCELL) vc);
413  }
414 
415  }
416  else {
417  N_put_array_2d_value_null(c, i, j);
418  }
419  }
420  }
421 
422  return c;
423 }
424 
425 /*!
426  * \brief Convert all null values to zero values
427  *
428  * The complete data array inclusively offsets is used.
429  * The array data types are automatically recognized.
430  *
431  * \param a N_array_2d *
432  * \return int - number of replaced values
433  * */
435 {
436  int i = 0, count = 0;
437 
438  G_debug(3, "N_convert_array_2d_null_to_zero: convert array of size %i",
439  a->cols_intern * a->rows_intern);
440 
441  if (a->type == CELL_TYPE)
442  for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
443  if (Rast_is_c_null_value((void *)&(a->cell_array[i]))) {
444  a->cell_array[i] = 0;
445  count++;
446  }
447  }
448 
449  if (a->type == FCELL_TYPE)
450  for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
451  if (Rast_is_f_null_value((void *)&(a->fcell_array[i]))) {
452  a->fcell_array[i] = 0.0;
453  count++;
454  }
455  }
456 
457 
458  if (a->type == DCELL_TYPE)
459  for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
460  if (Rast_is_d_null_value((void *)&(a->dcell_array[i]))) {
461  a->dcell_array[i] = 0.0;
462  count++;
463  }
464  }
465 
466 
467  if (a->type == CELL_TYPE)
468  G_debug(2,
469  "N_convert_array_2d_null_to_zero: %i values of type CELL_TYPE are converted",
470  count);
471  if (a->type == FCELL_TYPE)
472  G_debug(2,
473  "N_convert_array_2d_null_to_zero: %i valuess of type FCELL_TYPE are converted",
474  count);
475  if (a->type == DCELL_TYPE)
476  G_debug(2,
477  "N_convert_array_2d_null_to_zero: %i valuess of type DCELL_TYPE are converted",
478  count);
479 
480  return count;
481 }
482 
483 /* ******************** 3D ARRAY FUNCTIONS *********************** */
484 
485 /*!
486  * \brief Copy the source N_array_3d struct to the target N_array_3d struct
487  *
488  * The arrays must have the same size and the same offset.
489  *
490  * The array data types can be mixed, the values are automatically casted
491  * and the null values are set accordingly.
492  *
493  * If you copy a float array to a double array, the values are casted to DCELL and
494  * the null values are converted from FCELL-null to DCELL-null
495  *
496  * \param source N_array_3d *
497  * \param target N_array_3d *
498  * \return void
499  * */
501 {
502  int i;
503  int null;
504 
505  if (source->cols_intern != target->cols_intern)
506  G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
507 
508  if (source->rows_intern != target->rows_intern)
509  G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
510 
511  if (source->depths_intern != target->depths_intern)
512  G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
513 
514 
515  G_debug(3, "N_copy_array_3d: copy source array to target array size %i",
516  source->cols_intern * source->rows_intern *
517  source->depths_intern);
518 
519  for (i = 0;
520  i <
521  source->cols_intern * source->rows_intern * source->depths_intern;
522  i++) {
523  null = 0;
524  if (source->type == FCELL_TYPE) {
526  ((void *)&(source->fcell_array[i]), FCELL_TYPE))
527  null = 1;
528 
529  if (target->type == FCELL_TYPE) {
530  target->fcell_array[i] = source->fcell_array[i];
531  }
532  if (target->type == DCELL_TYPE) {
533  if (null)
534  Rast3d_set_null_value((void *)&(target->dcell_array[i]), 1,
535  DCELL_TYPE);
536  else
537  target->dcell_array[i] = (double)source->fcell_array[i];
538  }
539 
540  }
541  if (source->type == DCELL_TYPE) {
543  ((void *)&(source->dcell_array[i]), DCELL_TYPE))
544  null = 1;
545 
546  if (target->type == FCELL_TYPE) {
547  if (null)
548  Rast3d_set_null_value((void *)&(target->fcell_array[i]), 1,
549  FCELL_TYPE);
550  else
551  target->fcell_array[i] = (float)source->dcell_array[i];
552  }
553  if (target->type == DCELL_TYPE) {
554  target->dcell_array[i] = source->dcell_array[i];
555  }
556  }
557  }
558 
559  return;
560 }
561 
562 
563 /*!
564  * \brief Calculate the norm of the two input arrays
565  *
566  * The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
567  * All arrays must have equal sizes and offsets.
568  * The complete data array inclusively offsets is used for norm calucaltion.
569  * Only non-null values are used to calculate the norm.
570  *
571  * \param a N_array_3d *
572  * \param b N_array_3d *
573  * \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
574  * \return double the calculated norm
575  * */
576 double N_norm_array_3d(N_array_3d * a, N_array_3d * b, int type)
577 {
578  int i = 0;
579  double norm = 0.0, tmp = 0.0;
580  double v1 = 0.0, v2 = 0.0;
581 
582  if (a->cols_intern != b->cols_intern)
583  G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
584 
585  if (a->rows_intern != b->rows_intern)
586  G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
587 
588  if (a->depths_intern != b->depths_intern)
589  G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
590 
591  G_debug(3, "N_norm_array_3d: norm of a and b size %i",
592  a->cols_intern * a->rows_intern * a->depths_intern);
593 
594  for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern; i++) {
595  v1 = 0.0;
596  v2 = 0.0;
597 
598  if (a->type == FCELL_TYPE) {
599  if (!Rast3d_is_null_value_num((void *)&(a->fcell_array[i]), FCELL_TYPE))
600  v1 = (double)a->fcell_array[i];
601  }
602  if (a->type == DCELL_TYPE) {
603  if (!Rast3d_is_null_value_num((void *)&(a->dcell_array[i]), DCELL_TYPE))
604  v1 = (double)a->dcell_array[i];
605  }
606  if (b->type == FCELL_TYPE) {
607  if (!Rast3d_is_null_value_num((void *)&(b->fcell_array[i]), FCELL_TYPE))
608  v2 = (double)b->fcell_array[i];
609  }
610  if (b->type == DCELL_TYPE) {
611  if (!Rast3d_is_null_value_num((void *)&(b->dcell_array[i]), DCELL_TYPE))
612  v2 = (double)b->dcell_array[i];
613  }
614 
615  if (type == N_MAXIMUM_NORM) {
616  tmp = fabs(v2 - v1);
617  if ((tmp > norm))
618  norm = tmp;
619  }
620  if (type == N_EUKLID_NORM) {
621  norm += fabs(v2 - v1);
622  }
623  }
624 
625  return norm;
626 }
627 
628 /*!
629  * \brief Calculate basic statistics of the N_array_3d struct
630  *
631  * Calculates the minimum, maximum, sum and the number of
632  * non null values. The array offset can be included in the statistical calculation.
633  *
634  * \param a N_array_3d * - input array
635  * \param min double* - variable to store the computed minimum
636  * \param max double* - variable to store the computed maximum
637  * \param sum double* - variable to store the computed sum
638  * \param nonull int* - variable to store the number of non null values
639  * \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
640  * \return void
641  * */
642 void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max,
643  double *sum, int *nonull, int withoffset)
644 {
645  int i, j, k;
646  double val;
647 
648  *sum = 0.0;
649  *nonull = 0;
650 
651  if (withoffset == 1) {
652 
653  *min =
654  (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
655  0 - a->offset);
656  *max =
657  (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
658  0 - a->offset);
659 
660  for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
661  for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
662  for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
663  if (!N_is_array_3d_value_null(a, i, j, k)) {
664  val = (double)N_get_array_3d_d_value(a, i, j, k);
665  if (*min > val)
666  *min = val;
667  if (*max < val)
668  *max = val;
669  *sum += val;
670  (*nonull)++;
671  }
672  }
673  }
674  }
675  }
676  else {
677 
678  *min = (double)N_get_array_3d_d_value(a, 0, 0, 0);
679  *max = (double)N_get_array_3d_d_value(a, 0, 0, 0);
680 
681  for (k = 0; k < a->depths; k++) {
682  for (j = 0; j < a->rows; j++) {
683  for (i = 0; i < a->cols; i++) {
684  if (!N_is_array_3d_value_null(a, i, j, k)) {
685  val = (double)N_get_array_3d_d_value(a, i, j, k);
686  if (*min > val)
687  *min = val;
688  if (*max < val)
689  *max = val;
690  *sum += val;
691  (*nonull)++;
692  }
693  }
694  }
695  }
696  }
697 
698  G_debug(3,
699  "N_calc_array_3d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
700  *min, *max, *sum, *nonull);
701 
702  return;
703 }
704 
705 /*!
706  * \brief Perform calculations with two input arrays,
707  * the result is written to a third array.
708  *
709  * All arrays must have equal sizes and offsets.
710  * The complete data array inclusively offsets is used for calucaltions.
711  * Only non-null values are used. If one array value is null,
712  * the result array value will be null too.
713  * <br><br>
714  *
715  * If a division with zero is detected, the resulting arrays
716  * value will set to null and not to NaN.
717  * <br><br>
718  *
719  * The result array is optional, if the result arrays points to NULL,
720  * a new array will be allocated with the largest arrays data type
721  * (FCELL_TYPE or DCELL_TYPE) used by the input arrays.
722  * <br><br>
723  *
724  * the calculations are of the following form:
725  *
726  * <ul>
727  * <li>result = a + b -> N_ARRAY_SUM</li>
728  * <li>result = a - b -> N_ARRAY_DIF</li>
729  * <li>result = a * b -> N_ARRAY_MUL</li>
730  * <li>result = a / b -> N_ARRAY_DIV</li>
731  * </ul>
732  *
733  * \param a N_array_3d * - first input array
734  * \param b N_array_3d * - second input array
735  * \param result N_array_3d * - the optional result array
736  * \param type - the type of calculation
737  * \return N_array_3d * - the pointer to the result array
738  * */
740  N_array_3d * result, int type)
741 {
742  N_array_3d *c;
743  int i, j, k, setnull = 0;
744  double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
745 
746  /*Set the pointer */
747  c = result;
748 
749  /*Check the array sizes */
750  if (a->cols_intern != b->cols_intern)
751  G_fatal_error("N_math_array_3d: the arrays are not of equal size");
752  if (a->rows_intern != b->rows_intern)
753  G_fatal_error("N_math_array_3d: the arrays are not of equal size");
754  if (a->depths_intern != b->depths_intern)
755  G_fatal_error("N_math_array_3d: the arrays are not of equal size");
756  if (a->offset != b->offset)
757  G_fatal_error("N_math_array_3d: the arrays have different offsets");
758 
759  G_debug(3, "N_math_array_3d: mathematical calculations, size: %i",
760  a->cols_intern * a->rows_intern * a->depths_intern);
761 
762  /*if the result array is null, allocate a new one, use the
763  * largest data type of the input arrays*/
764  if (c == NULL) {
765  if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
766  c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
767  DCELL_TYPE);
768  G_debug(3, "N_math_array_3d: array of type DCELL_TYPE created");
769  }
770  else {
771  c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
772  FCELL_TYPE);
773  G_debug(3, "N_math_array_3d: array of type FCELL_TYPE created");
774  }
775  }
776  else {
777  /*Check the array sizes */
778  if (a->cols_intern != c->cols_intern)
780  ("N_math_array_3d: the arrays are not of equal size");
781  if (a->rows_intern != c->rows_intern)
783  ("N_math_array_3d: the arrays are not of equal size");
784  if (a->depths_intern != c->depths_intern)
786  ("N_math_array_3d: the arrays are not of equal size");
787  if (a->offset != c->offset)
789  ("N_math_array_3d: the arrays have different offsets");
790  }
791 
792  for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
793  for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
794  for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
795  if (!N_is_array_3d_value_null(a, i, j, k) &&
796  !N_is_array_3d_value_null(a, i, j, k)) {
797  /*we always calculate internally with double values */
798  va = (double)N_get_array_3d_d_value(a, i, j, k);
799  vb = (double)N_get_array_3d_d_value(b, i, j, k);
800  vc = 0;
801  setnull = 0;
802 
803  switch (type) {
804  case N_ARRAY_SUM:
805  vc = va + vb;
806  break;
807  case N_ARRAY_DIF:
808  vc = va - vb;
809  break;
810  case N_ARRAY_MUL:
811  vc = va * vb;
812  break;
813  case N_ARRAY_DIV:
814  if (vb != 0)
815  vc = va / vb;
816  else
817  setnull = 1;
818  break;
819  }
820 
821  if (c->type == FCELL_TYPE) {
822  if (setnull)
823  N_put_array_3d_value_null(c, i, j, k);
824  else
825  N_put_array_3d_f_value(c, i, j, k, (float)vc);
826  }
827  if (c->type == DCELL_TYPE) {
828  if (setnull)
829  N_put_array_3d_value_null(c, i, j, k);
830  else
831  N_put_array_3d_d_value(c, i, j, k, vc);
832  }
833  }
834  else {
835  N_put_array_3d_value_null(c, i, j, k);
836  }
837  }
838  }
839  }
840 
841  return c;
842 }
843 
844 /*!
845  * \brief Convert all null values to zero values
846  *
847  * The complete data array inclusively offsets is used.
848  *
849  * \param a N_array_3d *
850  * \return int - number of replaced null values
851  * */
853 {
854  int i = 0, count = 0;
855 
856  G_debug(3, "N_convert_array_3d_null_to_zero: convert array of size %i",
857  a->cols_intern * a->rows_intern * a->depths_intern);
858 
859  if (a->type == FCELL_TYPE)
860  for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
861  i++) {
862  if (Rast3d_is_null_value_num((void *)&(a->fcell_array[i]), FCELL_TYPE)) {
863  a->fcell_array[i] = 0.0;
864  count++;
865  }
866  }
867 
868  if (a->type == DCELL_TYPE)
869  for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
870  i++) {
871  if (Rast3d_is_null_value_num((void *)&(a->dcell_array[i]), DCELL_TYPE)) {
872  a->dcell_array[i] = 0.0;
873  count++;
874  }
875  }
876 
877 
878  if (a->type == FCELL_TYPE)
879  G_debug(3,
880  "N_convert_array_3d_null_to_zero: %i values of type FCELL_TYPE are converted",
881  count);
882 
883  if (a->type == DCELL_TYPE)
884  G_debug(3,
885  "N_convert_array_3d_null_to_zero: %i values of type DCELL_TYPE are converted",
886  count);
887 
888  return count;
889 }
#define CELL_TYPE
Definition: raster.h:11
void N_put_array_2d_value_null(N_array_2d *data, int col, int row)
Writes the null value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:458
#define N_MAXIMUM_NORM
Definition: N_pde.h:45
int depths_intern
Definition: N_pde.h:169
int rows
Definition: N_pde.h:134
int Rast_is_f_null_value(const FCELL *fcellVal)
To check if a FCELL raster value is set to NULL.
Definition: null_val.c:242
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.
void N_copy_array_3d(N_array_3d *source, N_array_3d *target)
Copy the source N_array_3d struct to the target N_array_3d 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
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: n_arrays.c:726
#define min(x, y)
Definition: draw2.c:31
void N_put_array_3d_f_value(N_array_3d *data, int col, int row, int depth, float value)
This function writes a float value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1148
double DCELL
Definition: gis.h:581
int depths
Definition: N_pde.h:168
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: n_arrays.c:72
int count
int cols
Definition: N_pde.h:134
int Rast_is_c_null_value(const CELL *cellVal)
To check if a CELL raster value is set to NULL.
Definition: null_val.c:209
int N_convert_array_3d_null_to_zero(N_array_3d *a)
Convert all null values to zero values.
#define NULL
Definition: ccmath.h:32
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:524
#define max(x, y)
Definition: draw2.c:32
void N_put_array_2d_f_value(N_array_2d *data, int col, int row, FCELL value)
Writes a FCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:554
int offset
Definition: N_pde.h:170
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
int offset
Definition: N_pde.h:136
void N_put_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function writes a null value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1076
int N_convert_array_2d_null_to_zero(N_array_2d *a)
Convert all null values to zero values.
CELL * cell_array
Definition: N_pde.h:137
void N_copy_array_2d(N_array_2d *source, N_array_2d *target)
Copy the source N_array_2d struct to the target N_array_2d struct.
Definition: n_arrays_calc.c:45
N_array_2d * N_math_array_2d(N_array_2d *a, N_array_2d *b, N_array_2d *result, int type)
Perform calculations with two input arrays, the result is written to a third array.
int cols_intern
Definition: N_pde.h:169
#define DCELL_TYPE
Definition: raster.h:13
void Rast_set_c_null_value(CELL *cellVals, int numVals)
To set a number of CELL raster values to NULL.
Definition: null_val.c:124
double b
Definition: r_raster.c:39
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
FCELL * fcell_array
Definition: N_pde.h:138
int Rast3d_is_null_value_num(const void *n, int type)
Definition: null.c:12
DCELL * dcell_array
Definition: N_pde.h:139
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
double N_norm_array_3d(N_array_3d *a, N_array_3d *b, int type)
Calculate the norm of the two input arrays.
void Rast_set_f_null_value(FCELL *fcellVals, int numVals)
To set a number of FCELL raster values to NULL.
Definition: null_val.c:138
const char * source
Definition: lz4.h:575
double N_norm_array_2d(N_array_2d *a, N_array_2d *b, int type)
Calculate the norm of the two input arrays.
float * fcell_array
Definition: N_pde.h:171
void Rast3d_set_null_value(void *c, int nofElts, int type)
Fills the vector pointed to by c with nofElts NULL-values of type.
Definition: null.c:35
void Rast_set_d_null_value(DCELL *dcellVals, int numVals)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:155
float FCELL
Definition: gis.h:582
double * dcell_array
Definition: N_pde.h:172
int type
Definition: N_pde.h:133
int CELL
Definition: gis.h:580
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
int cols_intern
Definition: N_pde.h:135
#define FCELL_TYPE
Definition: raster.h:12
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.
#define N_ARRAY_DIF
Definition: N_pde.h:49
N_array_3d * N_math_array_3d(N_array_3d *a, N_array_3d *b, N_array_3d *result, int type)
Perform calculations with two input arrays, the result is written to a third array.
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
#define N_EUKLID_NORM
Definition: N_pde.h:46
int type
Definition: N_pde.h:167
int rows_intern
Definition: N_pde.h:169
#define N_ARRAY_MUL
Definition: N_pde.h:50
int rows
Definition: N_pde.h:168
#define N_ARRAY_DIV
Definition: N_pde.h:51
int rows_intern
Definition: N_pde.h:135
#define N_ARRAY_SUM
Definition: N_pde.h:48
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
int Rast_is_d_null_value(const DCELL *dcellVal)
To check if a DCELL raster value is set to NULL.
Definition: null_val.c:261