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