1 #include <math.h>
2
3 #include <grass/gis.h>
4 #include <grass/raster.h>
5
6 #include "local_proto.h"
7
add_row_area(DCELL * top,DCELL * bottom,double sz,struct Cell_head * w,double * low,double * high)8 void add_row_area(DCELL * top, DCELL * bottom, double sz, struct Cell_head *w,
9 double *low, double *high)
10 {
11 double guess1, guess2, mag, tedge1[3], tedge2[3], crossp[3];
12 int col;
13
14 for (col = 0; col < w->cols - 1; col++) {
15
16 /*
17 For each cell**, we triangulate the four corners in
18 two different ways, 1) UppperLeft to LowerRight diagonal
19 and 2) LowerLeft to UpperRight diagonal. Then we add the
20 smaller of the two areas to "low" and the greater of
21 the two areas to "high".
22
23 ** here, the "cell" is actually the quadrangle formed by
24 the center point of four cells, since these are the
25 known elevation points.
26 */
27
28 /* If NAN go to next or we get NAN for everything */
29 if (Rast_is_d_null_value(&(bottom[col + 1])) ||
30 Rast_is_d_null_value(&(top[col])) ||
31 Rast_is_d_null_value(&(top[col + 1])) ||
32 Rast_is_d_null_value(&(bottom[col]))
33 )
34 continue;
35
36 /* guess1 --- ul to lr diag */
37 {
38 tedge1[X] = w->ew_res;
39 tedge1[Y] = -w->ns_res;
40 tedge1[Z] = sz * (bottom[col + 1] - top[col]);
41
42 /* upper */
43 tedge2[X] = 0.0;
44 tedge2[Y] = w->ns_res;
45 tedge2[Z] = sz * (top[col + 1] - bottom[col + 1]);
46
47 v3cross(tedge1, tedge2, crossp);
48 v3mag(crossp, &mag);
49 guess1 = .5 * mag;
50
51 /* lower */
52 tedge2[X] = -w->ew_res;
53 tedge2[Y] = 0.0;
54 tedge2[Z] = sz * (bottom[col] - bottom[col + 1]);
55
56 v3cross(tedge1, tedge2, crossp);
57 v3mag(crossp, &mag);
58 guess1 += .5 * mag;
59 }
60
61 /* guess2 --- ll to ur diag */
62 {
63 tedge1[X] = w->ew_res;
64 tedge1[Y] = w->ns_res;
65 tedge1[Z] = sz * (top[col + 1] - bottom[col]);
66
67 /* upper */
68 tedge2[X] = -w->ew_res;
69 tedge2[Y] = 0.0;
70 tedge2[Z] = sz * (top[col + 1] - top[col + 1]);
71
72 v3cross(tedge1, tedge2, crossp);
73 v3mag(crossp, &mag);
74 guess2 = .5 * mag;
75
76 /* lower */
77 tedge2[X] = 0.0;
78 tedge2[Y] = -w->ns_res;
79 tedge2[Z] = sz * (bottom[col + 1] - top[col + 1]);
80
81 v3cross(tedge1, tedge2, crossp);
82 v3mag(crossp, &mag);
83 guess2 += .5 * mag;
84 }
85 *low += (guess1 < guess2) ? guess1 : guess2;
86 *high += (guess1 < guess2) ? guess2 : guess1;
87
88 } /* ea col */
89
90 }
91
92 /* calculate the running area of null data cells */
add_null_area(DCELL * rast,struct Cell_head * region,double * area)93 void add_null_area(DCELL * rast, struct Cell_head *region, double *area)
94 {
95 int col;
96
97 for (col = 0; col < region->cols; col++) {
98 if (Rast_is_d_null_value(&(rast[col]))) {
99 *area += region->ew_res * region->ns_res;
100 }
101 }
102 }
103
104 /* return the cross product v3 = v1 cross v2 */
v3cross(double v1[3],double v2[3],double v3[3])105 void v3cross(double v1[3], double v2[3], double v3[3])
106 {
107 v3[X] = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
108 v3[Y] = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
109 v3[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
110 }
111
112 /* magnitude of vector */
v3mag(double v1[3],double * mag)113 void v3mag(double v1[3], double *mag)
114 {
115 *mag = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
116 }
117
conv_value(double value,int units)118 double conv_value(double value, int units)
119 {
120 if (units == U_UNDEFINED)
121 return value;
122
123 return value * G_meters_to_units_factor_sq(units);
124 }
125