Coverage for python/bellhop/readers.py: 91%
256 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-09-23 13:34 +0000
2import numpy as _np
3import pandas as _pd
4from bellhop.constants import _Strings, _Maps
5import bellhop.environment
7def read_env2d(fname):
8 """Read a 2D underwater environment from a BELLHOP .env file.
10 This function parses a BELLHOP .env file and returns a Python data structure
11 that is compatible with create_env2d(). This enables round-trip testing and
12 compatibility between file-based and programmatic environment definitions.
14 :param fname: path to .env file (with or without .env extension)
16 :returns: environment dictionary compatible with create_env2d()
18 The environment dictionary used in this code contains a large
19 number of parameters, documented here to keep the code later more concise:
21 ENV parameters
22 ---------------
24 :name: environment title/name
25 :type: '2D' (fixed for 2D environments)
26 :frequency: acoustic frequency in Hz
27 :soundspeed: sound speed profile (scalar for constant, array for depth-dependent)
28 :soundspeed_interp: interpolation method ('linear', 'spline', 'quadrilateral')
29 :bottom_soundspeed: bottom sediment sound speed in m/s
30 :bottom_soundspeed_shear: bottom sediment sound speed in m/s
31 :bottom_density: bottom sediment density in kg/m³
32 :bottom_absorption: bottom sediment absorption in dB/wavelength
33 :bottom_absorption_shear: bottom sediment absorption in dB/wavelength
34 :bottom_roughness: bottom roughness RMS in meters
35 :surface: surface altimetry profile (None if flat surface)
36 :surface_interp: surface interpolation method ('linear', 'curvilinear')
37 :top_boundary_condition: ('vacuum', 'acousto-elastic', 'rigid', 'from-file')
38 :volume_attenuation: ('none', 'thorp', 'francois-garrison', 'biological')
39 :attenuation_units: ('nepers per meter', 'frequency dependent', 'dB per meter', 'frequency scaled dB per meter', 'dB per wavelength', 'quality factor', 'loss parameter')
40 :tx_depth: transmitter depth(s) in meters
41 :tx_directionality: transmitter beam pattern (None if omnidirectional)
42 :rx_depth: receiver depth(s) in meters
43 :rx_range: receiver range(s) in meters
44 :depth: maximum water depth in meters
45 :depth_interp: bathymetry interpolation method ('linear', 'curvilinear')
46 :min_angle: minimum beam angle in degrees
47 :max_angle: maximum beam angle in degrees
48 :nbeams: number of beams (0 for automatic)
49 :step_size: (maximum) step size to trace rays in meters (0 for automatic)
50 :box_depth: box extent to trace rays in meters (auto-calculated based on max depth data if not specified)
51 :box_range: box extent to trace rays in meters (auto-calculated based on max receiver range if not specified)
52 :tx_type: point (default) or line
53 :beam_type: todo
54 :grid: rectilinear or irregular
56 **Supported ENV file formats:**
58 - Standard BELLHOP format with various boundary conditions
59 - Constant or depth-dependent sound speed profiles
60 - Compressed vector notation (e.g., "0.0 5000.0 /" for linearly spaced values)
61 - Comments (lines with ! are handled correctly)
62 - Different top/bottom boundary options (halfspace, file-based, etc.)
64 **Unit conversions performed:**
66 - Receiver ranges: km → m
67 - Bottom density: g/cm³ → kg/m³
68 - All other units preserved as in ENV file
70 **Examples:**
72 >>> import bellhop as bh
73 >>> env = bh.read_env2d('examples/Munk/MunkB_ray.env')
74 >>> print(env['name'])
75 'Munk profile'
76 >>> print(env['frequency'])
77 50.0
79 >>> # Use with existing functions
80 >>> checked_env = bh.check_env2d(env)
81 >>> rays = bh.compute_rays(env)
83 >>> # Round-trip compatibility
84 >>> env_orig = bh.create_env2d(name="test", frequency=100)
85 >>> # ... write to file via BELLHOP ...
86 >>> env_read = bh.read_env2d("test.env")
87 >>> assert env_read['frequency'] == env_orig['frequency']
89 **Limitations:**
91 - External files (.ssp, .bty, .ati, .sbp) are noted but not automatically loaded
92 - Some advanced BELLHOP features may not be fully supported
93 - Assumes standard 2D BELLHOP format (not BELLHOP3D)
94 """
95 import os
96 import re
98 # Add .env extension if not present
99 if not fname.endswith('.env'):
100 fname = fname + '.env'
102 if not os.path.exists(fname):
103 raise FileNotFoundError(f"Environment file not found: {fname}")
105 # Initialize environment with default values
106 env = bellhop.environment.new()
108 def _parse_quoted_string(line):
109 """Extract string from within quotes
111 Sometimes (why??) the leading quote was being stripped, so we also try to catch
112 this case with the regexp, stripping only a trailing '.
113 """
114 mtch = re.search(r"'([^']*)'", line)
115 mtch2 = re.search(r"([^']*)'$", line)
116 return mtch.group(1) if mtch else mtch2.group(1) if mtch2 else line.strip()
118 def _parse_line(line):
119 """Parse a line, removing comments and whitespace"""
120 # Remove comments (everything after !)
121 if '!' in line:
122 line = line[:line.index('!')].strip()
123 return line.strip()
125 def _parse_vector(f, dtype=float):
126 """Parse a vector that starts with count then values, ending with '/'"""
127 line = f.readline().strip()
128 if not line:
129 raise ValueError("Unexpected end of file while reading vector")
131 # First line is the count
132 linecount = int(_parse_line(line))
134 # Second line has the values
135 values_line = f.readline().strip()
136 values_line = _parse_line(values_line)
138 # Split by '/' and take only the first part (before the '/')
139 if '/' in values_line:
140 values_line = values_line.split('/')[0].strip()
142 parts = values_line.split()
143 values = [dtype(p) for p in parts]
145 # Note that we do not try to interpolate here, since Bellhop has its own routines
146 return _np.array(values), linecount
148 def _read_ssp_points(f):
149 """Read sound speed profile points until we find the bottom boundary line"""
150 ssp_points = []
152 while True:
153 line = f.readline().strip()
154 if not line:
155 break
157 # Check if this is a bottom boundary line (starts with quote)
158 if line.startswith("'"):
159 # This is the bottom boundary line, put it back
160 f.seek(f.tell() - len(line.encode()) - 1)
161 break
163 # Parse SSP point
164 line = _parse_line(line)
165 if line.endswith('/'):
166 line = line[:-1].strip()
168 parts = line.split()
169 if len(parts) >= 2:
170 try:
171 depth = float(parts[0])
172 speed = float(parts[1])
173 ssp_points.append([depth, speed])
174 except ValueError:
175 # This might be the end of SSP or a different format
176 # Put the line back and break
177 f.seek(f.tell() - len(line.encode()) - 1)
178 break
180 return _np.array(ssp_points) if ssp_points else None
182 with open(fname, 'r') as f:
183 # Line 1: Title
184 title_line = f.readline().strip()
185 env['name'] = _parse_quoted_string(title_line)
187 # Line 2: Frequency
188 freq_line = f.readline().strip()
189 env['frequency'] = float(_parse_line(freq_line))
191 # Line 3: NMedia (should be 1 for BELLHOP)
192 nmedia_line = f.readline().strip()
193 nmedia = int(_parse_line(nmedia_line))
194 if nmedia != 1:
195 raise ValueError(f"BELLHOP only supports 1 medium, found {nmedia}")
197 # Line 4: Top boundary options
198 topopt_line = f.readline().strip()
199 topopt = _parse_quoted_string(topopt_line)
201 # Parse SSP interpolation type from first character
202 def _invalid(opt):
203 raise ValueError(f"Interpolation option {opt!r} not available")
204 opt = topopt[0]
205 env["soundspeed_interp"] = _Maps.interp.get(opt) or _invalid(opt)
207 # Top boundary condition
208 def _invalid(opt):
209 raise ValueError(f"Top boundary condition option {opt!r} not available")
210 opt = topopt[1]
211 env["top_boundary_condition"] = _Maps.boundcond.get(opt) or _invalid(opt)
213 # Attenuation units
214 def _invalid(opt):
215 raise ValueError(f"Attenuation units option {opt!r} not available")
216 opt = topopt[2]
217 env["attenuation_units"] = _Maps.attunits.get(opt) or _invalid(opt)
219 # Volume attenuation
220 def _invalid(opt):
221 raise ValueError(f"Volume attenuation option {opt!r} not available")
222 env["volume_attenuation"] = 'none'
223 if len(topopt) > 3:
224 opt = topopt[3]
225 else:
226 opt = ""
227 env["volume_attenuation"] = _Maps.volatt.get(opt) or _invalid(opt)
229 if env["volume_attenuation"] == _Strings.francois_garrison:
230 fg_spec_line = f.readline().strip()
231 fg_parts = _parse_line(fg_spec_line).split()
232 env["fg_salinity"] = float(fg_parts[0])
233 env["fg_temperature"] = float(fg_parts[1])
234 env["fg_pH"] = float(fg_parts[2])
235 env["fg_depth"] = float(fg_parts[3])
238 if 'A' in topopt:
239 # Read halfspace parameters line
240 f.readline().strip()
241 # This line contains: depth, alphaR, betaR, rho, alphaI, betaI
242 # We skip this for now as it's not part of the standard env structure
244 # Check for surface altimetry (indicated by * in topopt)
245 if '*' in topopt:
246 # Surface altimetry file exists - would need to read .ati file
247 # For now, just note that surface is present
248 env['surface'] = _np.array([[0, 0], [1000, 0]]) # placeholder
250 # Line 5 or 6: SSP depth specification (format: npts sigma_z max_depth)
251 ssp_spec_line = f.readline().strip()
252 ssp_parts = _parse_line(ssp_spec_line).split()
253 env['depth_npts'] = int(ssp_parts[0])
254 if len(ssp_parts) > 1:
255 env['depth_sigmaz'] = float(ssp_parts[1])
256 if len(ssp_parts) > 2:
257 env['depth_max'] = float(ssp_parts[2])
258 env['depth'] = float(ssp_parts[2])
260 # Read SSP points
261 ssp_points = _read_ssp_points(f)
262 if ssp_points is not None and len(ssp_points) > 0:
263 if len(ssp_points) == 1:
264 # Single sound speed value
265 env['soundspeed'] = ssp_points[0, 1]
266 else:
267 # Multiple points - depth, sound speed pairs
268 env['soundspeed'] = ssp_points
269 env['ssp_env'] = ssp_points
271 # Bottom boundary options
272 print(" ")
273 print(fname)
275 line = f.readline()
276 bottom_line = line.strip()
277 bottom_parts = _parse_line(bottom_line).split()
278 botopt = _parse_quoted_string(bottom_parts[0])
279 def _invalid(opt):
280 raise ValueError(f"Bottom boundary condition option {opt!r} not available")
281 opt = botopt[0]
282 env["bottom_boundary_condition"] = _Maps.boundcond.get(opt) or _invalid(opt)
284 if len(botopt) > 1:
285 opt = botopt[1]
286 env["_bottom_bathymetry"] = _Maps.bottom.get(opt) or _invalid(opt)
287 if env["_bottom_bathymetry"] == _Strings.from_file:
288 print("TODO: automatically read bty file")
289 else:
290 pass # nothing needs to be done
292 if len(bottom_parts) >= 2:
293 env['bottom_roughness'] = float(bottom_parts[1])
294 if len(bottom_parts) >= 3:
295 env['bottom_beta'] = float(bottom_parts[2])
296 env['bottom_transition_freq'] = float(bottom_parts[3])
298 # Bottom properties (depth, sound_speed, density, absorption)
299 bottom_props_line = f.readline().strip()
300 bottom_props_line = _parse_line(bottom_props_line)
301 if bottom_props_line.endswith('/'):
302 bottom_props_line = bottom_props_line[:-1].strip()
304 bottom_props = bottom_props_line.split()
305 # fortran sources say: "z, alphaR, betaR, rhoR, alphaI, betaI"
306 # docs say:
307 # Syntax:
308 #
309 # ZB CPB CSB RHOB APB ASB
310 #
311 # Description:
312 #
313 # ZB: Depth (m).
314 # CPB: Bottom P-wave speed (m/s).
315 # CSB: Bottom S-wave speed (m/s).
316 # RHOB: Bottom density (g/cm3).
317 # APB: Bottom P-wave attenuation. (units as given by TOPOPT(3:3) )
318 # ASB: Bottom S-wave attenuation. ( " " " " " " )
319 if len(bottom_props) > 1:
320 env['bottom_soundspeed'] = float(bottom_props[1])
321 if len(bottom_props) > 2:
322 env['bottom_soundspeed_shear'] = float(bottom_props[2])
323 if len(bottom_props) > 3:
324 env['bottom_density'] = float(bottom_props[3]) * 1000 # convert from g/cm³ to kg/m³
325 if len(bottom_props) > 4:
326 env['bottom_absorption'] = float(bottom_props[4])
327 if len(bottom_props) > 5:
328 env['bottom_absorption_shear'] = float(bottom_props[5])
330 # Source depths
331 tx_depths, env['tx_ndepth'] = _parse_vector(f)
332 if len(tx_depths) == 1:
333 env['tx_depth'] = tx_depths[0]
334 else:
335 env['tx_depth'] = tx_depths
337 # Receiver depths
338 rx_depths, env['rx_ndepth'] = _parse_vector(f)
339 if len(rx_depths) == 1:
340 env['rx_depth'] = rx_depths[0]
341 else:
342 env['rx_depth'] = rx_depths
344 # Receiver ranges (in km, need to convert to m)
345 rx_ranges, env['rx_nrange'] = _parse_vector(f)
346 env['rx_range'] = rx_ranges * 1000 # convert km to m
348 # Task/run type (e.g., 'R', 'C', etc.)
349 task_line = f.readline().strip()
350 task_code = _parse_quoted_string(task_line)
351 env['task'] = task_code[0]
352 if len(task_code) > 1:
353 env['beam_type'] = _Maps.beam.get(task_code[1])
354 if len(task_code) > 3:
355 env['tx_type'] = _Maps.source.get(task_code[3])
356 if len(task_code) > 4:
357 env['grid'] = _Maps.grid.get(task_code[4])
359 # Check for source directionality (indicated by * in task code)
360 if '*' in task_code:
361 # Source directionality file exists - would need to read .sbp file
362 # For now, just note that directionality is present
363 env['tx_directionality'] = _np.array([[0, 0]]) # placeholder
365 # Number of beams
366 nbeams_line = f.readline().strip()
367 env['nbeams'] = int(_parse_line(nbeams_line))
369 # Beam angles (min_angle, max_angle)
370 angles_line = f.readline().strip()
371 angles_line = _parse_line(angles_line)
372 if angles_line.endswith('/'):
373 angles_line = angles_line[:-1].strip()
375 angle_parts = angles_line.split()
376 if len(angle_parts) >= 2:
377 env['min_angle'] = float(angle_parts[0])
378 env['max_angle'] = float(angle_parts[1])
380 # Ray tracing limits (step, max_depth, max_range) - last line
381 limits_line = f.readline().strip()
382 limits_parse = _parse_line(limits_line)
383 limits_parts = limits_parse.split()
384 env['step_size'] = float(limits_parts[0])
385 env['box_depth'] = float(limits_parts[1])
386 env['box_range'] = 1000*float(limits_parts[2])
388 return env
390def read_ssp(fname):
391 """Read a 2D sound speed profile (.ssp) file used by BELLHOP.
393 This function reads BELLHOP's .ssp files which contain range-dependent
394 sound speed profiles. The file format is:
395 - Line 1: Number of range profiles (NPROFILES)
396 - Line 2: Range coordinates in km (space-separated)
397 - Line 3+: Sound speed values, one line per depth point across all ranges
399 :param fname: path to .ssp file (with or without .ssp extension)
400 :returns: for single-profile files: numpy array with [depth, soundspeed] pairs;
401 for multi-profile files: pandas DataFrame with range-dependent sound speed data
403 **Return format:**
405 - **Single-profile files (1 range)**: Returns a 2D numpy array with [depth, soundspeed] pairs,
406 compatible with create_env2d() soundspeed parameter.
408 - **Multi-profile files (>1 ranges)**: Returns a pandas DataFrame where:
410 - **Columns**: Range coordinates (in meters, converted from km in file)
411 - **Index**: Depth indices (0, 1, 2, ... for each depth level in the file)
412 - **Values**: Sound speeds (m/s)
414 This DataFrame can be directly assigned to create_env2d() soundspeed parameter
415 for range-dependent acoustic modeling.
417 **Note on depths**: For multi-profile files, depth indices are used (0, 1, 2, ...)
418 since the actual depth coordinates come from the associated BELLHOP .env file.
419 Users can modify the DataFrame index if actual depth values are known.
421 **Examples:**
423 >>> import bellhop as bh
424 >>> # Single-profile file
425 >>> ssp1 = bh.read_ssp("single_profile.ssp") # Returns numpy array
426 >>> env = bh.create_env2d()
427 >>> env["soundspeed"] = ssp1
428 >>>
429 >>> # Multi-profile file
430 >>> ssp2 = bh.read_ssp("tests/MunkB_geo_rot/MunkB_geo_rot.ssp") # Returns DataFrame
431 >>> env = bh.create_env2d()
432 >>> env["soundspeed"] = ssp2 # Range-dependent sound speed
434 **File format example:**
436 ::
438 30
439 -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
440 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
441 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
442 """
443 import os
445 # Add .ssp extension if not present
446 if not fname.endswith('.ssp'):
447 fname = fname + '.ssp'
449 if not os.path.exists(fname):
450 raise FileNotFoundError(f"SSP file not found: {fname}")
452 with open(fname, 'r') as f:
453 # Read number of range profiles
454 nprofiles = int(f.readline().strip())
456 # Read range coordinates (in km)
457 range_line = f.readline().strip()
458 ranges = _np.array([float(x) for x in range_line.split()])
460 if len(ranges) != nprofiles:
461 raise ValueError(f"Expected {nprofiles} range profiles, but found {len(ranges)} ranges")
463 # Read sound speed data - read all remaining lines as a matrix
464 ssp_data = []
465 for line in f:
466 line = line.strip()
467 if line: # Skip empty lines
468 values = [float(x) for x in line.split()]
469 if len(values) == nprofiles:
470 ssp_data.append(values)
472 ssp_array = _np.array(ssp_data)
474 if ssp_array.size == 0:
475 raise ValueError("No sound speed data found in file")
477 if nprofiles == 1:
478 # Single profile - return as [depth, soundspeed] pairs for backward compatibility
479 # Create depth values - linearly spaced from 0 to number of depth points
480 ndepths = ssp_array.shape[0]
481 depths = _np.linspace(0, ndepths-1, ndepths, dtype=float)
482 return _np.column_stack([depths, ssp_array.flatten()])
483 else:
484 # Multiple ranges - return as pandas DataFrame for range-dependent modeling
485 # Convert ranges from km to meters (as expected by create_env2d)
486 ranges_m = ranges * 1000
488 # Create depth indices (actual depths would come from associated .env file)
489 ndepths = ssp_array.shape[0]
490 depths = _np.arange(ndepths, dtype=float)
492 # Create DataFrame with ranges as columns and depths as index
493 # ssp_array is [ndepths, nprofiles] which is the correct orientation
494 return _pd.DataFrame(ssp_array, index=depths, columns=ranges_m)
496def read_bty(fname):
497 """Read a bathymetry (.bty) file used by BELLHOP.
499 This function reads BELLHOP's .bty files which define the bottom depth
500 profile. The file format is:
501 - Line 1: Interpolation type ('L' for linear, 'C' for curvilinear)
502 - Line 2: Number of points
503 - Line 3+: Range (km) and depth (m) pairs
505 :param fname: path to .bty file (with or without .bty extension)
506 :returns: numpy array with [range, depth] pairs compatible with create_env2d()
508 The returned array can be assigned to env["depth"] for range-dependent bathymetry.
510 **Examples:**
512 >>> import bellhop as bh
513 >>> bty,bty_interp = bh.read_bty("tests/MunkB_geo_rot/MunkB_geo_rot.bty")
514 >>> env = bh.create_env2d()
515 >>> env["depth"] = bty
516 >>> env["depth_interp"] = bty_interp
517 >>> arrivals = bh.calculate_arrivals(env)
519 **File format example:**
521 ::
523 'L'
524 5
525 0 3000
526 10 3000
527 20 500
528 30 3000
529 100 3000
530 """
531 import os
533 # Add .bty extension if not present
534 if not fname.endswith('.bty'):
535 fname = fname + '.bty'
537 if not os.path.exists(fname):
538 raise FileNotFoundError(f"BTY file not found: {fname}")
540 with open(fname, 'r') as f:
541 # Read interpolation type (usually 'L' or 'C')
542 interp_type = f.readline().strip().strip("'\"")
544 # Read number of points
545 npoints = int(f.readline().strip())
547 # Read range,depth pairs
548 ranges = []
549 depths = []
551 for i in range(npoints):
552 line = f.readline().strip()
553 if line: # Skip empty lines
554 parts = line.split()
555 if len(parts) >= 2:
556 ranges.append(float(parts[0])) # Range in km
557 depths.append(float(parts[1])) # Depth in m
559 if len(ranges) != npoints:
560 raise ValueError(f"Expected {npoints} bathymetry points, but found {len(ranges)}")
562 # Convert ranges from km to m for consistency with bellhop env structure
563 ranges_m = _np.array(ranges) * 1000
564 depths_array = _np.array(depths)
566 # Return as [range, depth] pairs
567 return _np.column_stack([ranges_m, depths_array]), _Maps.bty_interp[interp_type]
569def read_refl_coeff(fname):
570 """Read a reflection coefficient (.brc) file used by BELLHOP.
572 This function reads BELLHOP's .brc files which define the reflection coefficient
573 data. The file format is:
574 - Line 1: Number of points
575 - Line 2+: THETA(j) RMAG(j) RPHASE(j)
577 Where:
578 - THETA(): Angle (degrees)
579 - RMAG(): Magnitude of reflection coefficient
580 - RPHASE(): Phase of reflection coefficient (degrees)
582 :param fname: path to .brc/.trc file (extension required)
583 :returns: numpy array with [theta, rmag, rphase] triplets compatible with create_env2d()
585 The returned array can be assigned to env["bottom_reflection_coefficient"] or env["top_reflection_coefficient"] .
587 **Example:**
589 >>> import bellhop as bh
590 >>> brc = bh.read_refl_coeff("tests/MunkB_geo_rot/MunkB_geo_rot.brc")
591 >>> env = bh.create_env2d()
592 >>> env["bottom_reflection_coefficient"] = brc
593 >>> arrivals = bh.calculate_arrivals(env)
595 **File format example:**
597 ::
599 3
600 0.0 1.00 180.0
601 45.0 0.95 175.0
602 90.0 0.90 170.0
603 """
604 import os
607 if not os.path.exists(fname):
608 raise FileNotFoundError(f"Reflection coefficient file not found: {fname}")
610 with open(fname, 'r') as f:
612 # Read number of points
613 npoints = int(f.readline().strip())
615 # Read range,depth pairs
616 theta = []
617 rmagn = []
618 rphas = []
620 for i in range(npoints):
621 line = f.readline().strip()
622 if line: # Skip empty lines
623 parts = line.split()
624 if len(parts) == 3:
625 theta.append(float(parts[0]))
626 rmagn.append(float(parts[1]))
627 rphas.append(float(parts[2]))
629 if len(theta) != npoints:
630 raise ValueError(f"Expected {npoints} bathymetry points, but found {len(theta)}")
632 # Return as [range, depth] pairs
633 return _np.column_stack([theta, rmagn, rphas])