1.M_MAPを使ってベクトルなどを描画(例:冬季の海氷密接度、海氷ベクトル、海面気圧、海上風の平均場)

※それぞれのバイナリは予め用意しておく。
M_MAPパッケージをダウンロードしてパスを通しておくこと。\MATLABのパス\toolbox\local\setup.mでaddpathしておく。
 

clear ; clc ; clf ; close all
%% make freezing season mean fied ice, Uice, SLP and UVwnd
%------- read latlon info. --------------- AMSRデータの緯度経度情報を読み込み
fid = fopen('/home/kmizobat/work/AMSR/AMSR-PS/IC0/latlon_low_NP','r');
latlon = fread(fid,[900,1800],'short');
lat = latlon(:,1:900)/100. ; lon = latlon(:,901:1800)/100.;

%% ------- read IBCAO flat file 10801x751 2-byte integer [65-90N, 0-360E]
ibf = ['/home/kmizobat/work/MAP/IBCAO/IBCAO_ver2_23_GEO_netCDF_2min.flat']; 予め用意したIBCAOのデータを読み込み
fidib=fopen(ibf,'r') ; ibd=fread(fidib,[10801,751],'short'); fclose(fidib);
ibd2=-ibd(5401:8101,:); % Extract
sprintf('%s','finish reading IBCAO')

iblon=zeros(2701,751) ; iblat=iblon;
for j=1:751
iblat(:,j)=90-(j-1)*((90-65)/751);
for i=1:2701
iblon(i,:)=-180+(i-1)*(360/10801);
end
end
for jk=1:751
for ik=1:2701
if(iblon(ik,jk) < 180 & iblon(ik,jk) >= 0)
ibd2(ik,jk)=NaN;
end
if(iblat(ik,jk) > 81 | iblat(ik,jk) < 65)
ibd2(ik,jk)=NaN;
end
end
end
sprintf('%s','finish modifying IBCAO')



%% loop start ここからForループの始まり
for yr=2002:2009
yr2=num2str(yr);
icef=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/' yr2 '.fallmean.flat']%海氷密接度ファイル指定
uicef=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/' yr2 '.fallmean.mve.mat']%氷速ベクトル東西成分のファイル指定
vicef=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/' yr2 '.fallmean.mvn.mat']%氷速ベクトル南北成分のファイル指定
uwndf=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/uwnd.' yr2 '.fallmean.flat']%風速ベクトル東西成分のファイル指定
vwndf=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/vwnd.' yr2 '.fallmean.flat']%風速ベクトル南北成分のファイル指定
slpf=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/slp.' yr2 '.fallmean.flat']%海面気圧のファイル指定
ofl=['/home/kmizobat/work/AMSR-kamo/winter/sim_raw/mat/FALLMEAN/',yr2,'.frzup.jpg']; ofl%出力ファイル名の指定

fid = fopen(icef,'r'); ice = fread(fid,[900,900],'float'); clear fid;
% ice =imrotate(ice,180);
load(uicef) ; load(vicef);
fid = fopen(uwndf,'r'); uwnd = fread(fid,[192,94],'float'); clear fid;
fid = fopen(vwndf,'r'); vwnd = fread(fid,[192,94],'float'); clear fid;
fid = fopen(slpf,'r'); slp = fread(fid,[144,72],'float'); clear fid;

%% modifying ice conc. data % 欠損データの除去
for jk=1:900 ; for ik=1:900
if(lon(ik,jk) < 180 & lon(ik,jk) > 0)
ice(ik,jk)=NaN;
end
if(lat(ik,jk) > 81)
ice(ik,jk)=NaN;
end

% if(ice(ik,jk) > 90)
% ice(ik,jk)=90;
% end
end ;end

%% modifying icevec. data
for jk=1:151 ; for ik=1:159
if(mve(ik,jk) < -8000 )
mve(ik,jk)=NaN;
end
if(mvn(ik,jk) < -8000 )
mvn(ik,jk)=NaN;
end

end ;end
%% modifying slp data
plon=zeros(144,72) ; plat=zeros(144,72);
for jj=1:72 ; for ii=1:144
plon(ii,jj)=0+(ii-1)*2.5;
plat(ii,jj)=-90+(jj-1)*2.5;
if(plon(ii,jj) < 180 & plon(ii,jj) > 0)
slp(ii,jj)=NaN;
end
if(plon(ii,jj) > -120 & plon(ii,jj) <= 0)
slp(ii,jj)=NaN;
end
end ;end
ind=find(plon >= 180) ; plon(ind)=plon(ind)-360; clear ind;

%% modifying uvwnd data
hlon=zeros(192,94) ; hlat=zeros(192,94);
for jj=1:94
for ii=1:192
hlon(ii,jj)=0+(ii-1)*1.875;
hlat(ii,jj)=-88.542+(jj-1)*1.914893617;
if(hlon(ii,jj) < 180 & hlon(ii,jj) > 0)
uwnd(ii,jj)=NaN;vwnd(ii,jj)=NaN;
end

end
end
ind=find(hlon >= 180) ; hlon(ind)=hlon(ind)-360; clear ind;

%% maping ice conc etc...
m_proj('lambert','lat',[65 81],'lon',[-180 -120]); %ランベルト図法に指定
% m_proj('stereographic','latitude',90,'radius',25, 'rotangle', 180);%ポーラーステレオに指定
hold on;

[cs,hs]=m_contourf(lon,lat,ice,[10:1:100]);
a=flipud(cptcmap('GMT_white-blue', 'mapping', 'showall','ncol',18));
colormap(a) ; h=colorbar; caxis([10 100]);
set(get(h,'ylabel'),'String','Sea ice concentration(%)','fontsize',12);
set(hs,'LineStyle','none')

%% landmask and coast line
% m_grid('fontsize',14,'xtick',10,'ytick',5,'tickdir','in',...
% 'linest','-','linewidth',3,'clip','point'); % stereo
%m_gshhs_h('save','cst');  % 描画領域の海岸線ファイルをセーブ ※ 毎回やる必要なし
m_usercoast('cst','patch',[.4 .4 .4]);% 海岸線の描画
% m_coast('patch',[.3 .3 .3],'edgecolor','k');

%% Contour --- IBCAO
[c2,h2]=m_contour(iblon,iblat,ibd2,[600:600:1800],'k','linewidth',1,'linest','-'); %海底地形のコンターラインを描画
% clabel(c2,h2,'fontsize',12); h22 = clabel(c2,h2);; %コンターラインのプロパティ
% set(h1,'FontName','Arial');
% set(h1,'Fontsize',12);
%[c2,h2]=m_contour(iblon,iblat,ibd2,[40:10:100],'k','linewidth',1,'linest','-');

%% Contour --- SLP
[c,h]=m_contour(plon,plat,slp,[1008:2:1030],'k-','linewidth',2); %海面気圧1008hPa以上を実線のコンターラインで描画
clabel(c,h,'fontsize',12); h1 = clabel(c,h);
set(h1,'FontName','Arial');
set(h1,'Fontsize',12);
[c,h]=m_contour(plon,plat,slp,[980:2:1006],'k-','linewidth',2,'linest','-.'); %海面気圧1006hPa以下を点線のコンターラインで描画
clabel(c,h,'fontsize',12); h1 = clabel(c,h);
set(h1,'FontName','Arial');
set(h1,'Fontsize',12);

%% Vector mapping
s=60 ;% ベクトルのスケールファクター
ik=2 ; jk=2; ; %ベクトルデータの間引き加減を決定
lon2=slon(1:ik:159,1:jk:151);
lat2=slat(1:ik:159,1:jk:151);
ve2=mve(1:ik:159,1:jk:151);
vn2=mvn(1:ik:159,1:jk:151);
% m_vec(s, -146, 68, 20, 0, 'g', 'key', '20 cm s^{-1}','EdgeColor','g');); % ベクトルのリファレンスを描画
% m_vec(s,lon2,lat2,ve2,vn2,'g','shaftwidth',0.7,'headangle',40); % ベクトル描画
hh=m_quiver(lon2,lat2,ve2,vn2,4,'linewidth',1.)); % 氷速ベクトル描画
set(hh,'Color',[0 0.6 1]);
hh=m_quiver(-150, 69, 0.4, 0, 4,'c','linewidth',1.);
set(hh,'Color',[0 0.6 1]);
htv5 = text(0.0, -0.075, 'Uice 20 cm s^{-1}');
set(htv5,'FontSize',8);
set(htv5,'Color',[0 0.6 1]);

s2=30;
ik=2 ; jk=2;
hlon3=imresize(hlon,[192,188],'bilinear');hlat3=imresize(hlat,[192,188],'bilinear'); % imresizeで行列サイズの変更
uwnd3=imresize(uwnd,[192,188],'bilinear');vwnd3=imresize(vwnd,[192,188],'bilinear');
iii=find(hlon3 >= -156 & hlon3 <= -124 & hlat3 >= 65 & hlat3 <= 69);
uwnd3(iii)=NaN;vwnd3(iii)=NaN;

hlon4=hlon3(1:ik:192,1:jk:188);
hlat4=hlat3(1:ik:192,1:jk:188);
uwnd4=uwnd3(1:ik:192,1:jk:188);
vwnd4=vwnd3(1:ik:192,1:jk:188);
m_vec(s2, -150, 67.7, 10, 0, 'r', 'key', 'Wind 10 m s^{-1}','EdgeColor','k'); % 風速ベクトルのリファレンスを描画
m_vec(s2,hlon4,hlat4,uwnd4,vwnd4,'r','shaftwidth',1.2,'headangle',45); % 風速ベクトル描画
% m_quiver(hlon4,hlat4,uwnd4,vwnd4); % 繝吶け繝医Ν謠冗判
% m_quiver(-150, 67.5, 0.4, 0, 4, 'b','linewidth',1.);

m_grid('box','fancy','tickdir','in');
xlabel('Ice(shaded)+Uice(cyan)+SLP(blk)+WND(red)','fontsize',14);
title([yr2 ' Frz-up season'],'fontsize',30,'fontweight','bold');


hold off;
print('-djpeg','-r200',ofl);
close all;
end