matlab画大气环流

%% 本脚本利用wgrib2先把grib2文件转换为nc文件；在进行读取；画出位势高度场
% 步骤如下:
% 1 下载wgrib2并“安装”（为什么这里的“安装”打了引号呢？因为严格意义上来wgrib2并不需要安装，只需要在环境变量中进行设置即可）
% （1）下载链接：https://pan.baidu.com/s/1EFOBbGvHMPFg-iZS25zAwg 提取码：8bwn
% （2）wgrib2设置windows下wgrib/wgrib2和python读取grib2数据_meryl2001的博客-CSDN博客。
% 2 在cmd中用wgrib2软件将grib2文件转成nc文件
% （1）同时按下windows和R键，输入cmd，回车进入这个黑黑的界面。
% (2）转换命令：wgrib2 grib2file -netcdf ncfile。（grib2file为grib2文件的名字，ncfile为转换出的nc文件的名字）

%%
%% 公众号：海洋与大气科学
%% 作者mapm
clear;clc;close all;% ty_name 指的是台风号；year指的是台风发生的年份；
% month指的是台风发生的月份；%area台风所在范围；% inter1、2 是画图xlabel间隔；
ty_name='1515';year='2015';month='08';area=[90 160;0 40];inter1=10;inter2=5;
days=[16;27];
%% 读取台风真实路径；
%% 本脚本读取台风路径：2018年24号”潭美“台风：
%% read txt
filename=['20151515.txt'];
fid=fopen(filename,'r');%open file
eng1='000000000';
row=0;%row at start is zero.
while ~feof(fid) %if eixt file ,then run
    row=row+1; % read file according row
    a=fgetl(fid);% data give a;
    num=length(a);
    eng1(row,1:num)=a;
end
fclose(fid);%close file ;
clear row a ans fid filename num;% 选取自己需要的；
%%
end1=eng1(2:end,:);
clear eng1;
time=str2num(end1(:,7:10));% read typhoon's time ;
tylon=str2num(end1(:,18:21))./10;% read typhoon's lon ;
tylat=str2num(end1(:,14:16))./10;% read typhoon's lat;
time1=str2num(end1(:,9:10));% read typhoon's time's hour
time2=str2num(end1(:,7:8));% read typhoon's time's day;
%% 对时间日期处理一下，补齐0；例如：1号补齐为01；小时6补齐为06；并且
%构造新的time_new,为完整日期，例如0706
for time1_i=1:length(time1)
    time111=time1(time1_i);
    time222=time2(time1_i);
    if length(num2str(time111))==1
        time1_new(time1_i,:)=['0',num2str(time111)];
    else
        time1_new(time1_i,:)=num2str(time111);
    end
    if length(num2str(time222))==1
        time2_new(time1_i,:)=['0',num2str(time222)];
    else
        time2_new(time1_i,:)=num2str(time222);
    end
end
time_new=[time2_new time1_new];
clear typhoon_name end1  time time1 time111 time1_i time1_new time2 time222 time2_new
%% 下面是fnl
filename='fnl_20150816_00_00.nc';
%% 读取风速和位势高度；500hpa;
lon = double(ncread(filename,'longitude'));
lat = double(ncread(filename,'latitude'));
GHT_500=double(ncread(filename,'HGT_500mb'));
U_500=double(ncread(filename,'UGRD_500mb'));
V_500=double(ncread(filename,'VGRD_500mb'));
%% 区域；
% 2015年15号台风天鹅；
% 0-40N；100-160；毕鑫鑫等超强台风天鹅路径突变过程机理研究
%         area_select=[100 160;0 40];
ln=find(lon>=area(1,1)&lon<=area(1,2));
la=find(lat>=area(2,1)&lat<=area(2,2));
lonx=lon(ln);laty=lat(la);
[x,y]=meshgrid(lonx,laty);
x=x';y=y';
%% 把区域变量找出来；
HGT=GHT_500(ln,la);
U=U_500(ln,la);
V=V_500(ln,la);
%% 画图；
%         min(min(HGT))
% 大于5880的数据；
[xj,yj]=size(HGT);
for xi=1:xj
    for yi=1:yj
        if HGT(xi,yi)>=5880
            HGT_5880(xi,yi)=HGT(xi,yi);
        else
            HGT_5880(xi,yi)=nan;
        end
    end
end
%end
%%
% 读取颜色；
color=colormore_46(100);
num=46;
fig_h=figure;
set(gcf,'position',[50,50,1000,800],'color','w')
m_proj('miller','lon',[area(1,1) area(1,2)],'lat',[area(2,1) area(2,2)]);
m_contourf(x,y,HGT_5880,50,'linestyle','none')% 大于5880的数据；
caxis([5.88e3 5.9e3])
hold on
[c,h]=m_contour(x,y,HGT,[5700:20:5920],'linecolor','k','ShowText','on','LineWidth',1.5);
clabel(c,h,'fontsize',14,'fontweight','bold','LabelSpacing',1000);
hold on
a=3;inter=2;
m_quiver(x(1:inter:end,1:inter:end),y(1:inter:end,1:inter:end),U(1:inter:end,1:inter:end),V(1:inter:end,1:inter:end),a,'AutoScale','off','color','k');
% 矩形并且白色背景；
hold on
colormap(color)%gaiyanse 改颜色；运行节就行了我数据泡好了
colorbar
ty_lon=tylon(9);ty_lat=tylat(9);
m_plot(ty_lon,ty_lat,'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k')%LineWidth——指定线宽
m_coast('linewidth',1.5,'color',[0.6,0.6,0.6]);
m_grid('ytick',[area(2,1):inter2:area(2,2)],'fontsize',16,'fontweight','bold','linewidth',3);
set(gca,'xtick',[area(1,1):inter1:area(1,2)],'xticklabel',[area(1,1):inter1:area(1,2)],'linewidth',2);
export_fig(fig_h,['colormore_',num2str(num),'台风500hpa位势高度和风场图'],'-dpng','-r300');