matlab最佳高度控制,用MATLAB畫925hPa位勢高度分布圖

用GFSdata數據畫單月平均位勢高度分布圖

function []=draw(y,m,d,file)

%畫單月平均

%GFSdata

%--載入toolbox--%

%setup_nctoolbox

%--計算經緯度對應的idx--%

[lon1,lon2,lat1,lat2]=findpos(60,140,15,50);%%%%%%%%%LINR caxis & axis

figure(1)

%--保存平均高度的矩陣--%

data=zeros(lat2-lat1+1,lon2-lon1+1);%%%%%%%%%%lat,lon

%--用於findmax()--%

%maxx=10000;

%--讀和處理數據文件--%

for i =1:d

i

if i<10

name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_00_000',file,y,m,i);

name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_06_000',file,y,m,i);

name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_12_000',file,y,m,i);

name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s0%d_18_000',file,y,m,i);

else

name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_00_000',file,y,m,i);

name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_06_000',file,y,m,i);

name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_12_000',file,y,m,i);

name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\data\\%s\\gfs05d_%s%s%d_18_000',file,y,m,i);

end

data=cutadd(name0,data,lon1,lon2,lat1,lat2);

data=cutadd(name6,data,lon1,lon2,lat1,lat2);

data=cutadd(name12,data,lon1,lon2,lat1,lat2);

data=cutadd(name18,data,lon1,lon2,lat1,lat2);

%{

maxx=findmax(name0,maxx);

maxx=findmax(name6,maxx);

maxx=findmax(name12,maxx);

maxx=findmax(name18,maxx);

%}

end

data=data./(4*d);

mydraw(data,name0,lon1,lon2,lat1,lat2);

%--設置顯示範圍--%

axis([60,140,15,50])%lon,lat

title(sprintf('20%s-%s',y,m));

xlabel('lon')

ylabel('lat')

%--設定colorbar範圍--%

caxis([600,900]);%650,800

%--保存圖片--%

print(1,'-dbmp',sprintf('C:\\Users\\kong\\KHIMatlab\\925\\picture\\GFS\\pic4\\SMALLuni20%s-%s',y,m));

end

function [data]=cutadd(name,data,lon1,lon2,lat1,lat2)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

%--選取數據範圍--%

data1=iso.data(1,4,lat1:lat2,lon1:lon2);%%%%%%%%%%lat,lon

%--二維化--%

data1=squeeze(data1);

data=data+data1;

da=[];

iso=[];

data1=[];

end

function []=mydraw(data,name,lon1,lon2,lat1,lat2)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

%--讀取經緯度--%

pos1=iso.grid_interop(1,4,:,:);

lat=pos1.lat(lat1:lat2);%%%%%%%%%%

lon=pos1.lon(lon1:lon2);%%%%%%%%%%

%--讀取地圖--%

gx=shaperead('C:\Users\kong\KHIMatlab\ChinaMap\province\bou2_4p.shp','usegeocoords',true);

[gca,h]=contourf(lon,lat,data);

%--畫等值線--%

clabel(gca,h);

%--設定colorbar文字--%

h=colorbar;

set(get(h,'Title'),'string','gpm');

geoshow(gx,'FaceColor','none','EdgeColor','w','LineWidth',1.5);

max(max(data))

min(min(data))

end

function [lon1,lon2,lat2,lat1]=findpos(lon01,lon02,lat01,lat02)

name='C:\Users\kong\KHIMatlab\925\data\1201\gfs05d_120101_18_000';

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

pos=iso.grid_interop(1,4,:,:);

lon1=find(pos.lon==lon01)

lon2=find(pos.lon==lon02)

lat1=find(pos.lat==lat01)

lat2=find(pos.lat==lat02)

end

function [maxx]=findmax(name,maxx)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

data1=iso.data(1,4,:,:);

data1=squeeze(data1);

if min(min(data1))

maxx= max(max(data1))

end

da=[];

iso=[];

data1=[];

end

用JRA25數據畫多年單月平均位勢高度分布圖

function []=drawncep(m,d)

%畫多年單月平均

%JAR

%--載入toolbox--%

%setup_nctoolbox

%--計算經緯度對應的idx--%

[lon1,lon2,lat1,lat2]=findpos(95,120,10,30);%%%%%%%%%LINR caxis & axis

figure(1)

%--保存平均高度的矩陣--%

data=zeros(lat2-lat1+1,lon2-lon1+1);%%%%%%%%%%lat,lon

%--用於findmax()--%

%maxx=10000;

%--讀和處理數據文件--%

for yy=6:11

if yy<10

y=sprintf('0%d',yy)

else

y=sprintf('%d',yy)

end

for i =1:d

i

if i<10

name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d00',y,m,y,m,i);

name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d06',y,m,y,m,i);

name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d12',y,m,y,m,i);

name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s0%d18',y,m,y,m,i);

else

name0=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d00',y,m,y,m,i);

name6=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d06',y,m,y,m,i);

name12=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d12',y,m,y,m,i);

name18=sprintf('C:\\Users\\kong\\KHIMatlab\\925\\JAR\\20%s%s\\anl_p.20%s%s%d18',y,m,y,m,i);

end

data=cutadd(name0,data,lon1,lon2,lat1,lat2);

data=cutadd(name6,data,lon1,lon2,lat1,lat2);

data=cutadd(name12,data,lon1,lon2,lat1,lat2);

data=cutadd(name18,data,lon1,lon2,lat1,lat2);

%{

maxx=findmax(name0,maxx);

maxx=findmax(name6,maxx);

maxx=findmax(name12,maxx);

maxx=findmax(name18,maxx);

%}

end

end

data=data./(4*6*d);

mydraw(data,name0,lon1,lon2,lat1,lat2);

%--設置顯示範圍--%

axis([95,120,10,30])%lon,lat

title(sprintf('%s',m));

xlabel('lon')

ylabel('lat')

%--設定colorbar範圍--%

caxis([650,850]);%650,800

%--保存圖片--%

print(1,'-dbmp',sprintf('C:\\Users\\kong\\KHIMatlab\\925\\picture\\JAR\\pic1\\uni%s',m));

end

function [data]=cutadd(name,data,lon1,lon2,lat1,lat2)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

%--選取數據範圍--%

data1=iso.data(1,22,lat1:lat2,lon1:lon2);%%%%%%%%%%lat,lon

%--二維化--%

data1=squeeze(data1);

data=data+data1;

da=[];

iso=[];

data1=[];

end

function []=mydraw(data,name,lon1,lon2,lat1,lat2)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

%--讀取經緯度--%

pos1=iso.grid_interop(1,22,:,:);

lat=pos1.lat(lat1:lat2);%%%%%%%%%%

lon=pos1.lon(lon1:lon2);%%%%%%%%%%

%--讀取地圖--%

gx=shaperead('C:\Users\kong\KHIMatlab\ChinaMap\province\bou2_4p.shp','usegeocoords',true);

[gca,h]=contourf(lon,lat,data);

%--畫等值線--%

clabel(gca,h);

%--設定colorbar文字--%

h=colorbar;

set(get(h,'Title'),'string','gpm');

%--畫中國地圖--%

geoshow(gx,'FaceColor','none','EdgeColor','w','LineWidth',1.5);

%--畫世界地圖--%

geoshow('landareas.shp', 'FaceColor','none','EdgeColor','w','LineWidth',1);

max(max(data))

min(min(data))

end

function [lon1,lon2,lat2,lat1]=findpos(lon01,lon02,lat01,lat02)

name='C:\Users\kong\KHIMatlab\925\JAR\200605\anl_p.2006050100';

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

pos=iso.grid_interop(1,22,:,:);

lon1=find(pos.lon==lon01)

lon2=find(pos.lon==lon02)

lat1=find(pos.lat==lat01)

lat2=find(pos.lat==lat02)

end

function [maxx]=findmax(name,maxx)

da=ncgeodataset(name);

iso=da.geovariable('Geopotential_height_isobaric');

data1=iso.data(1,22,:,:);

data1=squeeze(data1);

if min(min(data1))

maxx= max(max(data1))

end

da=[];

iso=[];

data1=[];

end