Coverage for src / gwtransport / surfacearea.py: 100%

1""" 

2Surface Area Calculations for Streamline-Based Transport Analysis. 

3 

4This module provides geometric utilities for computing surface areas and average 

5heights between streamlines in heterogeneous aquifer systems. These calculations 

6support direct estimation of pore volume distributions from streamline analysis, 

7providing an alternative to gamma distribution approximations. 

8 

9Available functions: 

10 

11- :func:`compute_average_heights` - Compute average heights of clipped trapezoids formed by 

12 streamlines. Trapezoids have vertical sides defined by x_edges and top/bottom edges defined 

13 by y_edges (2D array). Clipping bounds (y_lower, y_upper) restrict the integration domain. 

14 Returns area/width ratios representing average heights for use in pore volume calculations 

15 from streamline geometry. 

16 

17This file is part of gwtransport which is released under AGPL-3.0 license. 

18See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details. 

19""" 

20 

21import numpy as np 

22import numpy.typing as npt 

23 

24 

25def _positive_part_integral( 

26 a: npt.NDArray[np.floating], b: npt.NDArray[np.floating], w: npt.NDArray[np.floating] 

27) -> npt.NDArray[np.floating]: 

28 """ 

29 Integrate max(f(x), 0) from x=0 to x=w where f is linear from a to b. 

30 

31 Parameters 

32 ---------- 

33 a : ndarray 

34 Function values at x=0. 

35 b : ndarray 

36 Function values at x=w. 

37 w : ndarray 

38 Integration width. 

39 

40 Returns 

41 ------- 

42 ndarray 

43 Integral values. 

44 """ 

45 both_pos = (a > 0) & (b > 0) 

46 only_a_pos = (a > 0) & ~both_pos 

47 only_b_pos = (b > 0) & ~both_pos 

48 

49 abs_diff = np.abs(a - b) 

50 safe_diff = np.where(abs_diff > 0, abs_diff, 1.0) 

51 

52 excess = np.where(only_a_pos, a, np.where(only_b_pos, b, 0.0)) 

53 

54 return np.where( 

55 both_pos, 

56 w * (a + b) / 2, 

57 w * excess**2 / (2 * safe_diff), 

58 ) 

59 

60 

61def _clipped_linear_integral( 

62 a: npt.NDArray[np.floating], 

63 b: npt.NDArray[np.floating], 

64 w: npt.NDArray[np.floating], 

65 lo: float, 

66 hi: float, 

67) -> npt.NDArray[np.floating]: 

68 """ 

69 Integrate clip(f(x), lo, hi) from x=0 to x=w where f is linear from a to b. 

70 

71 Uses the identity ``clip(f) = f - max(f - hi, 0) + max(lo - f, 0)`` to 

72 compute the exact integral analytically. 

73 

74 Parameters 

75 ---------- 

76 a : ndarray 

77 Function values at x=0. 

78 b : ndarray 

79 Function values at x=w. 

80 w : ndarray 

81 Integration width. 

82 lo : float 

83 Lower clipping bound. 

84 hi : float 

85 Upper clipping bound. 

86 

87 Returns 

88 ------- 

89 ndarray 

90 Integral values. 

91 """ 

92 raw = w * (a + b) / 2 

93 excess_above = _positive_part_integral(a - hi, b - hi, w) 

94 deficit_below = _positive_part_integral(lo - a, lo - b, w) 

95 return raw - excess_above + deficit_below 

96 

97 

98def compute_average_heights( 

99 *, x_edges: npt.ArrayLike, y_edges: npt.ArrayLike, y_lower: float, y_upper: float 

100) -> npt.NDArray[np.floating]: 

101 """ 

102 Compute average heights of clipped trapezoids. 

103 

104 Trapezoids have vertical left and right sides, with corners at: 

105 - top-left: (y_edges[i, j], x_edges[j]) 

106 - top-right: (y_edges[i, j+1], x_edges[j+1]) 

107 - bottom-left: (y_edges[i+1, j], x_edges[j]) 

108 - bottom-right: (y_edges[i+1, j+1], x_edges[j+1]) 

109 

110 The area is computed as the exact integral of 

111 ``clip(top(x), y_lower, y_upper) - clip(bottom(x), y_lower, y_upper)`` 

112 where ``top(x)`` and ``bottom(x)`` are linear interpolants of the corner 

113 values, and ``clip`` restricts values to ``[y_lower, y_upper]``. 

114 

115 Parameters 

116 ---------- 

117 x_edges : numpy.ndarray 

118 1D array of x coordinates, shape (n_x,) 

119 y_edges : numpy.ndarray 

120 2D array of y coordinates, shape (n_y, n_x) 

121 y_lower : float 

122 Lower horizontal clipping bound 

123 y_upper : float 

124 Upper horizontal clipping bound 

125 

126 Returns 

127 ------- 

128 avg_heights : numpy.ndarray 

129 2D array of average heights (area/width) for each clipped trapezoid, 

130 shape (n_y-1, n_x-1) 

131 

132 See Also 

133 -------- 

134 gwtransport.advection.infiltration_to_extraction : Use computed areas in transport calculations 

135 

136 Examples 

137 -------- 

138 >>> import numpy as np 

139 >>> from gwtransport.surfacearea import compute_average_heights 

140 >>> # Create simple grid 

141 >>> x_edges = np.array([0.0, 1.0, 2.0]) 

142 >>> y_edges = np.array([[0.0, 0.5, 1.0], [1.0, 1.5, 2.0], [2.0, 2.5, 3.0]]) 

143 >>> # Compute average heights with clipping bounds 

144 >>> avg_heights = compute_average_heights( 

145 ... x_edges=x_edges, y_edges=y_edges, y_lower=0.5, y_upper=2.5 

146 ... ) 

147 >>> print(f"Shape: {avg_heights.shape}") 

148 Shape: (2, 2) 

149 """ 

150 x_edges = np.asarray(x_edges, dtype=float) 

151 y_edges = np.asarray(y_edges, dtype=float) 

152 

153 y_tl, y_tr = y_edges[:-1, :-1], y_edges[:-1, 1:] 

154 y_bl, y_br = y_edges[1:, :-1], y_edges[1:, 1:] 

155 widths = np.diff(x_edges) 

156 

157 top_integral = _clipped_linear_integral(y_tl, y_tr, widths, y_lower, y_upper) 

158 bottom_integral = _clipped_linear_integral(y_bl, y_br, widths, y_lower, y_upper) 

159 

160 areas = np.maximum(top_integral - bottom_integral, 0.0) 

161 

162 with np.errstate(invalid="ignore"): 

163 return areas / widths