MATLAB绘制台风风速剖面

MATLAB编码
提取利奇马台风持续期间每一时刻最小风速的经纬度，作为台风位置，将各个位置作为点标出，并连接起来。数据采用era5的nc文件，文件名为“era5_liqima.nc”,绘制利奇马台风风速剖面图，按每一个台风所处位置的前后最近位置与本位置互相之间的连线和正北方向的夹角的平均值，把这个平均值夹角画在此位置上，然后作垂线，以这个垂线为切面，绘制出以经度为x轴 最大风速为y轴的风速剖面图。在台风持续时间之内，每隔一小时绘制一张图。

需要MATLAB netCDF 支持。

简单写一段，假设你的nc文件有'longitude'、'latitude'、'time'、'wind_speed'，到时候你用的时候根据实际你的字段名去替换

% 加载nc文件
filename = 'era5_liqima.nc';
lon = ncread(filename,'longitude');
lat = ncread(filename,'latitude');
time = ncread(filename,'time');
wind_speed = ncread(filename,'wind_speed');

% 预处理数据
[~, num_time] = size(time);
min_speed_index = zeros(num_time, 2); % 存储每个时刻最小风速的经纬度索引

% 查找每一时刻最小风速的经纬度
for i = 1:num_time
    [~, min_index] = min(wind_speed(:,:,i), [], 'all', 'linear');
    [min_lat_index, min_lon_index] = ind2sub(size(wind_speed(:,:,i)), min_index);
    min_speed_index(i, :) = [min_lat_index, min_lon_index];
end

% 绘制台风路径
figure;
plot(lon(min_speed_index(:, 2)), lat(min_speed_index(:, 1)), '-o');
xlabel('Longitude');
ylabel('Latitude');
title('Typhoon Track');

% 按每一时刻，绘制风速剖面图
for i = 1:num_time
    figure;
    plot(lon, wind_speed(min_speed_index(i, 1), :, i));
    xlabel('Longitude');
    ylabel('Wind speed');
    title(['Wind speed profile at time ', num2str(time(i))]);
    pause(1); % 暂停一秒，可根据需要调整
end

% 导入nc文件
filename = 'era5_liqima.nc';
ncdisp(filename);  % 显示nc文件信息

% 读取经纬度和风速数据
lon = ncread(filename, 'longitude');
lat = ncread(filename, 'latitude');
time = ncread(filename, 'time');
wind_speed = ncread(filename, 'wind_speed');

% 初始化台风位置
typhoon_pos = zeros(length(time), 2);

% 计算每一时刻最小风速的经纬度
for i = 1:length(time)
    [min_wind, idx] = min(wind_speed(:,:,i), [], 'all', 'linear');
    [row, col] = ind2sub(size(wind_speed(:,:,i)), idx);
    typhoon_pos(i, :) = [lon(col), lat(row)];
end

% 绘制台风路径
figure;
plot(typhoon_pos(:,1), typhoon_pos(:,2), '-o');
xlabel('Longitude');
ylabel('Latitude');
title('Typhoon Path');

% 计算夹角并绘制风速剖面图
for i = 2:length(time)-1
    % 计算夹角
    angle = atan2(typhoon_pos(i+1,2)-typhoon_pos(i-1,2), typhoon_pos(i+1,1)-typhoon_pos(i-1,1));
    
    % 创建垂线
    perp_angle = angle + pi/2;
    x = typhoon_pos(i,1) + cos(perp_angle)*[-1, 1]*100;
    y = typhoon_pos(i,2) + sin(perp_angle)*[-1, 1]*100;
    
    % 绘制风速剖面图
    figure;
    plot(x, wind_speed(:,i));
    xlabel('Longitude');
    ylabel('Wind Speed');
    title(['Wind Speed Profile at Time ', num2str(time(i))]);
end