利用GFS多天的数据，做平均的边界层高度区域分布图，并对区域外进行白化

祈祷

常得益于气象家园，感谢大侠们的分享。

如下是综合论坛上学到的程序和自己的经验编写的小程序，供大家参考。

clear all   

c0=filesep;

x1=112;x2=115;y1=22;y2=25;

directory = strcat('E:',c0,’dataset’,c0,’GFS’);

file_name = importdata(strcat(directory,c0,'filename.txt'));

for i_name = 1:length(file_name)

Data_File = strcat(directory,c0,file_name{i_name,1})  

fid = open(Data_File);

if fid > 0

        grib = ncgeodataset(Data_File);

        utc_time=double(grib.data('time'));

        Data.lat = double(grib.data('lat'));

        Data.lon = double(grib.data('lon'));

        Data.PBLH = squeeze(double(grib.data('Planetary_Boundary_Layer_Height_surface')));

        clear ii1 jj1; [ii1,jj1] = find(Data.lat>y2&Data.lat<y1);        

        clear ii2 jj2; [ii2,jj2] = find(Data.lon>x1&Data.lon<x2);

        PBLH(:,:,i_name) = Data.PBLH(ii1,ii2);

    end

    fclose(fid);

end

PBLH_mean = nanmean(PBLH,3);

a1= worldmap([y1 y2],[x1 x2])

[lon lat]=meshgrid(data.lon(ii2),data.lat(ii1));

contourfm(lat,lon,PBLH_mean,50,'Linestyle','none')

shading interp

colorbar

caxis([-300,400])

l=colorbar;

set(l,'ytick',[-300:100:400]);

sf2 = shaperead(strcat(directory,c0,'Guangzhou.shp'),'usegeocoords',true,'attributes',{'Lat','Lon'});

clear aa bb;

aa = sf2.Lat;

bb = sf2.Lon;

aa = [aa 22 22 25 22 22];

bb = [bb 113.247075 112 112 115 115 113.247075];

h = patchm(aa,bb,[1 1 1]);

set(h,'EdgeColor','w')

sf1 = shaperead(strcat(directory,c0,'Guangzhou.shp'),'usegeocoords',true,'attributes',{'Lat','Lon'});

geoshow([sf1.Lat],[sf1.Lon],'Color', 'black')              setm(gca,'PLabelLocation',1,'grid','on','Gcolor','k','PLineLocation',1,'mLabelLocation',1,'grid','on','Gcolor','k','mLineLocation',1,'fontsize',fontsize0+2)

1.     读取GFS资料

（1）     filename.txt里面包括了该文件夹下你想要处理的文件名称。例如

gfs2015070100.grb2

gfs2015070112.grb2

gfs2015070200.grb2

gfs2015070212.grb2

gfs2015070300.grb2

gfs2015070312.grb2

（2）     filename.txt文件的生成

点开始—运行—键入cmd回车（或者win+r；cmd），窗口输入：dir E:\dataset\GFS /b>E:\dataset\GFS\filename.txt.

（3）     说明，要安装好nctoolbox才能用matlab读取grib2码格式的数据。具体方法论坛里面有。

2.     Patch方法

原本patch会将区域内白化，所以要在原来的数据中加入边界的几个点，使得地图数据和边界点组成新的区域，白化该区域。

论坛里已经有一些白化（掩膜）的方法，本方法提供多一种选择。

Lighting

其实楼主可以先给出效果图，程序作为附件单独收取贡献。

祈祷

Lighting 发表于 2016-10-26 11:35
其实楼主可以先给出效果图，程序作为附件单独收取贡献。

谢谢建议！发的少，没经验。
先这样吧，懒得改了，需要的同学就借鉴一下。收取贡献就是凑个热闹。

瓦伦宇

好东西，谢谢楼主。但是不知道楼主是否遇到过这样的情况。利用ncgeodataset函数读取grb文件，比如nc = ncgeodataset('D:\2008ECMWF\ecmf_20080102.grb');读完之后会生成两个同名文件，ecmf_20080102.grb.gbx9和ecmf_20080102.grb.ncx,直接删除的话说是文件正在被matlab使用，只能退出matlab才能删除掉，请问楼主这个问题怎么解决，能否不要生成这两个文件，或者是可以直接删除掉？

祈祷

瓦伦宇 发表于 2016-12-28 19:34
好东西，谢谢楼主。但是不知道楼主是否遇到过这样的情况。利用ncgeodataset函数读取grb文件，比如nc = ncge ...

我没遇到这样的情况啊。