简介
风场数据来源ERA5再分析数据集，下载方式网上有许多教程就不重复造轮子了 ：
Copernicus Climate Data Store | Copernicus Climate Data Store，也可以直接下载我后文分享的已经下载好的数据进行学习。
台风案例选择了闻名的2015年9号超强台风灿鸿，下图绘出了台风灿鸿的轨迹，数据来源海洋科学大数据中心的中国近海台风路径聚集数据集
风场数据时间选择2015-07-01到2015-07-14，时间分辨率1h，经度范围（110，160）°E，纬度范围（5，45）°N
绘出如下台风期间风场演变，可以清晰看到台风的结构和影响范围
固然也可以生成动画版本的，方便用于演示
原始数据
下面直接分享上面提到的台风轨迹、台风期间再分析数据，仅作学习用途，具体请参考相关网站
通过百度网盘分享的文件：灿鸿风场
链接：https://pan.baidu.com/s/1uJgGYYhqdQa0OKkIETo6nw?pwd=6666
提取码：6666
代码
接下来直接分析绘制风场及其动画的代码
起首读取路径下的数据
- % 读取ERA5数据
- % ncinfo("ERA5.nc") 用于了解NC文件结构
- lon = ncread("ERA5.nc",'longitude'); %longitude
- lat = ncread("ERA5.nc",'latitude'); %latitude
- date = ncread("ERA5.nc",'time'); %时间
- date = datetime(1900,1,1) + double(date(:))/24; %时间格式转化
- u = ncread("ERA5.nc",'u10'); %速度u分量
- v = ncread("ERA5.nc",'v10'); %速度v分量
- speed = sqrt(u.^2+v.^2); %速度
- % 生成m_map project
- m_proj('Equidistant Cylindrical','long',[110 160],'lat',[5 45]);
- % 计算坐标对应地图位置
- X=[];
- Y=[];
- for i=1:length(lon)
- for j=1:length(lat)
- [CX,CY]=m_ll2xy(lon(i),lat(j));
- X(i,j)=CX;
- Y(i,j)=CY;
- end
- end
- for i=1:6:length(date)
- figure;
- % 绘制风场
- hold on
- contourf(X,Y,(speed(:,:,i)),15,'linestyle','none');
- % 海岸线
- m_coast('patch',[.8 .8 .8]);
- % 绘制边框
- m_grid('box','on','tickdir','in','backcolor',[1,1,1],'xtick',[],'ytick',[],'LineWidth',1.5);
- hold off
- % 生成colorbar
- cb = colorbar;
- cb.FontSize=14;
- cb.Label.String='Wind speed (m/s)';
- % 颜色设置
- cr1 = flipud(slanCM(167));%彩色
- cr2 = slanCM(25);%灰色
- colormap([cr2(1:3:83,:);flipud(cr1(1:2:end-3,:))]);% 将低风速设置为灰色，突出高风速台风区域
- % 日期
- XLB = xlabel(string(date(i)),'FontSize',15,'FontWeight','bold');
- set(XLB,'Units','normalized','Position',[0.5,0],'HorizontalAlignment','center','VerticalAlignment','top')
- frame = getframe(gcf);
- im = frame2im(frame);
- [I,map] = rgb2ind(im,256);
- if i == 1
- imwrite(I,map,'windfield.gif','gif', 'Loopcount',inf,'DelayTime',1);
- else
- imwrite(I,map,'windfield.gif','gif','WriteMode','append','DelayTime',1);
- end
- close
- end
- winopen('windfield.gif')% 预览结果
最后是台风风场瞬时图
- %% 提取1:12日每天0点的数据
- index = hour(date) == 0& day(date) ~= 13;
- date = date(index);
- u = u(:,:,index);
- v = v(:,:,index);
- speed = speed(:,:,index);
- %% 绘制台风风场图
- % 分子图每幅图一个时刻
- t = tiledlayout(3,4);
- % 生成m_map project
- m_proj('Equidistant Cylindrical','long',[110 160],'lat',[5 45]);
- for i = 1:12
- % 生成子图
- nexttile(i);
- hold on
- % 绘制等高线图
- contourf(X,Y,speed(:,:,i),15,'LineStyle','none')
- caxis([0,32])% 颜色上下限
- % 海岸线图
- m_coast('patch',[.8 .8 .8]);
- % m_gshhs_i('patch',[.8 .8 .8]);
- % 绘制边框
- m_grid('box','on','tickdir','in','backcolor',[1,1,1],'xtick',[],'ytick',[],'LineWidth',1.5);
- hold off
- % xlabel
- XLB=xlabel(string(date(i)),'FontSize',15,'FontWeight','bold');
- set(XLB,'Units','normalized','Position',[0.5,0],'HorizontalAlignment','center','VerticalAlignment','top')
- end
- % 设置间隙
- t.Padding='none';
- t.TileSpacing="none";
- % 生成colorbar,实用Tiledlayout的设置
- cb = colorbar;
- cb.FontSize=14;
- cb.Layout.Tile='east';
- cb.Label.String='Wind speed (m/s)';
- % 颜色设置
- cr1 = flipud(slanCM(167));%彩色
- cr2 = slanCM(25);%灰色
- colormap([cr2(1:3:83,:);flipud(cr1(1:2:end-3,:))]);% 将低风速设置为灰色，突出高风速台风区域
- % 输出图片设置
- figWidth = 30; % 设置图片宽度
- figHeight = 20; % 设置图片高度
- set(gcf,'PaperUnits','centimeters'); % 图片尺寸所用单位
- set(gcf,'PaperPosition',[0 0 figWidth figHeight]);
- fileout = 'Windfield_Chan_hom'; % 输出图片的文件名
- print(gcf,[fileout,'.tif'],'-r300','-dtiff'); % 设置图片格式、分辨率
- winopen([fileout,'.tif'])
结尾
以上就是本次分享的全部内容了
如果觉得有用的话，
感谢您的点赞收藏和关注~
鼓励我继承更新吧！
接待评论区讨论
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。 |
