ADCIRC+SWAN耦合模型数据可视化实战：用Python处理fort.63/.64和swan输出文件并绘图

ADCIRCSWAN耦合模型数据可视化实战Python处理与专业图表生成风暴潮与波浪耦合模拟是海岸工程研究的重要工具但模拟结果的深度分析往往被繁琐的数据处理流程所阻碍。许多研究者花费大量时间成功运行ADCIRCSWAN模型后却在数据可视化阶段遇到瓶颈——如何从原始的fort.63、fort.64和swan输出文件中提取有价值的信息并转化为具有发表质量的科学图表这正是本文要解决的核心问题。1. 模型输出文件解析与Python读取技巧ADCIRCSWAN耦合模型会生成多种输出文件每种文件都承载着特定的海洋动力学信息。理解这些文件的内部结构是进行有效可视化的第一步。主要输出文件类型及内容文件格式对应变量数据维度典型应用场景fort.63风暴潮水位时间×空间增水分析、淹没范围fort.64流速矢量时间×空间流场动力学、物质输运swan_HS.63有效波高时间×空间波浪能量评估swan_DIR.63平均波向时间×空间波浪传播方向分析swan_TPS.63波峰周期时间×空间波浪频谱特性这些文件虽然扩展名不同但都遵循类似的ASCII格式结构。以fort.63为例其典型结构如下ADCIRC仿真结果 1 ! 文件类型标识 2023 07 15 00 00 00 ! 起始时间 3600.0 ! 时间步长(s) 2 ! 变量数量 1 ! 节点总数 ... TIME 0.0 ! 第一个时间步 1 23.456 ! 节点编号 水位值(m) 2 24.567 ... TIME 3600.0 ! 第二个时间步 1 24.123 2 25.234 ...在Python中我们可以使用numpy和pandas高效读取这些文件import numpy as np import pandas as pd def read_adcirc_output(filename): with open(filename, r) as f: header [next(f) for _ in range(4)] # 读取前4行头文件 nodes int(header[3].strip()) # 获取节点总数 data [] time_marks [] current_time None for line in f: if line.startswith(TIME): current_time float(line.split()[1].strip()) time_marks.append(current_time) continue parts line.strip().split() if len(parts) 2: # 水位或波高数据 node_id, value map(float, parts) data.append([current_time, node_id, value]) df pd.DataFrame(data, columns[time, node, value]) return df, header注意fort.64和fort.74文件包含矢量数据U/V分量读取时需要特别处理两列数值数据。2. 时空数据重构与网格映射技术原始输出数据是按时间步和节点编号排列的一维序列需要将其映射回计算网格才能进行空间分析。这需要结合网格文件fort.14中的拓扑信息。网格映射关键步骤解析fort.14文件获取节点坐标和元素连接关系将节点数据插值到规则网格可选构建时空数据立方体x, y, time, valueimport matplotlib.tri as tri def create_triangulation(fort14_path): with open(fort14_path, r) as f: f.readline() # 跳过标题行 node_count, elem_count map(int, f.readline().split()[:2]) # 读取节点坐标 nodes [] for _ in range(node_count): parts f.readline().split() nodes.append((float(parts[1]), float(parts[2]))) # 读取三角形元素 triangles [] for _ in range(elem_count): parts f.readline().split() triangles.append([int(parts[2])-1, int(parts[3])-1, int(parts[4])-1]) x np.array([n[0] for n in nodes]) y np.array([n[1] for n in nodes]) return tri.Triangulation(x, y, triangles)时空数据重构示例def reshape_to_space_time(df, time_steps): 将DataFrame转换为时空矩阵 node_count df[node].nunique() space_time np.zeros((len(time_steps), node_count)) for i, t in enumerate(time_steps): subset df[df[time] t] space_time[i,:] subset[value].values return space_time提示对于大型数据集建议使用xarray库处理多维数组它提供了更高效的内存管理和便捷的切片操作。3. 专业级海洋数据可视化技术有了结构化的时空数据我们可以创建各种专业图表。以下是几种典型可视化场景的实现方法。3.1 风暴潮增水分布图import cartopy.crs as ccrs import cartopy.feature as cfeature def plot_water_level(triangulation, values, bbox, title): fig plt.figure(figsize(12, 8)) ax fig.add_subplot(111, projectionccrs.PlateCarree()) # 设置地图范围 ax.set_extent(bbox, crsccrs.PlateCarree()) # 添加地图要素 ax.add_feature(cfeature.LAND, facecolorlightgray) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:) # 绘制水位等高线 levels np.linspace(values.min(), values.max(), 20) contour ax.tricontourf(triangulation, values, levelslevels, cmapviridis, transformccrs.PlateCarree()) # 添加色标 cbar plt.colorbar(contour, axax, orientationvertical, shrink0.7) cbar.set_label(Water Level (m)) # 添加标题和网格 ax.set_title(title, fontsize14) ax.gridlines(draw_labelsTrue, linestyle--) return fig3.2 流场矢量动态可视化fort.64文件包含流速的U/V分量适合制作动态流场图from matplotlib.animation import FuncAnimation def create_flow_animation(triangulation, u_series, v_series, times): fig, ax plt.subplots(figsize(12, 8)) # 初始帧 quiver ax.quiver(triangulation.x, triangulation.y, u_series[0], v_series[0], scale20, colorblue) title ax.set_title(fFlow Field at t{times[0]:.1f}s) def update(frame): quiver.set_UVC(u_series[frame], v_series[frame]) title.set_text(fFlow Field at t{times[frame]:.1f}s) return quiver, title ani FuncAnimation(fig, update, frameslen(times), interval100, blitTrue) return ani3.3 波浪参数联合分析图SWAN输出通常包含多个波浪参数可以制作联合图表展示它们的关系def plot_wave_parameters(hs, dir, tps, triangulation): fig, (ax1, ax2, ax3) plt.subplots(1, 3, figsize(18, 6)) # 有效波高 cont1 ax1.tricontourf(triangulation, hs, levels20, cmaprainbow) fig.colorbar(cont1, axax1, labelSignificant Wave Height (m)) ax1.set_title(Significant Wave Height) # 波向玫瑰图 ax2 fig.add_subplot(132, projectionpolar) dir_rad np.deg2rad(dir) hist, bins np.histogram(dir_rad, bins36, range(0, 2*np.pi)) width (2*np.pi) / 36 bars ax2.bar(bins[:-1], hist, widthwidth, bottom0.0) ax2.set_title(Wave Direction Distribution) # 波高-周期散点图 sc ax3.scatter(hs, tps, cdir, cmaphsv, alpha0.6) fig.colorbar(sc, axax3, labelWave Direction (deg)) ax3.set_xlabel(Significant Wave Height (m)) ax3.set_ylabel(Peak Period (s)) ax3.set_title(Hs-Tp Relationship) plt.tight_layout() return fig4. 高级技巧与性能优化处理大规模海洋模型输出时性能和内存管理至关重要。以下是几个实用技巧NetCDF格式转换将ASCII输出转换为NetCDF可以大幅提高后续处理效率import xarray as xr def convert_to_netcdf(df, triangulation, output_path): # 创建xarray Dataset ds xr.Dataset( { water_level: ([time, node], df.pivot(indextime, columnsnode, valuesvalue).values) }, coords{ time: df[time].unique(), lon: (node, triangulation.x), lat: (node, triangulation.y) } ) # 添加元数据 ds.attrs[description] ADCIRC water level output ds[water_level].attrs[units] meters # 保存为NetCDF ds.to_netcdf(output_path)分块处理大型文件对于超出内存的数据可以使用分块处理策略def process_large_file(filename, chunk_size1000000): results [] for chunk in pd.read_csv(filename, chunksizechunk_size, delim_whitespaceTrue, headerNone): # 处理每个数据块 processed process_chunk(chunk) results.append(processed) # 合并结果 return pd.concat(results)并行计算加速利用多核处理器加速数据处理from multiprocessing import Pool def parallel_process_time_steps(time_steps, func, n_workers4): with Pool(n_workers) as pool: results pool.map(func, time_steps) return np.stack(results)在完成一系列风暴潮模拟项目后我发现最耗时的往往不是模型计算本身而是后续的数据分析和可视化阶段。通过建立标准化的处理流程和可复用的代码库可以节省大量时间。例如将常用可视化函数封装成类可以方便地在不同项目间共享class AdcircVisualizer: def __init__(self, grid_file): self.tri self._load_grid(grid_file) self.crs ccrs.PlateCarree() def _load_grid(self, path): # 实现网格加载逻辑 pass def plot_field(self, values, **kwargs): # 实现通用场绘图逻辑 pass def animate_field(self, time_series, **kwargs): # 实现动画生成逻辑 pass