MATLAB基于全球潮汐模型TPXO绘制某经纬度海域同潮图

1.全球潮汐模型TPXO。网站：OSUTPXOTide Models(www.tpxo.net/home)

2.需要下载MATLAB工具包并设置路径：TMD Matlab Toolbox。网址：https://ww2.mathworks.cn/matlabcentral/fileexchange/75599-tide-model-driver-tmd-version-2-5-toolbox-for-matlab

3.绘制同潮图

（1）.提取调和常数

（2）.振幅--颜色填充 pcolor

（3）.迟角—等值线 contour

4.成果图展示（以中国南海K1、M2分潮为例）

5.中国南海K1分潮同潮图：

clc; clear; close all;
% addpath('D:\matlab\toolbox\TMD_Matlab_Toolbox\TMD')
Model = 'C:\Users\Mayn\Documents\MATLAB\TPXO_model\TPXO_model\Model_Ind_2016';
% 设定高分辨率网格
lon_grid = 105:0.01:125; 
lat_grid = 3:0.01:25; 
[lon_mesh, lat_mesh] = meshgrid(lon_grid, lat_grid); 
% 提取 K1 分潮数据
type = 'z'; 
Cid_K1 = 5; % 对应 K1 分潮
[amp_K1, Gph_K1, ~, ~] =tmd_extract_HC(Model, lat_mesh(:), lon_mesh(:), type, Cid_K1);
amp_K1 = reshape(amp_K1, size(lon_mesh));
Gph_K1 = reshape(Gph_K1, size(lat_mesh));
% 插值
F_amp_K1 = scatteredInterpolant(lon_mesh(:), lat_mesh(:), amp_K1(:), 'linear', 'none');
F_Gph_K1 = scatteredInterpolant(lon_mesh(:), lat_mesh(:), Gph_K1(:), 'linear', 'none');
% 高分辨率经纬度网格
[lon_fine, lat_fine] = meshgrid(105:0.02:125, 3:0.02:25);
% 插值
amp_K1_fine = F_amp_K1(lon_fine, lat_fine);
Gph_K1_fine = F_Gph_K1(lon_fine, lat_fine);
% 将 NaN 值替换为最小值
amp_K1_fine(isnan(amp_K1_fine)) = min(amp_K1_fine(:));
% 绘制 K1 分潮图
figure;
m_proj('Mercator', 'lon', [105 125], 'lat', [3 25]); 
hold on;
% 绘制 K1 振幅
m_pcolor(lon_fine, lat_fine, amp_K1_fine);
shading interp; 
colormap(jet); % 使用 jet 颜色映射
colorbar;
caxis([0 0.5]);
% 绘制 K1 迟角等值线
[C_K1, h_K1] = m_contour(lon_fine, lat_fine, Gph_K1_fine, 0:30:360, 'k');
clabel(C_K1, h_K1, 'fontsize', 8, 'color', 'k');
% 添加海岸线并设置绿色陆地
m_gshhs_h('patch', [0.63 0.8 0.32], 'edgecolor', 'k'); 
m_grid('box', 'on', 'tickdir', 'out'); 
m_gshhs('ib', 'color', 'k', 'linestyle','--','linewidth', 2.5);
m_gshhs('fr', 'color', 'b', 'linewidth', 0.6);
title('南海K1分潮同潮图', 'fontsize', 10);
xlabel('经度 (°E)');
ylabel('纬度 (°N)');
hold off;
% 保存为 PNG 文件
print(gcf, '-dpng', 'K1.png');

中国南海区域的M2分潮同潮图：

clc; clear; close all;
% addpath('D:\matlab\toolbox\TMD_Matlab_Toolbox\TMD')
Model = 'C:\Users\Mayn\Documents\MATLAB\TPXO_model\TPXO_model\Model_Ind_2016';
% 设定高分辨率网格
lon_grid = 105:0.01:125; 
lat_grid = 3:0.01:25; 
[lon_mesh, lat_mesh] = meshgrid(lon_grid, lat_grid); 
% 提取 M2 分潮数据
type = 'z'; 
Cid_M2 = 1; % 对应 M2 分潮
[amp_M2, Gph_M2, ~, ~] = tmd_extract_HC(Model, lat_mesh(:), lon_mesh(:), type, Cid_M2);
amp_M2 = reshape(amp_M2, size(lon_mesh));
Gph_M2 = reshape(Gph_M2, size(lat_mesh));
% 插值
F_amp_M2 = scatteredInterpolant(lon_mesh(:), lat_mesh(:), amp_M2(:), 'linear', 'none');
F_Gph_M2 = scatteredInterpolant(lon_mesh(:), lat_mesh(:), Gph_M2(:), 'linear', 'none');
% 高分辨率经纬度网格
[lon_fine, lat_fine] = meshgrid(105:0.02:125, 3:0.02:25);
% 插值
amp_M2_fine = F_amp_M2(lon_fine, lat_fine);
Gph_M2_fine = F_Gph_M2(lon_fine, lat_fine);
% 将 NaN 值替换为最小值
amp_M2_fine(isnan(amp_M2_fine)) = min(amp_M2_fine(:));
% 绘制 M2 分潮图
figure;
m_proj('Mercator', 'lon', [105 125], 'lat', [3 25]); 
hold on;
% 绘制 M2 振幅
m_pcolor(lon_fine, lat_fine, amp_M2_fine);
shading interp; 
colormap(jet); % 使用 jet 颜色映射
colorbar;
caxis([0 0.8]);
% 绘制 M2 迟角等值线
[C_M2, h_M2] = m_contour(lon_fine, lat_fine, Gph_M2_fine, 0:30:360, 'k');
clabel(C_M2, h_M2, 'fontsize', 8, 'color', 'k');
% 添加海岸线并设置绿色陆地
m_gshhs_h('patch', [0.63 0.8 0.32], 'edgecolor', 'k'); 
m_grid('box', 'on', 'tickdir', 'out'); 
m_gshhs('ib', 'color', 'k', 'linestyle','--','linewidth', 2.5);
m_gshhs('fr', 'color', 'b', 'linewidth', 0.8);
title('南海M2分潮同潮图', 'fontsize', 10);
xlabel('经度 (°E)');
ylabel('纬度 (°N)');
hold off;
print(gcf,'-dpng','M2.png')

2022年1月份站点流速数据（8个分潮，长度1年，时间间隔1h）数据图：

% addpath('"D:\matlab\toolbox\old_tmd_toolbox"')
% 设定 TPXO 模型路径
Model = 'C:\Users\Mayn\Documents\MATLAB\TPXO_model\TPXO_model\Model_Ind_2016';;
% 设定开始和结束时间
t0 = datetime(2022, 1, 1, 0, 0, 0); % 开始时间
t1 = datetime(2022, 12, 31, 23, 0, 0); % 结束时间
SDtime = datenum(t0):1/24:datenum(t1); % 转换为 MATLAB 序列时间
% 站点信息（南海某站点）
lat = 22.7; % 站点纬度
lon = 113.7; % 站点经度
% 8 个主要分潮（M2, S2, N2, K1, O1, P1, Q1, K2）
Cid = [1,2,3,4,5,6,7,8]; 
% 计算 U 和 V 方向的潮流流速
U = tide_pred(Model, SDtime, lat, lon, 'u', Cid);
V = tide_pred(Model, SDtime, lat, lon, 'v', Cid);
% 计算流速大小
velocity = sqrt(U.^2 + V.^2);
% 选取第一个月的数据
t_start = datenum(datetime(2022, 1, 1, 0, 0, 0));
t_end = datenum(datetime(2022, 1, 31, 23, 0, 0));
idx = (SDtime >= t_start) & (SDtime <= t_end);
% 绘制潮流流速时间序列
figure;
subplot(211);
plot(datetime(SDtime(idx), 'ConvertFrom', 'datenum'), V(idx), 'b', 'LineWidth', 1.5);
xlabel('Time (days)');
ylabel('V (m/s)');
title('南海站点 潮流流速V时间序列 (2022年1月)','fontsize',8);
grid on;
datetick('x', 'dd', 'keeplimits'); % X 轴显示日期
subplot(212);
plot(datetime(SDtime(idx), 'ConvertFrom', 'datenum'), U(idx), 'b', 'LineWidth', 1.5);
xlabel('Time (days)');
ylabel('U (m/s)');
title('南海站点 潮流流速U时间序列 (2022年1月)','fontsize',8);
grid on;
datetick('x', 'dd', 'keeplimits'); % X 轴显示日期
sgtitle('南海站点 潮流流速时间序列 (2022年1月)','fontsize',10)
print(gcf,'-dpng','-r1080','南海站点流速.png')

以上代码及实验过程仅供参考。