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

1""" 

2Surface area computation utilities for groundwater transport analysis. 

3 

4This module provides functions for computing surface areas and average heights 

5of geometric shapes used in groundwater transport calculations. 

6""" 

7 

8import numpy as np 

9 

10 

11def compute_average_heights(x_edges, y_edges, y_lower, y_upper): 

12 """ 

13 Compute average heights of clipped trapezoids. 

14 

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

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

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

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

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

20 

21 Parameters 

22 ---------- 

23 x_edges : numpy.ndarray 

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

25 y_edges : numpy.ndarray 

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

27 y_lower : float 

28 Lower horizontal clipping bound 

29 y_upper : float 

30 Upper horizontal clipping bound 

31 

32 Returns 

33 ------- 

34 avg_heights : numpy.ndarray 

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

36 shape (n_y-1, n_x-1) 

37 """ 

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

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

40 widths = np.diff(x_edges) 

41 

42 # Handle complete outside cases 

43 all_corners = np.stack([y_tl, y_tr, y_bl, y_br], axis=-1) 

44 outside_mask = np.all(all_corners >= y_upper, axis=-1) | np.all(all_corners <= y_lower, axis=-1) 

45 

46 # Vectorized area calculation using shoelace formula for clipped quadrilaterals 

47 y_tl_c, y_tr_c = np.clip(y_tl, y_lower, y_upper), np.clip(y_tr, y_lower, y_upper) 

48 y_bl_c, y_br_c = np.clip(y_bl, y_lower, y_upper), np.clip(y_br, y_lower, y_upper) 

49 

50 # Use exact shoelace formula for quadrilateral with vertices: 

51 # (0, y_bl_c), (width, y_br_c), (width, y_tr_c), (0, y_tl_c) 

52 areas = 0.5 * widths * np.abs(y_bl_c + y_br_c - y_tl_c - y_tr_c) 

53 

54 # Enhanced correction for cases where edges cross clipping bounds 

55 # This handles the specific geometry of sloped edge intersections 

56 top_crosses = ((y_tl - y_upper) * (y_tr - y_upper)) < 0 

57 correction = np.where( 

58 top_crosses & (y_tl > y_upper) & (y_tr < y_upper), 

59 widths * (y_tl - y_upper) ** 2 / (2 * np.maximum(y_tl - y_tr, 1e-15)), 

60 0, 

61 ) 

62 areas += correction 

63 areas = np.where(outside_mask, 0, areas) 

64 

65 return areas / widths