Coverage for C:\src\imod-python\imod\evaluate\streamfunction.py: 76%

45 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-04-08 13:27 +0200

1import imod 

2 

3 

4def _streamfunction(sfrf, sfff, dx, dy): 

5 # Scale total flux to length of cell traversed 

6 # --> no, this gives really akward patterns 

7 scale_factor = ( 

8 1.0 # (sfrf.dx ** 2 + sfff.dy ** 2) ** 0.5 / (dx ** 2 + dy ** 2) ** 0.5 

9 ) 

10 sfrfi = sfrf * -scale_factor # FRF is negative eastward, so flip 

11 sfffi = sfff * scale_factor # FFF is positive northward 

12 

13 # calculate projection of flux onto cross section line 

14 # Fi_proj = (Fi dotproduct Si) / |Si| 

15 # keep sign 

16 Fi_dot_si = sfrfi * sfrf.dx + sfffi * sfff.dy 

17 Fi_proj = Fi_dot_si / Fi_dot_si.ds 

18 

19 # calculate stream function by cumsumming from bottom upwards 

20 # method taken from https://olsthoorn.readthedocs.io/en/latest/07_stream_lines.html 

21 layindex = Fi_proj.layer.values 

22 Fi_proj = ( 

23 Fi_proj.reindex(layer=layindex[::-1]) 

24 .cumsum(dim="layer") 

25 .reindex(layer=layindex) 

26 ) 

27 return Fi_proj 

28 

29 

30def streamfunction_line(frf, fff, start, end): 

31 """Obtain the streamfunction for a line cross section through 

32 a three-dimensional flow field. The streamfunction is obtained 

33 by first projecting the horizontal flow components onto the provided 

34 cross-section. The streamfunction can be contoured to visualize stream lines. 

35 Stream lines are an efficient way to visualize groundwater flow. 

36 

37 Note, however, that the streamfunction is only defined in 2D, non-diverging, 

38 steady-state flow without sources and sinks. These assumption are violated even 

39 in a 2D model, but even more so in the approach followed here. Flow perpendicular 

40 to the cross-section will not be visualized. It is up to the user to choose 

41 cross-sections as perpendicular to the main flow direction as possible. 

42 

43 The 2D streamfunction and stream line visualization is based on work of Em. Prof. Olsthoorn. 

44 

45 Parameters 

46 ---------- 

47 frf: `xarray.DataArray` 

48 Three- (or higher) dimensional dataarray of flow component along the rows (FLOW RIGHT FACE). 

49 fff: `xarray.DataArray` 

50 Three- (or higher) dimensional dataarray of flow component along the columns (FLOW FRONT FACE). 

51 start: (2, ) array_like 

52 A latitude-longitude pair designating the start point of the cross 

53 section. 

54 end: (2, ) array_like 

55 A latitude-longitude pair designating the end point of the cross 

56 section. 

57 

58 Returns 

59 ------- 

60 `xarray.DataArray` 

61 The streamfunction projected on the cross-section between start and end coordinate, 

62 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

63 the "ds" coordinate. 

64 """ 

65 # interpolate frf and fff to cell center 

66 frf = 0.5 * frf + 0.5 * frf.shift({"x": -1}) 

67 fff = 0.5 * fff + 0.5 * fff.shift({"y": -1}) 

68 # get cross section of frf and fff 

69 sfrf = imod.select.cross_section_line(frf, start, end) 

70 sfff = imod.select.cross_section_line(fff, start, end) 

71 return _streamfunction(sfrf, sfff, frf.dx, fff.dy) 

72 

73 

74def streamfunction_linestring(frf, fff, linestring): 

75 """Obtain the streamfunction for a linestring cross section through 

76 a three-dimensional flow field. The streamfunction is obtained 

77 by first projecting the horizontal flow components onto the provided 

78 cross-section. The streamfunction can be contoured to visualize stream lines. 

79 Stream lines are an efficient way to visualize groundwater flow. 

80 

81 Note, however, that the streamfunction is only defined in 2D, non-diverging, 

82 steady-state flow without sources and sinks. These assumption are violated even 

83 in a 2D model, but even more so in the approach followed here. Flow perpendicular 

84 to the cross-section will not be visualized. It is up to the user to choose 

85 cross-sections as perpendicular to the main flow direction as possible. 

86 

87 The 2D streamfunction and stream line visualization is based on work of Em. Prof. Olsthoorn. 

88 

89 Parameters 

90 ---------- 

91 frf: `xarray.DataArray` 

92 Three- (or higher) dimensional dataarray of flow component along the rows (FLOW RIGHT FACE). 

93 fff: `xarray.DataArray` 

94 Three- (or higher) dimensional dataarray of flow component along the columns (FLOW FRONT FACE). 

95 

96 linestring : shapely.geometry.LineString 

97 Shapely geometry designating the linestring along which to sample the 

98 cross section. 

99 

100 Returns 

101 ------- 

102 `xarray.DataArray` 

103 The streamfunction projected on the cross-section defined by provided linestring, 

104 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

105 the "ds" coordinate. 

106 """ 

107 

108 # interpolate frf and fff to cell center 

109 frf = 0.5 * frf + 0.5 * frf.shift({"x": -1}) 

110 fff = 0.5 * fff + 0.5 * fff.shift({"y": -1}) 

111 # get cross section of frf and fff 

112 sfrf = imod.select.cross_section_linestring(frf, linestring) 

113 sfff = imod.select.cross_section_linestring(fff, linestring) 

114 return _streamfunction(sfrf, sfff, frf.dx, fff.dy) 

115 

116 

117def _quiver(sfrf, sfff, sflf): # , dx, dy, dz): 

118 # Scale total flux to length of cell traversed 

119 # --> no, this gives really akward patterns -> also for quivers? here it might seem like a good idea 

120 scale_factor = ( 

121 1.0 # (sfrf.dx ** 2 + sfff.dy ** 2) ** 0.5 / (dx ** 2 + dy ** 2) ** 0.5 

122 ) 

123 sfrfi = sfrf * -scale_factor # FRF is negative eastward, so flip 

124 sfffi = sfff * scale_factor # FFF is positive northward 

125 sflfi = sflf * scale_factor # FLF is positive upward 

126 

127 # calculate projection of flux onto cross section plane 

128 # cross section plane is given by? u: <dx, dy, 0>, v: <0, 0, 1> 

129 # Fi_proj_u = (Fi dotproduct Si) / |Si| * (Si / |Si|) 

130 # unitvector not necessary: return u and v as scalars 

131 # keep sign 

132 Fi_dot_u = sfrfi * sfrf.dx + sfffi * sfff.dy # + sflfi * 0. 

133 Fi_proj_u = Fi_dot_u / Fi_dot_u.ds 

134 Fi_proj_v = sflfi # no projection necessary for v 

135 

136 return Fi_proj_u, Fi_proj_v 

137 

138 

139def quiver_line(frf, fff, flf, start, end): 

140 """Obtain the u and v components for quiver plots for a line cross section 

141 through a three-dimensional flux field. The u and v components are obtained 

142 by first projecting the threedimensional flux components onto the provided 

143 cross-section. 

144 

145 Parameters 

146 ---------- 

147 frf: `xarray.DataArray` 

148 Three- (or higher) dimensional dataarray of flow component along the rows (FLOW RIGHT FACE). 

149 fff: `xarray.DataArray` 

150 Three- (or higher) dimensional dataarray of flow component along the columns (FLOW FRONT FACE). 

151 flf: `xarray.DataArray` 

152 Three- (or higher) dimensional dataarray of flow component along the layers (FLOW LOWER FACE). 

153 start: (2, ) array_like 

154 A latitude-longitude pair designating the start point of the cross 

155 section. 

156 end: (2, ) array_like 

157 A latitude-longitude pair designating the end point of the cross 

158 section. 

159 

160 Returns 

161 ------- 

162 u: `xarray.DataArray` 

163 The u component (x-component) of the flow projection on the cross-section between start and end coordinate, 

164 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

165 the "ds" coordinate. 

166 v: `xarray.DataArray` 

167 The v component (y-component) of the flow projection on the cross-section between start and end coordinate, 

168 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

169 the "ds" coordinate. 

170 

171 Notes 

172 ----- 

173 Use imod.evaluate.flow_velocity() first to obtain groundwater velocities 

174 as input for this function. Velocity in x direction is, however, inverted and must 

175 be re-inverted before using as input here. 

176 """ 

177 # interpolate frf and fff to cell center 

178 frf = 0.5 * frf + 0.5 * frf.shift({"x": -1}) 

179 fff = 0.5 * fff + 0.5 * fff.shift({"y": -1}) 

180 # get cross section of frf, fff and flf 

181 sfrf = imod.select.cross_section_line(frf, start, end) 

182 sfff = imod.select.cross_section_line(fff, start, end) 

183 sflf = imod.select.cross_section_line(flf, start, end) 

184 return _quiver(sfrf, sfff, sflf) 

185 

186 

187def quiver_linestring(frf, fff, flf, linestring): 

188 """Obtain the u and v components for quiver plots for a linestring cross section 

189 through a three-dimensional flow field. The u and v components are obtained 

190 by first projecting the threedimensional flow components onto the provided 

191 cross-section. 

192 

193 Parameters 

194 ---------- 

195 frf: `xarray.DataArray` 

196 Three- (or higher) dimensional dataarray of flow component along the rows (FLOW RIGHT FACE). 

197 fff: `xarray.DataArray` 

198 Three- (or higher) dimensional dataarray of flow component along the columns (FLOW FRONT FACE). 

199 flf: `xarray.DataArray` 

200 Three- (or higher) dimensional dataarray of flow component along the layers (FLOW LOWER FACE). 

201 linestring : shapely.geometry.LineString 

202 Shapely geometry designating the linestring along which to sample the 

203 cross section. 

204 

205 Returns 

206 ------- 

207 u: `xarray.DataArray` 

208 The u component (x-component) of the flow projection on the cross-section between start and end coordinate, 

209 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

210 the "ds" coordinate. 

211 v: `xarray.DataArray` 

212 The v component (y-component) of the flow projection on the cross-section between start and end coordinate, 

213 with new dimension "s" along the cross-section. The cellsizes along "s" are given in 

214 the "ds" coordinate. 

215 

216 Notes 

217 ----- 

218 Use imod.evaluate.flow_velocity() first to obtain groundwater velocities 

219 as input for this function. Velocity in x direction is, however, inverted and must 

220 be re-inverted before using as input here. 

221 """ 

222 # interpolate frf and fff to cell center 

223 frf = 0.5 * frf + 0.5 * frf.shift({"x": -1}) 

224 fff = 0.5 * fff + 0.5 * fff.shift({"y": -1}) 

225 # get cross section of frf, fff and flf 

226 sfrf = imod.select.cross_section_linestring(frf, linestring) 

227 sfff = imod.select.cross_section_linestring(fff, linestring) 

228 sflf = imod.select.cross_section_linestring(flf, linestring) 

229 return _quiver(sfrf, sfff, sflf)