简介DHI MATLAB工具箱面向水利、海洋及环境领域需要处理MIKE系列模型输出数据的科研与工程人员提供了一套可直接运行的MATLAB脚本用于读写和转换DHI格式文件如dfs0、dfs2、dfsu等尤其适合有MATLAB基础、希望自动化批量处理DHI数据的开发者。压缩包共108个文件以m脚本为主体辅以C源码和mexw64/mexw32编译组件配合bat构建脚本完成Mex接口的重建与工具链配置同时包含示例数据、PDF与docx说明文档、license许可及res1d/res11/mesh等模型网格文件。整体仅5.41MB轻量易部署目前已有116人学习。借助内置的水位、浓度等样例数据读者可快速验证读写脚本的正确性并在此基础上扩展DHI结果的后处理流程通过bat与c/cs文件还能理解从源码编译到调用的完整链路适合需要在MATLAB中深度集成DHI工具链的中高级用户。1. 先认识DHI文件为什么MATLAB原生读不了做水利、海洋、环境模拟这行的对DHI这个名字应该都不陌生。丹麦水力研究所DHI出品的MIKE系列软件从MIKE 11、MIKE 21到MIKE 3几乎算是水动力模拟领域事实上的标准工具。这些软件跑完模拟之后输出的结果文件就是DHI自己定义的二进制格式常见的包括.dfs0时间序列、.dfs1一维断面、.dfs2二维平面场、.dfs3三维体数据还有MIKE 11的.res11、MIKE 21的.res21这些结果文件。问题就出在这里。这类文件是高度自定义的二进制结构MATLAB自带的fopen、fread根本不知道里面存了什么。你把.dfs2文件拖进MATLAB里直接读得到的是一堆乱码数字完全没法用。以前我见过不少同事的做法是先在MIKE软件里把结果导出成.asc或者.txt再拿回MATLAB处理。这种方式有两个硬伤——第一MIKE软件本身价格不菲不可能每个处理数据的人都装一套第二巨型模型的输出文件动辄几个GB导成文本文件之后体积膨胀好几倍读写慢得让人怀疑人生。所以DHI官方和第三方社区折腾出了好几套方案。DHI官方在GitHub上有开源的DHI MATLAB Toolbox核心目的就是让MATLAB用户能直接读写.dfs0、.dfs1、.dfs2、.dfs3以及res系列文件不需要安装MIKE软件也不需要中间转格式。这个工具箱本质上是把MATLAB和DHI底层文件格式之间搭了一座桥把所有二进制解析细节封装成了一套看起来像普通MATLAB函数的接口。这个标题里说的DHI MATLAB工具箱包含用于处理DHI文件的MATLAB脚本指的就是这类工具。我在实际项目里用这套东西处理过几十个GB的二维水动力模拟结果深有体会——搞懂它的使用逻辑能省掉大量重复劳动。这篇文章就围绕这套工具箱把它的文件格式背景、脚本功能、实际用法和踩坑经历一次讲清楚。2. DHI文件结构的核心逻辑network、item和time step在动手写代码之前有个概念必须先建立起来否则你看着工具箱的接口函数根本不明白为什么参数要这么传。DHI的文件格式虽然二进制层面很复杂但在逻辑上抽象成了三个层级network网格/网络、item变量项和time step时间步。network描述的是数据的空间结构。对于.dfs2文件来说network就是一张矩形网格包含网格原点坐标、网格尺寸、行列数、投影坐标系信息对于.dfs1来说network是一组沿断面的点.dfs3则是三维体网格。item描述的是文件里存了哪些物理量比如水位、流速U分量、流速V分量、含沙量、温度等。每个item有自己的名字、数据类型和单位。time step就好理解了就是每个时间帧对应的数据。你拿到一个.dfs2文件本质上就是一堆二维数组的集合每个item在每个time step上都对应一个二维数组行×列。这跟MATLAB里处理一个普通矩阵的思维方式完全不同因为数据在磁盘上不是连续存放的而是按item × time的组织方式分块存储的。% 使用DHI工具箱读取文件的基本流程 % 第一步打开文件获取文件的全局信息 dfs2_file dfsTSO(open, flow_result.dfs2); % 第二步读取item信息 item_info dfs2_file.ItemInfo; % 返回一个结构体数组 % 每个元素包含Name变量名、DataType、Quantity含单位 % 第三步读取某个item在某个time step下的数据 data dfs2_file.readItemTimeStep(item_index, time_step_index); % 返回的data有两个字段 % data.data —— 实际的二维数据矩阵 % data.time —— 当前时间步的时间戳这个工具箱把底层一切复杂的东西都藏起来了。你不需要知道文件里哪些字节是网格坐标、哪些字节是缩放因子只需要按照打开—查询—读取—关闭这个流程走就行。但要注意一个关键点读取的索引从1开始这是工具箱为了适应MATLAB习惯做的处理跟很多从0开始的底层接口不一样新手经常在这里踩坑。还有一个很重要的设计readItemTimeStep每次只读取一个item的一个时间步。这意味着如果文件里有10个item、100个时间步你要完整读完就得调用1000次。对于大型文件这个读取效率值得单独说道说道后面我会详细展开。3. 工具箱的核心脚本拆解每个函数是干什么的DHI MATLAB工具箱包含的脚本数量不少但真正高频使用的核心函数掰着手指头数也就那么几个。我按功能把它们分成四类帮大家快速建立认知框架。3.1 文件打开与关闭类dfs0 dfsTSO(open, rainfall.dfs0); % 打开dfs0时序文件 dfs2 dfsTSO(open, depth.dfs2); % 打开dfs2网格文件 % 用完必须关闭 close(dfs0); close(dfs2);注意一个细节dfsTSO这个类叫Time Series Object但实际操作中对dfs1、dfs2、dfs3文件也是用它。工具箱内部会根据文件扩展名自动判断文件类型所以你不需要手动区分。文件打开之后所有后续操作都基于返回的这个对象进行。3.2 信息查询类% 获取item列表 items dfs2.ItemInfo; for i 1:length(items) fprintf(Item %d: %s, 单位: %s\n, i, items(i).Name, items(i).Quantity); end % 获取时间轴信息 time_axis dfs2.TimeAxis; % 包含 startTime起始时间、TimeStep时间步长、NumberOfTimeSteps总步数 % 获取网格/network信息 grid_info dfs2.GridInfo; % 包含 originX、originY、gridStepX、gridStepY、numberOfGridPointsX、numberOfGridPointsY这些查询函数返回的都是MATLAB结构体数组字段命名比较直观基本看名字就能理解含义。这也是工具箱做得好的一点——虽然底层是C实现的DHI文件处理库但在MATLAB这层封装得很MATLAB化。3.3 数据读取类数据读取是使用频率最高、也是最容易出问题的一环。前面提到过的基本读取方式之外工具箱还提供了一个更高效的批量读取接口% 一次性读取某个item的所有时间步 all_data dfs2.readItem(item_index);readItem和readItemTimeStep的区别在于前者一次性分配好内存、连续读取全部时间步的数据适合中等规模的数据集后者每次只读一帧适合大文件按需处理。这个选择直接影响内存占用和运行速度我在第5部分会给出具体的数据对比。3.4 数据写入与创建类工具箱不只是能读还能写。你可以用dfsTSO(create, ...)创建一个新文件然后向里面写入数据。这在做模型结果后处理、格式转换时非常有用。% 创建一个新的dfs2文件 dfs2_new dfsTSO(create, output.dfs2, 文件标题, ... grid, grid_info, ... items, item_def, ... time, time_axis); % 写入数据 dfs2_new.writeItemTimeStep(time_step_index, item_index, data_matrix); % 关闭文件 close(dfs2_new);需要提醒的是writeItemTimeStep的写入顺序并不要求按时间步递增你可以跳着写。但必须保证每个item的每个时间步都被写入否则生成的文件在这个时间步上就是空数据用MIKE软件打开时会报错。这个坑我踩过后面细说。3.5 坐标系与地理参考的坑DHI工具箱处理的是带地理参考的网格数据。.dfs2文件的网格信息里包含投影坐标系定义常见的有UTM投影通用横轴墨卡托投影、经纬度坐标等。工具箱在打开文件时并不会自动检查坐标系是否匹配——如果你同时打开两个不同投影系的文件做计算MATLAB不会报错但算出来的结果在空间上完全对不上。一个实际案例我曾经处理一个项目需要用MIKE模拟的水位结果和实测的潮位站数据做对比。实测数据是经纬度坐标WGS84模拟结果是UTM Zone 50N投影坐标。直接读取模拟结果的网格数据再和自己手写的经纬度距离公式计算站点的位置结果差了十几公里。查了半天才发现是坐标转换没有做。正确的做法是先用工具箱读取网格信息然后用mfwdtran或projfwdMapping Toolbox把经纬度转换成UTM坐标再去网格里定位对应的行列号。如果用的是新版工具箱还可以直接用mapproj函数处理方便很多。4. 从零搭建一个水动力结果提取流程一个可落地的实例前面铺垫了那么多概念和函数现在用一个完整的实例把这些串起来。假设场景是我有一个MIKE 21模拟输出的.dfs2文件包含24小时的水位场结果时间步长10分钟共145个时间步网格是500行×400列。现在需要提取某个坐标点处的逐时水位序列并画成过程线。4.1 步骤一读取网格信息并定位目标点% 打开文件 dfs2 dfsTSO(open, storm_surge.dfs2); % 获取网格参数 g dfs2.GridInfo; fprintf(Origin: (%f, %f)\n, g.originX, g.originY); fprintf(Grid size: %f x %f\n, g.gridStepX, g.gridStepY); fprintf(Grid points: %d x %d\n, g.numberOfGridPointsX, g.numberOfGridPointsY); % 目标点经纬度/UTM坐标 target_x 580000; % UTM坐标 东向 target_y 4270000; % UTM坐标 北向 % 计算目标点在网格中的行列号向上取整 col round((target_x - g.originX) / g.gridStepX) 1; row round((target_y - g.originY) / g.gridStepY) 1;这里有个小细节网格列的编号是X方向东向行的编号是Y方向北向和MATLAB矩阵的下标行在前、列在后恰好相反。所以读取数据时矩阵索引应该是data(row, col)而不是data(col, row)。方向搞反的话提取出来的位置会偏移在异方差地形区域结果完全不可用。4.2 步骤二批量读取时间序列确定了目标点之后循环读取每个时间步的数据% 找到水位对应的item索引 items dfs2.ItemInfo; % 假设水位在item 1 num_steps dfs2.TimeAxis.NumberOfTimeSteps; water_level NaN(num_steps, 1); time_vector NaT(num_steps, 1); for t 1:num_steps frame dfs2.readItemTimeStep(1, t); water_level(t) frame.data(row, col); time_vector(t) frame.time; end % 关闭文件 close(dfs2);对于145个时间步这种循环速度很快几乎是瞬间完成。但如果遇到几千个时间步的大型文件这种逐个读取的方式效率就会下降。替代方案是用readItem一次读入整个itemall_data dfs2.readItem(1); % 此时all_data.data是一个 num_steps × row × col 的三维数组 water_level squeeze(all_data.data(:, row, col));这段代码在数据量较大时能感觉到明显的速度差异。readItem底层是连续读入的内存块而readItemTimeStep需要每次做文件定位和读取准备工作。耗时对比上读取一个2000步的dfs2文件readItem大概只要几秒而循环读取可能要几十秒甚至更久。4.3 步骤三绘图和输出figure; plot(time_vector, water_level, b-, LineWidth, 1.5); datetick(x, HH:MM, keepticks); xlabel(时间); ylabel(水位 (m)); title(目标点逐时水位变化); grid on; % 输出到csv方便后续使用 T table(time_vector, water_level); writetable(T, water_level_at_station.txt, Delimiter, tab);整个流程最耗时的环节其实不在提取而在初始的坐标转换和位置定位。如果工具箱使用频繁建议把这些逻辑封装成函数输入文件路径和经纬度坐标直接返回该位置的时间序列以后所有项目都能复用。5. 实测中的那些坑版本兼容、单位换算与批处理陷阱5.1 版本兼容性——旧代码可能跑不通新文件DHI的文件格式在历次版本中有过调整工具箱本身也在持续更新。一个常见的问题是用旧版工具箱打开新版MIKE软件输出的文件可能出现文件打不开、item字段读取异常等情况。解决思路是版本匹配。DHI官方在GitHub的仓库里会标注支持的MATLAB版本和对应的MIKE版本。我个人的经验是能装新版工具箱就不用旧版因为新版本不仅支持新格式还修复了很多老版本的内存泄漏问题。但要注意新版工具箱对MATLAB的最低版本有要求需要提前确认。5.2 单位陷阱——读出来的数据跟你想象的可能不一样DFS文件里每个item都有自己的Quantity定义包含基本单位和导出单位。工具箱读取时默认返回的是文件里存的原始单位不一定是你要的单位。比如MIKE在水深文件中可能以m为单位但在流速文件中可能是m/s有时也可能是cm/s取决于建模时怎么设置的。有一个项目我印象很深模拟工程师在MIKE里设置的流速单位是cm/s读取数据后我直接用这个数值去计算流量发现结果比经验估算大了100倍。排查半天才查到单位问题。所以每次读取数据之后一定要先检查item_info.Quantity和单位字段确认无误再进入后续计算。% 检查item单位避免单位不一致问题 for i 1:length(items) disp([Item: , items(i).Name, , unit: , items(i).Quantity]); end5.3 批处理中的文件句柄泄漏如果要一次性处理成百上千个文件一个很隐蔽的坑是文件句柄没有正确释放。dfsTSO对象在MATLAB中属于外部库对象垃圾回收机制不一定能及时释放底层文件句柄。我在一台服务器上批量处理数据跑到几百个文件之后突然报Too many open files错误查了好久才发现是对象没有显式关闭。解决方法是每处理完一个文件立刻close或者用try...catch...end确保出错时也能关闭文件file_list dir(results/*.dfs2); for k 1:length(file_list) dfs2 dfsTSO(open, fullfile(file_list(k).folder, file_list(k).name)); try % 处理数据 data dfs2.readItem(1); % ... 其他操作 ... catch ME warning(处理文件 %s 时出错: %s, file_list(k).name, ME.message); end close(dfs2); % 确保每次循环都关闭文件 end5.4 缺失值的表示方式很多DHI文件里会有缺测值land points这些位置的数据在文件里存储的是特殊的大负数通常是-1.e-30或-999.0。直接拿这些数值做分析会污染结果。工具箱的文档里会给出缺测值的具体数值通常可以这样处理data(data -999) NaN; % 把缺测值替换为NaN % 但注意不同版本的DHI文件缺测值可能不同不要盲目套用 % 可以先查看文件的DeleteValue属性或者直接分析数据分布我在处理水深数据时发现文件中除了水域格点外陆地格点存的都是-1.e-30这种接近0的负值如果不处理最小值、均值等统计量全都会被带偏。用NaN替换后后续分析和可视化才正常。建议在读取每个文件后先打印一下数据的极值和分布看看有没有异常值形成习惯。5.5 网格方向的错误假设DFS文件网格的方向和MATLAB矩阵的默认方向不一样这是最基础但也最容易反复犯的错。在DFS2中第1行是Y坐标最大值的那一行即最北行第1列是X坐标最小值的那一列最西列。而MATLAB的矩阵arr(1,1)是左上角元素。如果直接把data当作矩阵画imagesc图像上下就是颠倒的。正确的绘图方式imagesc(g.originX (0:g.numberOfGridPointsX-1)*g.gridStepX, ... g.originY (0:g.numberOfGridPointsY-1)*g.gridStepY, ... data); axis xy; % 注意这个命令让Y轴朝上否则图像是倒着的 set(gca, YDir, normal); % 与axis xy同理但直接控制Y轴方向这个细节虽然小但直接影响图和实际地理位置的对应关系。我第一次画结果图的时候没注意图是倒置的整个海域南北方向反了过来如果不是后续叠加岸线时发现不对就拿着错误的结果去汇报了。6. 在大型项目里我是怎么组合这套工具箱的效率优化与二次封装当数据处理规模大起来之后单独调用几个核心函数就有点不够用了。我在一个沿海城市风暴潮风险评估项目里要处理的是一整年的MIKE 21模拟结果每天1个文件每个文件2.5GB左右共有365个文件。所有文件包含12个item、1440个时间步10分钟间隔。如果用最基础的循环读取方式一个文件处理就要半个多小时整个数据集处理完不现实。最终的做法是做了一层二次封装核心是解决两个问题只读需要的数据和并行化处理。6.1 按需读取不读整个文件很多模拟输出的dfs2文件里水深、初始水位这些item基本是常数只有少数几个item是变化的。按需读取可以避免不必要的数据搬运。比如我们项目里真正需要的是水位、流速U、流速V三个item那就只读取这三个其他item直接跳过。% 建立item名称到索引的映射只取出感兴趣的变量 function idx findItemIndex(items, targetName) idx []; for i 1:length(items) if strcmpi(items(i).Name, targetName) idx [idx, i]; end end end % 主处理逻辑只循环需要的item不遍历全部6.2 用parfor并行处理多个文件由于每个文件相互独立多个文件之间的处理流程完全可以用MATLAB并行计算工具箱来加速parpool(local, 4); % 根据自己的核心数调整 parfor k 1:length(file_list) dfs2 dfsTSO(open, fullfile(file_list(k).folder, file_list(k).name)); try % 只读取关键变量 items dfs2.ItemInfo; idx_wl findItemIndex(items, Surface elevation); idx_u findItemIndex(items, U velocity); idx_v findItemIndex(items, V velocity); data_wl dfs2.readItem(idx_wl); data_u dfs2.readItem(idx_u); data_v dfs2.readItem(idx_v); % 处理数据比如提取最大值、累积量等 % ... % 保存结果到单独的mat文件 save(sprintf(processed_%d.mat, k), data_wl, data_u, data_v); catch ME warning(处理文件 %s 出错: %s, file_list(k).name, ME.message); end close(dfs2); end我从4核并行到8核365个文件总共花了大概7个小时跑完比原来单线程的预期时间缩短了三分之二。如果手动处理这活儿可能要加班干一整周。6.3 中间结果用MAT文件中转处理过程中建议把中间结果保存成.mat文件而不是再次用dfs2格式保存。原因很简单.mat文件读取速度比重新解析dfs2快得多尤其是处理完一轮数据后后续的统计分析和可视化都用.mat做中转能明显提升效率。6.4 封装一个完整的读取函数最终我把整套流程封装成了一个函数输入是文件路径和坐标点输出是该点的全部变量时间序列。这样团队里其他人不需要了解dfs2格式的细节拿着经纬度坐标就能直接取得数据。这也是工具类代码真正发挥价值的地方——把复杂的事情封装成简单接口。function result extractPointData(file_path, target_x, target_y) dfs2 dfsTSO(open, file_path); try g dfs2.GridInfo; col round((target_x - g.originX) / g.gridStepX) 1; row round((target_y - g.originY) / g.gridStepY) 1; items dfs2.ItemInfo; num_items length(items); num_steps dfs2.TimeAxis.NumberOfTimeSteps; result struct(); result.time NaT(num_steps, 1); for i 1:num_items result.(matlab.lang.makeValidName(items(i).Name)) NaN(num_steps, 1); end for t 1:num_steps for i 1:num_items frame dfs2.readItemTimeStep(i, t); name matlab.lang.makeValidName(items(i).Name); result.(name)(t) frame.data(row, col); end result.time(t) frame.time; end catch ME rethrow(ME); end close(dfs2); end注意一个细节matlab.lang.makeValidName用来处理item名称中可能存在的空格和特殊字符避免字段名不合法导致赋值报错。7. 能扩展的方向从读取到写回、从点数据到场数据工具箱的价值不只是读还有一个容易被低估的能力是写回和创建文件。如果想把MATLAB里的计算成果比如做数据同化之后的修正场、或者用遥感反演成果替代某个区域的结果重新写成一个合法的dfs2文件拿给MIKE软件继续做下一阶段模拟这个工作流是完全能走通的。我在项目里做过一次类似的操作用MATLAB跑完粒子追踪后把每个时间步的粒子浓度场写成dfs2文件直接送给MIKE的Advection-Dispersion模块继续算。整个过程里最关键的是创建文件时把网格信息、item信息、时间轴完全照搬原文件不要自己重新定义否则输出文件的网格范围和原点对不上模型直接罢工。% 从已有文件复制定义到新文件 template dfsTSO(open, template.dfs2); grid template.GridInfo; items template.ItemInfo; time template.TimeAxis; % 创建新文件沿用原文件空间和时间定义 new_file dfsTSO(create, lagrangian_output.dfs2, ... Particle concentration, ... grid, grid, items, items, time, time); % 写入处理结果 for t 1:time.NumberOfTimeSteps new_file.writeItemTimeStep(t, 1, concentration_field(:,:,t)); end close(new_file); close(template);工具箱还有一个配套的dfsutil工具包包含若干命令行小工具比如在多个文件中批量查询item定义、或者压缩转存大规模文件。不过这些工具的使用场景比较窄等真正遇到需求时去翻工具箱自带的help文件也不迟。最后说一个个人体会这类专门格式的读取工具最忌讳的就是拿过来直接用而不做测试。每拿到一个新的DHI文件第一次总是先用工具箱读一下item定义和网格信息打印出来确认和自己预期的数据结构一致再跑正式的批处理流程。这套先探路、再赶路的习惯帮我避开了很多因为文件格式版本、坐标系设置、单位定义差异带来的坑。如果你也经常跟MIKE数据打交道建议在这个工具箱上花点时间把基础流程搞清楚后面真的能省出好几周的功夫。本文还有配套的精品资源点击获取