Coverage for python / aubellhop / readers.py: 89%

1"""Reader functions for aubellhop. 

2 

3This file provides a collection of reading methods for both input 

4and output Bellhop files. 

5 

6Reading the environment file has substantial logic when relies on 

7the state of the environment information it is passed. Therefore, this is 

8managed by class `EnvironmentReader()`, used internally by the `Environment` 

9class via its class method `Environment.from_file()`. 

10 

11Subsequent files (ssp, bty, ati, brc, trc) are read automatically, but public 

12interfaces are provided for these for the purposes of mix and matching data 

13if desired. 

14 

15Finally, readers are also provides for output files for rays, arrivals, and 

16binary shd files. Again, these are used automatically by `aubellhop` after 

17executing a computation, but public methods are provides to help with 

18testing and verification. 

19 

20""" 

21 

22from __future__ import annotations 

23 

24import os 

25from struct import unpack as _unpack 

26from pathlib import Path 

27from typing import Any, TextIO, cast 

28from numpy.typing import NDArray 

29 

30import numpy as np 

31import pandas as pd 

32from aubellhop.constants import BHStrings, FlagMaps, FileExt, MiscDefaults 

33from aubellhop.environment import Environment 

34 

35 

36class EnvironmentReader: 

37 """Read and parse Bellhop environment files. 

38 

39 Although this class is only used for one task, 

40 the use of a class provides the clearest code interface compared 

41 to nested functions, which either implicitly set 

42 dict parameters, or have many repeated and superfluous 

43 arguments as dicts are passed in and returned at 

44 each stage. 

45 """ 

46 

47 def __init__(self, env: Environment, fname: str): 

48 """Initialize reader with filename. 

49 

50 Args: 

51 fname: Path to .env file (with or without extension) 

52 """ 

53 self.env = env 

54 self.fname, self.fname_base = _prepare_filename(fname, FileExt.env, "Environment") 

55 

56 def read(self) -> Environment: 

57 """Do the reading...""" 

58 with open(self.fname, 'r') as f: 

59 self._read_header(f) 

60 self._read_top_boundary(f) 

61 next_line = self._read_sound_speed_profile(f) 

62 self._read_bottom_boundary(f, next_line) 

63 self._read_sources_receivers_task(f) 

64 self._read_beams_limits(f) 

65 self._read_gaussian_params(f) 

66 return self.env 

67 

68 def _read_header(self, f: TextIO) -> None: 

69 """Read environment file header""" 

70 self.env['name'] = self._unquote_string(_read_next_valid_line(f)) 

71 self.env['frequency'] = _parse_vector(_read_next_valid_line(f)) 

72 self.env["_num_media"] = _parse_line_int(_read_next_valid_line(f)) 

73 

74 def _read_top_boundary(self, f: TextIO) -> None: 

75 """Read environment file top boundary options (multiple lines)""" 

76 

77 # Line 4: Top boundary options 

78 topopt_line = _read_next_valid_line(f) 

79 topopt = self._unquote_string(topopt_line) + " " 

80 self.env["soundspeed_interp"] = self._opt_lookup("Interpolation", topopt[0], FlagMaps.soundspeed_interp) 

81 self.env["surface_boundary_condition"] = self._opt_lookup("Top boundary condition", topopt[1], FlagMaps.surface_boundary_condition) 

82 self.env["attenuation_units"] = self._opt_lookup("Attenuation units", topopt[2], FlagMaps.attenuation_units) 

83 self.env["volume_attenuation"] = self._opt_lookup("Volume attenuation", topopt[3], FlagMaps.volume_attenuation) 

84 self.env["_altimetry"] = self._opt_lookup("Altimetry", topopt[4], FlagMaps._altimetry) 

85 self.env["_single_beam"] = self._opt_lookup("Single beam", topopt[5], FlagMaps._single_beam) 

86 if self.env["_altimetry"] == BHStrings.from_file: 

87 self.env["surface_depth"], self.env["surface_interp"] = read_ati(self.fname_base) 

88 

89 if self.env["volume_attenuation"] == BHStrings.francois_garrison: 

90 fg_spec_line = _read_next_valid_line(f) 

91 fg_parts = _parse_line(fg_spec_line) 

92 self.env["_fg_salinity"] = _float(fg_parts[0]) 

93 self.env["_fg_temperature"] = _float(fg_parts[1]) 

94 self.env["_fg_pH"] = _float(fg_parts[2]) 

95 self.env["_fg_depth"] = _float(fg_parts[3]) 

96 

97 # Line 4a: Boundary condition params 

98 if self.env["surface_boundary_condition"] == BHStrings.acousto_elastic: 

99 surface_props_line = _read_next_valid_line(f) 

100 surface_props = _parse_line(surface_props_line, none_pad=6) 

101 self.env['_surface_min'] = _float(surface_props[0]) 

102 self.env['surface_soundspeed'] = _float(surface_props[1]) 

103 self.env['_surface_soundspeed_shear'] = _float(surface_props[2]) 

104 self.env['surface_density'] = _float(surface_props[3], scale=1000) # convert from g/cm³ to kg/m³ 

105 self.env['surface_attenuation'] = _float(surface_props[4]) 

106 self.env['_surface_attenuation_shear'] = _float(surface_props[5]) 

107 

108 # Line 4b: Biological layer properties 

109 if self.env["volume_attenuation"] == BHStrings.biological: 

110 self.env['biological_layer_parameters'] = self._read_biological_layers(f) 

111 

112 def _read_biological_layers(self, f: TextIO) -> pd.DataFrame: 

113 """Read biological layer parameters for attenuation due to fish.""" 

114 next_line = _read_next_valid_line(f) 

115 npoints = int(next_line) 

116 z1 = [] 

117 z2 = [] 

118 f0 = [] 

119 QQ = [] 

120 a0 = [] 

121 for i in range(npoints): 

122 line = _read_next_valid_line(f) 

123 parts = _parse_line(line) 

124 if len(parts) == 5: 

125 z1.append(_float(parts[0])) 

126 z2.append(_float(parts[1])) 

127 f0.append(_float(parts[2])) 

128 QQ.append(_float(parts[3])) 

129 a0.append(_float(parts[4])) 

130 if len(z1) != npoints: 

131 raise ValueError(f"Expected {npoints} points, but found {len(z1)}") 

132 return pd.DataFrame({"z1": z1, "z2": z2, "f0": f0, "Q": QQ, "a0": a0}) 

133 

134 def _read_sound_speed_profile(self, f: TextIO) -> str: 

135 """Read environment file sound speed profile""" 

136 

137 # SSP depth specification 

138 ssp_spec_line = _read_next_valid_line(f) 

139 ssp_spec_line = ssp_spec_line.replace(",", " ") 

140 ssp_parts = _parse_line(ssp_spec_line, none_pad=3) 

141 self.env['_mesh_npts'] = _int(ssp_parts[0]) 

142 self.env['_depth_sigma'] = _float(ssp_parts[1]) 

143 self.env['_depth_max'] = _float(ssp_parts[2]) 

144 self.env['bottom_depth'] = self.env['_depth_max'] # TODO: interrogate this, looks unsafe 

145 

146 # Read SSP points and from file if applicable 

147 ssp_lines, next_line = self._read_until_quote(f) 

148 self.env['soundspeed'] = self._read_ssp_points(ssp_lines) 

149 if self.env["soundspeed_interp"] == BHStrings.quadrilateral: 

150 self.env['soundspeed'] = read_ssp(self.fname_base, self.env['soundspeed'].index) 

151 elif self.env["soundspeed_interp"] == BHStrings.hexahedral: 

152 self.env['soundspeed'] = read_ssp_3d(self.fname_base) 

153 return next_line 

154 

155 def _read_until_quote(self, f: TextIO) -> tuple[list[str],str]: 

156 """Read lines until one starts with ' character.""" 

157 lines: list[str] = [] 

158 while True: 

159 line = f.readline() 

160 if not line: 

161 raise EOFError("File ended during env file reading of SSP points.") 

162 line = line.strip() 

163 if not line: # completely empty line 

164 continue 

165 if line.startswith("'"): # Check if this is a bottom boundary line (starts with quote) 

166 return lines, line 

167 lines.append(line) 

168 

169 def _read_ssp_points(self, lines: list[str]) -> pd.DataFrame: 

170 """Read sound speed profile points until we find the bottom boundary line 

171 

172 Default values are according to 'EnvironmentalFile.htm'.""" 

173 

174 ssp_depth: list[float] = [] 

175 ssp_speed: list[float] = [] 

176 ssp_shear: list[float] = [] 

177 ssp_density: list[float] = [] 

178 ssp_atten: list[float] = [] 

179 ssp_att_shear: list[float] = [] 

180 ssp = dict(depth=0.0, speed=MiscDefaults.sound_speed, speed_shear=0.0, density=MiscDefaults.density, atten=0.0, att_shear=0.0) 

181 

182 MAX_SSP_COLS: int = 6 

183 sound_speed_param_count: int = 0 

184 for line in lines: 

185 raw_parts = _parse_line(line) 

186 parts = (raw_parts + [None] * MAX_SSP_COLS)[:MAX_SSP_COLS] 

187 if parts[0] is None: # empty line after stripping comments 

188 continue 

189 sound_speed_param_count = max(len(raw_parts), sound_speed_param_count) 

190 for k, v in zip(ssp, parts): 

191 if v is not None: 

192 ssp[k] = cast(float, _float(v)) # "fill in the blanks" if empty entries on any lines 

193 ssp_depth.append(ssp["depth"]) 

194 ssp_speed.append(ssp["speed"]) 

195 ssp_shear.append(ssp["speed_shear"]) 

196 ssp_density.append(ssp["density"] * 1000.0) # units scaling 

197 ssp_atten.append(ssp["atten"]) 

198 ssp_att_shear.append(ssp["att_shear"]) 

199 

200 if len(ssp_speed) == 0: 

201 raise ValueError("No SSP points were found in the env file.") 

202 elif len(ssp_speed) == 1: 

203 raise ValueError("Only one SSP point found but at least two required (top and bottom)") 

204 

205 df = pd.DataFrame({ 

206 "speed": ssp_speed, 

207 "shear_speed": ssp_shear, 

208 "density": ssp_density, 

209 "attenuation": ssp_atten, 

210 "shear_attenuation": ssp_att_shear 

211 }, index=pd.Index(ssp_depth)) 

212 df = df.iloc[:, :(sound_speed_param_count-1)] # Keep only the max number of columns read 

213 df.index.name = "depth" 

214 return df 

215 

216 def _read_bottom_boundary(self, f: TextIO, bottom_line: str) -> None: 

217 """Read environment file bottom boundary condition""" 

218 bottom_parts = _parse_line(bottom_line,none_pad=3) 

219 botopt = self._unquote_string(cast(str,bottom_parts[0])) + " " # cast() => I promise this is a str :) 

220 self.env["bottom_boundary_condition"] = self._opt_lookup("Bottom boundary condition", botopt[0], FlagMaps.bottom_boundary_condition) 

221 self.env["_bathymetry"] = self._opt_lookup("Bathymetry", botopt[1], FlagMaps._bathymetry) 

222 self.env['bottom_roughness'] = _float(bottom_parts[1]) 

223 self.env['bottom_beta'] = _float(bottom_parts[2]) 

224 self.env['bottom_transition_freq'] = _float(bottom_parts[3]) 

225 if self.env["_bathymetry"] == BHStrings.from_file: 

226 self.env["bottom_depth"], self.env["bottom_interp"] = read_bty(self.fname_base) 

227 

228 # Bottom properties (depth, sound_speed, density, absorption) 

229 if self.env["bottom_boundary_condition"] == BHStrings.acousto_elastic: 

230 bottom_props_line = _read_next_valid_line(f) 

231 bottom_props = _parse_line(bottom_props_line,none_pad=6) 

232 self.env['_bottom_depth'] = _float(bottom_props[0]) 

233 self.env['bottom_soundspeed'] = _float(bottom_props[1]) 

234 self.env['_bottom_soundspeed_shear'] = _float(bottom_props[2]) 

235 self.env['bottom_density'] = _float(bottom_props[3], 1000) # convert from g/cm³ to kg/m³ 

236 self.env['bottom_attenuation'] = _float(bottom_props[4]) 

237 self.env['_bottom_attenuation_shear'] = _float(bottom_props[5]) 

238 elif self.env["bottom_boundary_condition"] == BHStrings.grain: 

239 bottom_props_line = _read_next_valid_line(f) 

240 bottom_props = _parse_line(bottom_props_line, none_pad=6) 

241 self.env['_bottom_depth'] = _float(bottom_props[0]) 

242 self.env['bottom_grain_size'] = _float(bottom_props[1]) 

243 

244 def _read_sources_receivers_task(self, f: TextIO) -> None: 

245 """Read environment file sources, receivers, and task. 

246 

247 Bellhop and Bellhop3D have different numbers of variables specified before 

248 the task line. Luckily we can detect that reliably by looking for a line which 

249 starts with `'` to mark the definition of the tasks.""" 

250 

251 sr_lines = [] 

252 while True: 

253 next_line = _read_next_valid_line(f) 

254 if next_line.startswith("'"): 

255 break 

256 sr_lines.append(next_line) 

257 

258 self._parse_src_rcv(sr_lines) 

259 self._parse_task(next_line) 

260 

261 def _parse_src_rcv(self, sr_lines: list[str]) -> None: 

262 """Parse the N lines read defining sources and receivers and assign the corresponding variables.""" 

263 nlines = len(sr_lines) 

264 if nlines == 3: # sometimes see a shorthand version like: ['1 25.0 /', '10 100.0 1000.0 /', '1, 80.0'] 

265 self.env['_dimension'] = 2 

266 val0 = np.asarray(_parse_vector(sr_lines[0])) 

267 val1 = np.asarray(_parse_vector(sr_lines[1])) 

268 val2 = np.asarray(_parse_vector(sr_lines[2])) 

269 self.env['source_ndepth'] = int(val0[0]) 

270 self.env['receiver_nrange'] = int(val1[0]) 

271 self.env['receiver_ndepth'] = int(val2[0]) 

272 self.env['source_depth'] = val0[1:] 

273 self.env['receiver_depth'] = val1[1:] 

274 self.env['receiver_range'] = val2[1:] * 1000.0 # convert km to m 

275 elif nlines == 6: 

276 self.env['_dimension'] = 2 

277 self.env['source_ndepth'] = _parse_line_int(sr_lines[0]) 

278 self.env['receiver_ndepth'] = _parse_line_int(sr_lines[2]) 

279 self.env['receiver_nrange'] = _parse_line_int(sr_lines[4]) 

280 self.env['source_depth'] = _parse_vector(sr_lines[1]) 

281 self.env['receiver_depth'] = _parse_vector(sr_lines[3]) 

282 self.env['receiver_range'] = _parse_vector(sr_lines[5]) * 1000.0 # convert km to m 

283 elif nlines == 12: 

284 self.env['_dimension'] = 3 

285 self.env['source_nrange'] = _parse_line_int(sr_lines[0]) 

286 self.env['source_ncrossrange'] = _parse_line_int(sr_lines[2]) 

287 self.env['source_ndepth'] = _parse_line_int(sr_lines[4]) 

288 self.env['receiver_ndepth'] = _parse_line_int(sr_lines[6]) 

289 self.env['receiver_nrange'] = _parse_line_int(sr_lines[8]) 

290 self.env['receiver_nbearing'] = _parse_line_int(sr_lines[10]) 

291 self.env['source_range'] = _parse_vector(sr_lines[1]) * 1000.0 # convert km to m 

292 self.env['source_cross_range'] = _parse_vector(sr_lines[3]) * 1000.0 # convert km to m 

293 self.env['source_depth'] = _parse_vector(sr_lines[5]) 

294 self.env['receiver_depth'] = _parse_vector(sr_lines[7]) 

295 self.env['receiver_range'] = _parse_vector(sr_lines[9]) * 1000.0 # convert km to m 

296 self.env['receiver_bearing'] = _parse_vector(sr_lines[11]) 

297 else: 

298 print("SCANNED SRC/RCV LINES:") 

299 print(sr_lines) 

300 raise RuntimeError( 

301 "The python parsing of Bellhop's so-called 'list-directed IO' is not robust." 

302 f"Expected to read 6 or 12 lines (2D or 3D cases); found: {nlines}") 

303 

304 def _parse_task(self, task_line: str) -> None: 

305 """Parse the 'task' line.""" 

306 task_code = self._unquote_string(task_line) + " " 

307 self.env['task'] = FlagMaps.task.get(task_code[0]) 

308 self.env['beam_type'] = FlagMaps.beam_type.get(task_code[1]) 

309 self.env['_sbp_file'] = FlagMaps._sbp_file.get(task_code[2]) 

310 self.env['source_type'] = FlagMaps.source_type.get(task_code[3]) 

311 self.env['grid_type'] = FlagMaps.grid_type.get(task_code[4]) 

312 if self.env['_dimension'] == 2: 

313 self.env['dimension'] = BHStrings.two_d 

314 else: 

315 self.env['dimension'] = FlagMaps.dimension.get(task_code[5]) 

316 

317 if self.env["_sbp_file"] == BHStrings.from_file: 

318 self.env["source_directionality"] = read_sbp(self.fname_base) 

319 

320 def _read_beams_limits(self, f: TextIO) -> None: 

321 """Read environment file beams and limits""" 

322 

323 # Number of beams 

324 beam_num_line = _read_next_valid_line(f) 

325 beam_num_parts = _parse_line(beam_num_line, none_pad=1) 

326 self.env['beam_num'] = int(beam_num_parts[0] or 0) 

327 if self.env["_single_beam"]: # defensive in case there is a spurious value in here 

328 self.env['single_beam_index'] = _int(beam_num_parts[1]) 

329 

330 # Beam angles (beam_angle_min, beam_angle_max) 

331 angles_line = _read_next_valid_line(f) 

332 angle_parts = _parse_line(angles_line, none_pad=2) 

333 self.env['beam_angle_min'] = _float(angle_parts[0]) 

334 self.env['beam_angle_max'] = _float(angle_parts[1]) 

335 

336 if self.env['_dimension'] == 3: 

337 # Beam bearing fan 

338 beam_num_line = _read_next_valid_line(f) 

339 beam_num_parts = _parse_line(beam_num_line, none_pad=1) 

340 self.env['beam_angle_num'] = int(beam_num_parts[0] or 0) 

341 angles_line = _read_next_valid_line(f) 

342 angle_parts = _parse_line(angles_line, none_pad=2) 

343 self.env['beam_bearing_min'] = _float(angle_parts[0]) 

344 self.env['beam_bearing_max'] = _float(angle_parts[1]) 

345 

346 # Ray tracing limits (step, max_depth, max_range) - last line 

347 limits_line = _read_next_valid_line(f) 

348 limits_parts = _parse_line(limits_line) 

349 self.env['step_size'] = _float(limits_parts[0]) 

350 if self.env['_dimension'] == 2: 

351 self.env['simulation_depth'] = _float(limits_parts[1]) 

352 self.env['simulation_range'] = _float(limits_parts[2], 1000.0) # convert km to m 

353 else: 

354 self.env['simulation_range'] = _float(limits_parts[1], 1000.0) # convert km to m 

355 self.env['simulation_cross_range'] = _float(limits_parts[2], 1000.0) # convert km to m 

356 self.env['simulation_depth'] = _float(limits_parts[3]) 

357 

358 def _read_gaussian_params(self, f: TextIO) -> None: 

359 """Read parameters for Cerveny Gaussian Beams, if applicable""" 

360 if self.env['beam_type'] not in (BHStrings.gaussian_simple, BHStrings.ray): 

361 return None 

362 line = _read_next_valid_line(f) 

363 parts = _parse_line(line, none_pad=3) 

364 assert isinstance(parts[0],str) 

365 self.env['beam_width_type'] = self._unquote_string(parts[0]) 

366 self.env['beam_epsilon_multipler'] = _float(parts[1]) 

367 self.env['beam_range_loop'] = _float(parts[2],1000) 

368 

369 line = _read_next_valid_line(f) 

370 parts = _parse_line(line, none_pad=3) 

371 self.env['beam_images_num'] = _int(parts[0]) 

372 self.env['beam_window'] = _int(parts[1]) 

373 self.env['beam_component'] = self._unquote_string(parts[2]) if parts[2] is not None else " " 

374 

375 def _opt_lookup(self, name: str, opt: str, _map: dict[str, BHStrings]) -> str | None: 

376 opt_str = _map.get(opt) 

377 if opt_str is None: 

378 raise ValueError(f"{name} option {opt!r} not available") 

379 return opt_str 

380 

381 def _unquote_string(self, line: str) -> str: 

382 """Extract string from within single quotes, possibly with commas too.""" 

383 return line.strip().strip(",'") 

384 

385######################################## 

386 

387def read_ssp(fname: str, 

388 depths: list[float] | NDArray[np.float64] | pd.DataFrame | None = None 

389 ) -> NDArray[np.float64] | pd.DataFrame: 

390 """Read a 2D sound speed profile (.ssp) file used by BELLHOP. 

391 

392 This function reads BELLHOP's .ssp files which contain range-dependent 

393 sound speed profiles. The file format is: 

394 - Line 1: Number of range profiles (NPROFILES) 

395 - Line 2: Range coordinates in km (space-separated) 

396 - Line 3+: Sound speed values, one line per depth point across all ranges 

397 

398 Parameters 

399 ---------- 

400 fname : str 

401 Path to .ssp file (with or without .ssp extension) 

402 

403 Returns 

404 ------- 

405 numpy.ndarray or pandas.DataFrame 

406 For single-profile files: numpy array with [depth, soundspeed] pairs; 

407 for multi-profile files: pandas DataFrame with range-dependent sound speed data 

408 

409 Notes 

410 ----- 

411 **Return format:** 

412 

413 - **Single-profile files (1 range)**: Returns a 2D numpy array with [depth, soundspeed] pairs, 

414 compatible with the `Environment()` soundspeed parameter. 

415 

416 - **Multi-profile files (>1 ranges)**: Returns a pandas DataFrame where: 

417 

418 - **Columns**: Range coordinates (in meters, converted from km in file) 

419 - **Index**: Depth indices (0, 1, 2, ... for each depth level in the file) 

420 - **Values**: Sound speeds (m/s) 

421 

422 This DataFrame can be directly assigned to the `Environment()` soundspeed parameter 

423 for range-dependent acoustic modeling. 

424 

425 **Note on depths**: For multi-profile files, depth indices are used (0, 1, 2, ...) 

426 since the actual depth coordinates come from the associated BELLHOP .env file. 

427 Users can modify the DataFrame index if actual depth values are known. 

428 

429 Examples 

430 -------- 

431 >>> import aubellhop as bh 

432 >>> # Single-profile file 

433 >>> ssp1 = bh.read_ssp("single_profile.ssp") # Returns numpy array 

434 >>> env = bh.Environment() 

435 >>> env["soundspeed"] = ssp1 

436 >>> 

437 >>> # Multi-profile file 

438 >>> ssp2 = bh.read_ssp("tests/MunkB_geo_rot/MunkB_geo_rot.ssp") # Returns DataFrame 

439 >>> env = bh.Environment() 

440 >>> env["soundspeed"] = ssp2 # Range-dependent sound speed 

441 

442 **File format example:** 

443 

444 :: 

445 

446 30 

447 -50 -5 -1 -.8 -.75 -.6 -.4 -.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 10.0 

448 1500 1500 1548.52 1530.29 1526.69 1517.78 1509.49 1504.30 1501.38 1500.14 1500.12 1501.02 1502.57 1504.62 1507.02 1509.69 1512.55 1515.56 1518.67 1521.85 1525.10 1528.38 1531.70 1535.04 1538.39 1541.76 1545.14 1548.52 1551.91 1551.91 

449 1500 1500 1548.52 1530.29 1526.69 1517.78 1509.49 1504.30 1501.38 1500.14 1500.12 1501.02 1502.57 1504.62 1507.02 1509.69 1512.55 1515.56 1518.67 1521.85 1525.10 1528.38 1531.70 1535.04 1538.39 1541.76 1545.14 1548.52 1551.91 1551.91 

450 """ 

451 

452 fname, _ = _prepare_filename(fname, FileExt.ssp, "SSP") 

453 with open(fname, 'r') as f: 

454 nranges = int(_read_next_valid_line(f)) 

455 range_line = _read_next_valid_line(f) 

456 ranges = np.array([_float(x) for x in _parse_line(range_line)]) 

457 ranges_m = ranges * 1000 # Convert ranges from km to meters (as expected by Environment()) 

458 

459 if len(ranges) != nranges: 

460 raise ValueError(f"Expected {nranges} ranges, but found {len(ranges)}") 

461 

462 # Read sound speed data - read all remaining lines as a matrix 

463 ssp_data = [] 

464 line_num = 0 

465 for line in f: 

466 line_num += 1 

467 line = line.replace(","," ").strip() 

468 if line: # Skip empty lines 

469 values = [_float(x) for x in line.split()] 

470 if len(values) != nranges: 

471 raise ValueError(f"SSP line {line_num} has {len(values)} range values, expected {nranges}") 

472 ssp_data.append(values) 

473 

474 ssp_array = np.array(ssp_data) 

475 ndepths = ssp_array.shape[0] 

476 

477 # Create depth indices (actual depths would normally come from associated .env file) 

478 if depths is None: 

479 depths = np.arange(ndepths, dtype=float) 

480 

481 if ndepths == 0 or len(depths) != ndepths: 

482 raise ValueError("Wrong number of depths found in sound speed data file" 

483 f" (expected {ndepths}, found {ssp_array.shape[0]})") 

484 

485 df = pd.DataFrame(ssp_array, index=pd.Index(depths), columns=ranges_m) 

486 df.index.name = "depth" 

487 return df 

488 

489def read_ssp_3d(fname: str) -> dict[str, np.ndarray]: 

490 """Read a 3D sound speed profile (.ssp) file used by BELLHOP3D and return a 3D Numpy array. 

491 

492 This function reads BELLHOP's .ssp files which contain range-dependent 

493 sound speed profiles. The file format is: 

494 - Line 1: Nx: Number of coordinates in the range (x) direction 

495 - Line 2: Range coordinates in km (space-separated) 

496 - Line 3: Ny: Number of coordinates in the crossrange (y) direction 

497 - Line 4: Crossrange coordinates in km (space-separated) 

498 - Line 5: Nz: Number of coordinates in the depth (z) direction 

499 - Line 6: Depth coordinates in km (space-separated) 

500 - Line 7+: Sound speed values, one line per crossrange and depth across all ranges 

501 (i.e., each line has Nx number of values and there are Ny*Nz lines) 

502 """ 

503 

504 fname, _ = _prepare_filename(fname, FileExt.ssp, "SSP") 

505 with open(fname, 'r') as f: 

506 

507 nranges = int(_read_next_valid_line(f)) 

508 range_line = _read_next_valid_line(f) 

509 ncrossranges = int(_read_next_valid_line(f)) 

510 crossrange_line = _read_next_valid_line(f) 

511 ndepths = int(_read_next_valid_line(f)) 

512 depth_line = _read_next_valid_line(f) 

513 

514 ranges = np.array([_float(x) for x in _parse_line(range_line)]) 

515 crossranges = np.array([_float(x) for x in _parse_line(crossrange_line)]) 

516 depths = np.array([_float(x) for x in _parse_line(depth_line)]) 

517 

518 # Convert ranges from km to meters (as expected by Environment()) 

519 ranges_m = ranges * 1000 

520 crossranges_m = crossranges * 1000 

521 depths_m = depths * 1000 

522 

523 if len(ranges) != nranges: 

524 raise ValueError(f"Expected {nranges} ranges, but found {len(ranges)}") 

525 if len(crossranges) != ncrossranges: 

526 raise ValueError(f"Expected {ncrossranges} ranges, but found {len(crossranges)}") 

527 if len(depths) != ndepths: 

528 raise ValueError(f"Expected {ndepths} ranges, but found {len(depths)}") 

529 

530 # Read sound speed data - read all remaining lines as a matrix 

531 ssp_data = [] 

532 line_num = 0 

533 for line in f: 

534 line_num += 1 

535 line = line.replace(","," ").strip() 

536 if line: # Skip empty lines 

537 values = [_float(x) for x in line.split()] 

538 if len(values) != nranges: 

539 raise ValueError(f"SSP line {line_num} has {len(values)} range values, expected {nranges}") 

540 ssp_data.append(values) 

541 

542 ssp_array = np.array(ssp_data) 

543 narray = ssp_array.shape[0] 

544 if narray != ncrossranges * ndepths: 

545 raise ValueError("Wrong number of entries found in sound speed data file" 

546 f" (expected {ncrossranges} * {ndepths} = {ncrossranges * ndepths} lines, found {narray})") 

547 

548 ssp_3d = ssp_array.reshape(ndepths, ncrossranges, nranges, order="C") 

549 ssp_xyz = np.transpose(ssp_3d, (2, 1, 0)) 

550 

551 return {"ssp": ssp_xyz, "ranges": ranges_m, "crossranges": crossranges_m, "depths": depths_m} 

552 

553def read_ati(fname: str) -> tuple[NDArray[np.float64], str]: 

554 """Read an altimetry file used by Bellhop.""" 

555 fname, _ = _prepare_filename(fname, FileExt.ati, "ATI") 

556 return _read_ati_bty(fname) 

557 

558def read_bty(fname: str) -> tuple[NDArray[np.float64], str]: 

559 """Read a bathymetry file used by Bellhop.""" 

560 fname, _ = _prepare_filename(fname, FileExt.bty, "BTY") 

561 return _read_ati_bty(fname) 

562 

563def _read_ati_bty(fname: str) -> tuple[NDArray[np.float64], str]: 

564 """Read an altimetry (.ati) or bathymetry (.bty) file used by BELLHOP. 

565 

566 This function reads BELLHOP's .bty files which define the bottom depth 

567 profile. The file format is: 

568 - Line 1: Interpolation type ('L' for linear, 'C' for curvilinear) 

569 - Line 2: Number of points 

570 - Line 3+: Range (km) and depth (m) pairs 

571 

572 Parameters 

573 ---------- 

574 fname : str 

575 Path to .ati/.bty file (with extension) 

576 

577 Returns 

578 ------- 

579 numpy.ndarray 

580 Numpy array with [range, depth] pairs compatible with Environment() 

581 

582 Notes 

583 ----- 

584 The returned array can be assigned to env["bottom_depth"] for range-dependent bathymetry. 

585 

586 **Examples:** 

587 

588 >>> import aubellhop as bh 

589 >>> bty,bty_interp = bh.read_bty("tests/MunkB_geo_rot/MunkB_geo_rot.bty") 

590 >>> env = bh.Environment() 

591 >>> env["bottom_depth"] = bty 

592 >>> env["bottom_interp"] = bty_interp 

593 >>> arrivals = bh.calculate_arrivals(env) 

594 

595 **File format example:** 

596 

597 :: 

598 

599 'L' 

600 5 

601 0 3000 

602 10 3000 

603 20 500 

604 30 3000 

605 100 3000 

606 """ 

607 

608 with open(fname, 'r') as f: 

609 # Read interpolation type ('L' or 'C') 

610 interp_type = _read_next_valid_line(f).strip("'\"") 

611 nvalues = 2 

612 if len(interp_type) > 1: 

613 nvalues = 7 if interp_type[1] == "L" else 2 

614 interp_type = interp_type[0] 

615 npoints = int(_read_next_valid_line(f)) 

616 ranges = [] 

617 depths = [] 

618 wave_speed = [] 

619 wave_attenuation = [] 

620 density = [] 

621 shear_speed = [] 

622 shear_attenuation = [] 

623 for i in range(npoints): 

624 try: 

625 line = _read_next_valid_line(f) 

626 except EOFError: 

627 break 

628 parts = _parse_line(line) 

629 if nvalues == 2 and len(parts) >= 2: 

630 ranges.append(_float(parts[0])) # Range in km 

631 depths.append(_float(parts[1])) # Depth in m 

632 elif nvalues == 7 and len(parts) == 7: 

633 ranges.append(_float(parts[0])) # Range in km 

634 depths.append(_float(parts[1])) # Depth in m 

635 wave_speed.append(_float(parts[2])) # 

636 wave_attenuation.append(_float(parts[3])) # 

637 density.append(_float(parts[4])) # 

638 shear_speed.append(_float(parts[5])) # 

639 shear_attenuation.append(_float(parts[6])) # 

640 

641 if len(ranges) != npoints: 

642 raise ValueError(f"Expected {npoints} altimetry/bathymetry points, but found {len(ranges)}") 

643 

644 # Convert ranges from km to m for consistency with bellhop env structure 

645 ranges_m = np.array(ranges) * 1000 

646 depths_array = np.array(depths) 

647 if nvalues == 2: 

648 val_array = [ranges_m, depths_array] 

649 elif nvalues == 7: 

650 wave_speed_array = np.array(wave_speed) 

651 wave_attenuation_array = np.array(wave_attenuation) 

652 density_array = np.array(density) 

653 shear_speed_array = np.array(shear_speed) 

654 shear_attenuation_array = np.array(shear_attenuation) 

655 val_array = [ 

656 ranges_m, 

657 depths_array, 

658 wave_speed_array, 

659 wave_attenuation_array, 

660 density_array, 

661 shear_speed_array, 

662 shear_attenuation_array, 

663 ] 

664 return np.column_stack(val_array), FlagMaps.bottom_interp[interp_type] 

665 

666def read_ati_3d(fname: str) -> dict[str, np.ndarray]: 

667 """Read an altimetry file used by Bellhop.""" 

668 fname, _ = _prepare_filename(fname, FileExt.ati, "ATI") 

669 return _read_ati_bty3d(fname) 

670 

671def read_bty_3d(fname: str) -> dict[str, np.ndarray]: 

672 """Read a bathymetry file used by Bellhop.""" 

673 fname, _ = _prepare_filename(fname, FileExt.bty, "BTY") 

674 return _read_ati_bty3d(fname) 

675 

676def _read_ati_bty3d(fname: str) -> dict[str, NDArray[np.float64]]: 

677 """Read an altimetry (.ati) or bathymetry (.bty) file used by BELLHOP. 

678 

679 This function reads BELLHOP's .bty files which define the bottom depth 

680 profile. The file format is: 

681 - Line 1: Interpolation type ('R' always for Bellhop3D) 

682 - Line 2: Nx: Number of range points 

683 - Line 3: Range coordinates (x) 

684 - Line 4: Ny: Number of crossrange points 

685 - Line 5: Crossrange points (y) 

686 - Line 6+: 2D array of depths. Ny lines, one per crossrange, with Nx depths per line (space separated) 

687 

688 Parameters 

689 ---------- 

690 fname : str 

691 Path to .ati/.bty file (with extension) 

692 

693 Returns 

694 ------- 

695 dict 

696 Dictionary of ranges, crossranges, and 2D depth data 

697 

698 Notes 

699 -----