Coverage for C:\src\imod-python\imod\flow\cap.py: 100%

50 statements  

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

1import os 

2import pathlib 

3 

4import jinja2 

5 

6from imod.flow.pkgbase import Package 

7 

8 

9class MetaSwap(Package): 

10 """ 

11 The MetaSWAP package (CAP), provides the input to be converted to a 

12 MetaSWAP model, which is an external model code used to simulate the 

13 unsaturated zone. 

14 

15 Note that only two-dimensional DataArrays with dimensions ``("y", "x")`` 

16 should be supplied to this package. In the current implementation 

17 time-related files are provided as external files ("extra files"). Similar 

18 to the iMODFLOW implementation of the projectfile. For now these need to be 

19 provided as a path. 

20 

21 MetaSWAP is developed by Alterra, Wageningen as part of the SIMGRO model 

22 code. The SIMGRO framework is intended for regions with an undulating 

23 topography and unconsolidated sediments in the (shallow) subsoil. Both 

24 shallow and deep groundwater levels can be modelled by MetaSWAP. This 

25 model is based on a simplification of ‘straight Richards’, meaning that no 

26 special processes like hysteresis, preferential flow and bypass flow are 

27 modelled. Snow is not modelled, and neither the influence of frost on the 

28 soil water conductivity. A perched watertable can be present in the SVAT 

29 column model, but interflow is not modelled. Processes that are typical for 

30 steep slopes are not included. The code contains several parameterized 

31 water management schemes, including irrigation and water level management. 

32 

33 References: 

34 

35 * Van Walsum, P. E. V., 2017a. SIMGRO V7.3.3.2, Input and Output reference 

36 manual. Tech. Rep. Alterra-Report 913.3, Alterra, Wageningen. 98 pp. 

37 

38 * Van Walsum, P. E. V., 2017b. SIMGRO V7.3.3.2, Users Guide. Tech. Rep. 

39 Alterra-Report 913.2, Alterra, Wageningen. 111 pp. 

40 

41 * Van Walsum, P. E. V. and P. Groenendijk, 2008. "Quasi Steady-State 

42 Simulation on of the Unsaturated Zone in Groundwater Modeling of Lowland 

43 Regions." Vadose Zone Journal 7: 769-778. 

44 

45 * Van Walsum, P. E. V., A. A. Veldhuizen and P. Groenendijk, 2016. SIMGRO 

46 V7.2.27, Theory and model implementation. Tech. Rep. Alterra-Report 913.1, 

47 Alterra, Wageningen. 93 pp 491. 

48 

49 Parameters 

50 ---------- 

51 boundary : int or xr.DataArray of ints 

52 2D boundary, used to specify active MetaSWAP elements, similar to 

53 ibound in the Boundary package 

54 

55 landuse : int or xr.DataArray of ints 

56 Landuse codes, referred to in the lookup table file luse_mswp.inp 

57 

58 rootzone_thickness : float or xr.DataArray of floats 

59 Rootzone thickness in cm (min. value is 10 centimeter). 

60 

61 soil_physical_unit : int or xr.DataArray of ints 

62 Soil Physical Unit, referred to in the lookup table file fact_mswp.inp. 

63 

64 meteo_station_number : float or xr.DataArray of ints 

65 Meteo station number, referred to by mete_mswp.inp. 

66 

67 surface_elevation : float or xr.DataArray of floats 

68 Surface Elevation (m+MSL) 

69 

70 sprinkling_type : int or xr.DataArray of ints 

71 Sprinkling type ("Artificial Recharge Type" in iMOD manual): 

72 

73 * 0 = no occurrence 

74 * 1 = from groundwater 

75 * 2 = from surface water 

76 

77 sprinkling_layer : int or xr.DataArray of ints 

78 Number of modellayer from which water is extracted ("Artificial 

79 Recharge Location" in iMOD manual) 

80 

81 sprinkling_capacity : float or xr.DataArray of floats 

82 Sprinkling capacity (mm/d) sets the maximum amount extracted for 

83 sprinkling ("Artificial Recharge Capacity" in iMOD manual) 

84 

85 wetted_area : float or xr.DataArray of floats 

86 Total area (m2) occupied by surface water elements. Values will be 

87 truncated by maximum cellsize. 

88 

89 urban_area : float or xr.DataArray of floats 

90 Total area (m2) occupied by urban area. Values will be truncated by 

91 maximum cellsize. 

92 

93 ponding_depth_urban : float or xr.DataArray of floats 

94 Ponding Depth Urban Area (m), specifying the acceptable depth of the 

95 ponding of water on the surface in the urban area before surface runoff 

96 occurs. 

97 

98 ponding_depth_rural : float or xr.DataArray of floats 

99 Ponding Depth Rural Area (m), specifying the acceptable depth of the 

100 ponding of water on the surface in the rural area before surface runoff 

101 occurs. 

102 

103 runoff_resistance_urban : float or xr.DataArray of floats 

104 Runoff Resistance Urban Area (day), specifying the resistance surface 

105 flow encounters in the urban area. The minimum value is equal to the 

106 model time period. 

107 

108 runoff_resistance_rural : float or xr.DataArray of floats 

109 Runoff Resistance Rural Area (day), specifying the resistance surface 

110 flow encounters in the rural area. The minimum value is equal to the 

111 model time period. 

112 

113 runon_resistance_urban : float or xr.DataArray of floats 

114 Runon Resistance Urban Area (day), specifying the resistance surface 

115 flow encounters to a model cell from an adjacent cell in the urban 

116 area. The minimum value is equal to the model time period. 

117 

118 runon_resistance_rural : float or xr.DataArray of floats 

119 Runon Resistance Rural Area (day), specifying the resistance surface 

120 flow encounters to a model cell from an adjacent cell in the rural 

121 area. The minimum value is equal to the model time period. 

122 

123 infiltration_capacity_urban : float or xr.DataArray of floats 

124 the infiltration capacity (m/d) of the soil surface in the urban area. 

125 The range is 0-1000 m/d. The NoDataValue -9999 indicates unlimited 

126 infiltration is possible. 

127 

128 infiltration_capacity_rural : float or xr.DataArray of floats 

129 the infiltration capacity (m/d) of the soil surface in the urban area. 

130 The range is 0-1000 m/d. The NoDataValue -9999 indicates unlimited 

131 infiltration is possible. 

132 

133 perched_water_table : float or xr.DataArray of floats 

134 Depth of the perched water table level (m) 

135 

136 soil_moisture_factor : float 

137 Soil Moisture Factor to adjust the soil moisture coefficient. This 

138 factor may be used during calibration. Default value is 1.0. 

139 

140 conductivity_factor : float 

141 Conductivity Factor to adjust the vertical conductivity. This factor 

142 may be used during calibration. Default value is 1.0. 

143 

144 lookup_and_forcing_files : list of pathlib.Path or str 

145 List of paths to files required by MetaSWAP. This a list of 

146 lookup tables and meteorological information that is required by 

147 MetaSwap. Note that MetaSwap looks for files with a specific name, so 

148 calling "luse_svat.inp" something else will result in errors. To view 

149 the files required, you can call: ``print(MetaSwap()._required_extra)`` 

150 

151 """ 

152 

153 _pkg_id = "cap" 

154 _variable_order = [ 

155 "boundary", 

156 "landuse", 

157 "rootzone_thickness", 

158 "soil_physical_unit", 

159 "meteo_station_number", 

160 "surface_elevation", 

161 "sprinkling_type", 

162 "sprinkling_layer", 

163 "sprinkling_capacity", 

164 "wetted_area", 

165 "urban_area", 

166 "ponding_depth_urban", 

167 "ponding_depth_rural", 

168 "runoff_resistance_urban", 

169 "runoff_resistance_rural", 

170 "runon_resistance_urban", 

171 "runon_resistance_rural", 

172 "infiltration_capacity_urban", 

173 "infiltration_capacity_rural", 

174 "perched_water_table", 

175 "soil_moisture_factor", 

176 "conductivity_factor", 

177 ] 

178 

179 _template_projectfile = jinja2.Template( 

180 "0001, ({{pkg_id}}), 1, {{name}}, {{variable_order}}\n" 

181 '{{"{:03d}".format(variable_order|length)}}, {{"{:03d}".format(n_entry)}}\n' 

182 "{%- for variable in variable_order%}\n" # Preserve variable order 

183 "{%- for layer, value in package_data[variable].items()%}\n" 

184 "{%- if value is string %}\n" 

185 # If string then assume path: 

186 # 1 indicates the layer is activated 

187 # 2 indicates the second element of the final two elements should be read 

188 # 1.000 is the multiplication factor 

189 # 0.000 is the addition factor 

190 # -9999 indicates there is no data, following iMOD usual practice 

191 '1, 2, {{"{:03d}".format(layer)}}, 1.000, 0.000, -9999., {{value}}\n' 

192 "{%- else %}\n" 

193 # Else assume a constant value is provided 

194 '1, 1, {{"{:03d}".format(layer)}}, 1.000, 0.000, {{value}}, ""\n' 

195 "{%- endif %}\n" 

196 "{%- endfor %}\n" 

197 "{%- endfor %}\n" 

198 # Section for EXTRA FILES comes below 

199 '{{"{:03d}".format(lookup_and_forcing_files|length)}},extra files\n' 

200 "{%- for file in lookup_and_forcing_files %}\n" 

201 "{{file}}\n" 

202 "{%- endfor %}\n" 

203 ) 

204 

205 # TODO: Check which of these actually are required. 

206 _required_extra = [ 

207 "fact_svat.inp", 

208 "luse_svat.inp", 

209 "mete_grid.inp", 

210 "para_sim.inp", 

211 "tiop_sim.inp", 

212 "init_svat.inp", 

213 "comp_post.inp", 

214 "sel_key_svat_per.inp", 

215 ] 

216 

217 def __init__( 

218 self, 

219 boundary, 

220 landuse, 

221 rootzone_thickness, 

222 soil_physical_unit, 

223 meteo_station_number, 

224 surface_elevation, 

225 sprinkling_type, 

226 sprinkling_layer, 

227 sprinkling_capacity, 

228 wetted_area, 

229 urban_area, 

230 ponding_depth_urban, 

231 ponding_depth_rural, 

232 runoff_resistance_urban, 

233 runoff_resistance_rural, 

234 runon_resistance_urban, 

235 runon_resistance_rural, 

236 infiltration_capacity_urban, 

237 infiltration_capacity_rural, 

238 perched_water_table, 

239 lookup_and_forcing_files, 

240 soil_moisture_factor=1.0, 

241 conductivity_factor=1.0, 

242 ): 

243 super().__init__() 

244 self.dataset["boundary"] = boundary 

245 self.dataset["landuse"] = landuse 

246 self.dataset["rootzone_thickness"] = rootzone_thickness 

247 self.dataset["soil_physical_unit"] = soil_physical_unit 

248 self.dataset["meteo_station_number"] = meteo_station_number 

249 self.dataset["surface_elevation"] = surface_elevation 

250 self.dataset["sprinkling_type"] = sprinkling_type 

251 self.dataset["sprinkling_layer"] = sprinkling_layer 

252 self.dataset["sprinkling_capacity"] = sprinkling_capacity 

253 self.dataset["wetted_area"] = wetted_area 

254 self.dataset["urban_area"] = urban_area 

255 self.dataset["ponding_depth_urban"] = ponding_depth_urban 

256 self.dataset["ponding_depth_rural"] = ponding_depth_rural 

257 self.dataset["runoff_resistance_urban"] = runoff_resistance_urban 

258 self.dataset["runoff_resistance_rural"] = runoff_resistance_rural 

259 self.dataset["runon_resistance_urban"] = runon_resistance_urban 

260 self.dataset["runon_resistance_rural"] = runon_resistance_rural 

261 self.dataset["infiltration_capacity_urban"] = infiltration_capacity_urban 

262 self.dataset["infiltration_capacity_rural"] = infiltration_capacity_rural 

263 self.dataset["perched_water_table"] = perched_water_table 

264 self.dataset["soil_moisture_factor"] = soil_moisture_factor 

265 self.dataset["conductivity_factor"] = conductivity_factor 

266 self.lookup_and_forcing_files = lookup_and_forcing_files 

267 

268 def _force_absolute_path(self, f): 

269 """Force absolute path, because projectfile cannot handle relative paths""" 

270 return str(pathlib.Path(f).resolve()) 

271 

272 def _render_projectfile(self, **kwargs): 

273 """ 

274 Returns 

275 ------- 

276 rendered : str 

277 The rendered projfectfile part, 

278 if part of PkgGroup: for a single boundary condition system. 

279 """ 

280 lookup_and_forcing_files = [ 

281 self._force_absolute_path(file) for file in self.lookup_and_forcing_files 

282 ] 

283 kwargs["lookup_and_forcing_files"] = lookup_and_forcing_files 

284 return self._template_projectfile.render(**kwargs) 

285 

286 def check_lookup_and_forcing_files(self): 

287 """Check for presence of required MetaSWAP input files.""" 

288 filenames = set([os.path.basename(f) for f in self.lookup_and_forcing_files]) 

289 missing_files = set(self._required_extra) - filenames 

290 if len(missing_files) > 0: 

291 raise ValueError(f"Missing MetaSWAP input files {missing_files}") 

292 

293 def _pkgcheck(self, active_cells=None): 

294 # Dataset.dims does not return a tuple, like DataArray does. 

295 # http://xarray.pydata.org/en/stable/generated/xarray.Dataset.dims.html 

296 # Frozen(SortedKeysDict).keys() does not preserve ordering in keys 

297 # Therefore convert to set and check for symmetric difference, 

298 # to see if any dimensions are missing or unexpected dimensions are included. 

299 dims = set(self.dataset.dims.keys()) 

300 

301 if len(dims.symmetric_difference(("y", "x"))) > 0: 

302 raise ValueError(f'Dataset dims not ("y", "x"), instead got {dims}') 

303 self.check_lookup_and_forcing_files()