matlab读取NCEP细网格资料及计算散度、涡度平流、温度平流、K指数等要素

F117_ren_0

首先感谢气象家园各位元老，学习和借鉴了很多有用的东西，具体名字就不再点出来了，再次表示感谢！读取FNL数据首先安装read_grib工具箱和M_map工具箱，具体安装方法请百度！
1,读取FNL数据
grib_struct=read_grib('F:\FNL\fnl_19990801_00_00_c',-1); %%%读取NCEP再分析资料
    read_grib('F:\FNL\fnl_19990801_00_00_c','inv'); %生成变量目录表文档
    grib_struct(13).gds;
    Hight_500hPa=grib_struct(213).fltarray;
%Hight_1000_500hPa=Hight_1000hPa_01;
for i=1:91;
    for j=1:181;
        lon(i,j)=181-j;
        lat(i,j)=91-i;
        Hight(i,j)=Hight_500hPa((i-1)*360+j)/10.0;
    end
end

这里是grib1格式，好像是07年以后是grib2格式，但读法一样，只是从排列发生变化，具体可以见%生成变量目录表文档%想读取的数据 Hight_500hPa=grib_struct(213).fltarray里，只要改动括号里的数字就ok了

2,计算散度
clc
clear all
%%%%%%%%%%%%%%%%%%%%读取NFL
    grib_struct=read_grib('F:\FNL\fnl_20100311_00_00_c',-1); %%%读取NCEP再分析资料
    read_grib('F:\FNL\fnl_20100311_00_00_c','inv'); %生成变量目录表文档
    %grib_struct(13).gds;  
    U850=grib_struct(157).fltarray;
    V850=grib_struct(158).fltarray;
for i=34:50;
    for j=70:100;
        lon(i,j)=100-j;
        lat(i,j)=50-i;
        u850(i,j)=U850((i-1)*360+j);
        v850(i,j)=V850((i-1)*360+j);
    end
end
u850=u850(34:end,70:end);
v850=v850(34:end,70:end);
%计算散度
   R0=6371004;  pi=3.1415926;
   lat0=50;  lon0=70;
    for lati=50:-1:34
      latm(50-lati+1)= distance(lat0,lon0,lati,lon0)/180*pi*R0;
    end

   for loni=70:1:100
      lonm(loni-69)= distance(lat0,lon0,lat0,loni)/180*pi*R0;
   end

   [xx,yy]=meshgrid(lonm,latm);
   div=divergence(xx,yy,u850,v850);  
  save D:\matt\div850 div xx yy
850hPa散度计算完成，这里的经纬度我选的是新疆的范围，如有需要可自行改动。

                               
登录/注册后可看大图

3、计算涡度平流和温度平流
首先读取出UV，涡度，温度
lon=70:100;
lat=34:50;
PI=3.1415926;
R0=6371004;
[dtx850,dty850]=gradient(w850,1,1);

dx=gradient(lon,1)*PI/180;
dy=gradient(lat,1)*PI/180;
wp850=-((u850.*dtx850)./(dx.'*cos(lat*PI/180))'+v850.*(dty850'*diag(1./dy))')/R0;

以上w850为850hPa的涡度，wp850为计算出来的850hPa的涡度平流
[dtx850,dty850]=gradient(t850,1,1);
tp850=-((u850.*dtx850)./(dx.'*cos(lat.*PI/180))'+v850.*(dty850'*diag(1./dy))')/R0;
以上t850为850hPa的温度，tp850为计算出来的850hPa的温度平流。

                               
登录/注册后可看大图

                               
登录/注册后可看大图

4、计算K指数
clc
clear all
%%%%%%%%%%%%%%%%%%K指数
% K=(T850-T500)+Td850-(T700-Td700)
%%%%%%%%%%%%%%%%%%
load D:\matt\Temp2013041712
load D:\matt\RH2013041712
E0=611.2;
%%%%t>0  a=7.5; b=237.3;
%%%%t<=0 a=9.5 b=265.5
for i=1:17;
    for j=1:31;
if t850(i,j)>0;
    a=7.5;
    b=237.3;
else
    a=9.5;
    b=265.5;
end
qq1(i,j)=a*t850(i,j);
qq2(i,j)=b+t850(i,j);
    end
end
qq=qq1./qq2;
Es=E0*(10.^qq);%%%%%%%%饱和水汽压

E=Es.*rh850.*0.01;  %%%该相对湿度下水汽压

T1=b*log10(E./Es);
T2=a-log10(E./Es);

Td850=T1./T2;
save D:\matt\Td850 Td850
clear all
load D:\matt\Temp2013041712
load D:\matt\RH2013041712
E0=611.2;
%%%%t>0  a=7.5; b=237.3;
%%%%t<=0 a=9.5 b=265.5
for i=1:17;
    for j=1:31;
if t700(i,j)>0;
    a=7.5;
    b=237.3;
else
    a=9.5;
    b=265.5;
end
qq1(i,j)=a*t700(i,j);
qq2(i,j)=b+t700(i,j);
    end
end
qq=qq1./qq2;
Es=E0*(10.^qq);%%%%%%%%饱和水汽压

E=Es.*rh700.*0.01;  %%%该相对湿度下水汽压

T1=b*log10(E./Es);
T2=a-log10(E./Es);

Td700=T1./T2;
save D:\matt\Td700 Td700
clear all
clc
load D:\matt\Temp2013041712
load D:\matt\Td850
load D:\matt\Td700
K=(t850-t500)+Td850-(t700-Td700);
save D:\matt\Kzhishu2013041712 K

首先要读出湿度和温度，我以上程序中load出来的就是存好的数据，后面就是计算了。

                               
登录/注册后可看大图

5、画图
load D:\matt\wdpll2010031312;
fnshp_P='D:\Program Files\toolbox\m_map\bou2_4m\bou2_4p.shp';% 省界
S = shaperead(fnshp_P, 'UseGeoCoords', true,...  
  'Selector',{@(NAME) strcmp(NAME,'新疆维吾尔自治区'), 'NAME'});

figure(1);
set (gca,'Position',[0.1,0.1,0.85,0.85], 'color','w')

lo=70:100;
la=34:50;
[lo,la]=meshgrid(lo,la);
wp850=wp850.*10.^10;
wp850(isnan(wp850)==1)= 0;
[x,y]=meshgrid(70:0.1:100,34:0.1:50);
KK=griddata(lo(:),la(:),wp850(:),x,y,'v4');
m_proj('Equidistant Cylindrical','lon',[70,100],'lat',[34,50]);
[cs,h]=m_contour(x,y,KK,'LineWidth',2);
set(h,'ShowText','on','LevelList',[-500 -350 -200 -100 -50 0 50 100 200 350 500]);
clabel(cs,h,'fontsize',5);
hold on;
shading flat;
m_grid('box','fancy');
%axis equal
hold on
m_plot(S.Lon,S.Lat,'-k','linewidth',2);
hold on
title('850涡度平流');
close(1)

其他的也一样，就不一一表示了。

Gongql

赞，好贴！
计算温度/涡度平流时，我觉得楼主的公式存在计算符号的问题，矩阵与矩阵相乘为*，要求前面的列数与后面矩阵的行数相等，而.*要求两个矩阵完全一样，所以我觉得应该为
tp850=-((u850.*dtx850)./(dx*cos(lat'*PI/180))+v850.*(dty850*diag(1./dy)))/R0;
希望楼主看到可以批评指正！

F117_ren_0

图片不让我上传，没图没真相啊，呵呵

F117_ren_0

print出来的图片不让贴出来？奇怪!

传说中的谁

楼主好牛逼的样子，竟然知道grib数据的内部结构，真了不起

Aires

为啥不能加图呢？

F117_ren_0

传说中的谁 发表于 2014-8-4 17:21
楼主好牛逼的样子，竟然知道grib数据的内部结构，真了不起

还是笨笨牛b，哈哈

jiajian0512

貌似楼主很牛逼的样子，狂支持，下载了！

cherven

这么详细啊 学学

F117_ren_0

Aires 发表于 2014-8-4 17:24
为啥不能加图呢？

不知道，也许新疆这边不让上图吧，我已经让广州笨笨帮我上图了

传说中的谁

高了半天终于插进来了，好费劲