用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
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
maxx= max(max(data1))
end
da=[];
iso=[];
data1=[];
end