【Delft3D FM数据后处理系列】2. 绘制模型地形图

        平时在学习Delft3D FM的时候，发现网上很多资料主要是讲如何运行模型，鲜有讲如何对结果进行可视化的。虽然Delft3D有quickplot工具，但对于科研绘图来说还远远不够，个人感觉只适合于简单的查看数据。此外，Deltares的OpenEarthTools虽然也有处理Delft3D数据的程序，但是那个我个人目前还没成功使用过，程序报错后很难处理。因此，我尝试自己写处理Delft3D数据的程序，方便后面出图与修改。主要用到的工具是Matlab和Surfer。

        注意：这个系列面向的是刚接触Delft3D FM的朋友，文章用到的都是十分基础的命令及操作，如果大家对Delft3D之类的东西已经十分熟悉的话请降低期待值观看。

        这一节主要介绍如何绘制模型地形图。在论文的资料及方法这一章节中通常需要绘制模型地形图。我见过的模型地形图有3类，第一类是将模型网格与地形结合起来的（如图1），第二类是直接输出模型地形的（如图2），第三类是模型地形中还加上一些等值线的。下面分别介绍这三类该怎么画。

第一类：需要用到Delft3D FM的地形文件（通常为xyz格式），方法跟第一篇文章类似，也是用matlab的patch函数。

代码如下：

% 用来处理Delft3D输出的Map.nc文件，画网格地形
clear,clc;

%% 读取数据
% ncdisp('Boluo2map.nc')
mapfilename = 'Boluo2map.nc';
lon = ncread(mapfilename,'mesh2d_node_x');
lat = ncread(mapfilename,'mesh2d_node_y');
FaceConnect = ncread(mapfilename,'mesh2d_face_nodes');

sn1 = find(isnan(FaceConnect(4,:)));
sn2 = find(~isnan(FaceConnect(4,:)));

data = dlmread('2016_depth.xyz');
F = scatteredInterpolant(data(:,1),data(:,2),data(:,3));
z = F(lon,lat);

%% 画图

patch('Faces',FaceConnect(1:4,sn2)','Vertices',[lon lat],'CData',z,...
    'EdgeColor','interp','FaceColor','none')
patch('Faces',FaceConnect(1:3,sn1)','Vertices',[lon lat],'CData',z,...
    'EdgeColor','interp','FaceColor','none')

结果如下：

第二类：方法也是一样，用patch函数（其实map.nc文件里面数据名称带mesh2d_face的都可以用patch函数去画）

clear,clc;
addpath E:\Delft3D

%% 读取数据并确定范围
% ncdisp('PRDModel2016_Dec_001_merged_map.nc')
mapfilename = 'PRDModel2016_Dec_001_merged_map.nc';
lon = ncread(mapfilename,'mesh2d_node_x');
lat = ncread(mapfilename,'mesh2d_node_y');
FaceConnect = ncread(mapfilename,'mesh2d_face_nodes');
t = ncread(mapfilename,'time');
bl = ncread(mapfilename,'mesh2d_flowelem_bl');


sn1 = find(isnan(FaceConnect(4,:)));
sn2 = find(~isnan(FaceConnect(4,:)));

patch(lon(FaceConnect(1:4,sn2)),lat(FaceConnect(1:4,sn2)),bl(sn2))
hold on
patch(lon(FaceConnect(1:3,sn1)),lat(FaceConnect(1:3,sn1)),bl(sn1))
h = colorbar;

第三类：patch方法在画mesh2d_face类型的图时确实很方便，但是也有个缺陷，就是不能画等值线图（也可能可以，但我还不会）。所以如果需要画带等值线的网格地形图时，可以用到surfer软件进行绘图，此时仅用到地形数据文件（*.xyz）。

这里有可能有读者不太懂白化的意思，我的理解是，surfer会先将xyz数据插值成一个3D数组，也就是说在空间上每个点都会被插值出一个水深值，包括陆地上的点。但实际上陆地上的点是不在水动力模型范围内的，此时就需要把陆地上的点值设置成空值，即白化。白化需要*.bln文件，这个可以自己对着地图抠或者找师兄师姐要，后续也可能会出bln文件的制作方法。

关于surfer图片的调整，不熟悉的读者可以私信我交流（虽然我也不一定懂）。

以上就是本节的内容啦！大家可以在评论区踊跃交流！感谢你读到这里~