【Matlab】GIM TEC的變化圖 - colormap, colorbar, caxis 的應用

GIMTEC濃度變化

tags: Matlab 高層大氣科學

目標

1. 下載2014年10月19日之後任一天GIM資料
   • 注意檔名結尾，應為 I.Z
   • 檔名範例：CODG0010.15I.Z
2. 畫出每小時GIM TEC的變化圖，共25張)
3. 每張圖片上標示日期與時間

GIM 資料 載點

• 連結：ftp://ftp.aiub.unibe.ch/CODE

GIM 檔名解釋

程式碼範例：只能繪製1小時的 GIM 資料

副程式：data_read.m

%======================================
% Read the GIM data every one row line
%======================================
function [data_count,rdline]=data_read(data,fseg,data_count);
data_count=data_count+1;
line_index=(fseg(data_count)+1):(fseg(data_count+1)-1);
rdline='                                                                                                      ';
rdline(1:length(line_index))=data(line_index);

主程式

%=======================
% Read the GIM TEC data
%=======================
clear all
close all

fn='CORG0850.18I';

fid=fopen(fn,'r');
data=fread(fid,inf);
fclose(fid);
fseg=find(data==10);
fseg=[0 fseg'];
data=char(data');
data_count=0;
headeryn=1;
count=1;

%----- Setting GIM parameters -----
gim_lat=87.5:-2.5:-87.5;
gim_lon=-180:5:180;

while data_count<(length(fseg)-1)
    [data_count,rdline]=data_read(data,fseg,data_count);
    %----- Read the GIM file header -----
    if headeryn==1
        if rdline(61:77)=='# OF MAPS IN FILE'
            map_num=str2num(rdline(1:6));
        end
        if rdline(61:73)=='END OF HEADER'
            headeryn=0;
        end
    else
        %----- Read GIM TEC data -----
        if rdline(61:71)=='END OF FILE'
            break
        end
        if rdline(61:76)=='START OF RMS MAP'
            break
        end
        if rdline(61:80)=='LAT/LON1/LON2/DLON/H'
            count2=1;
            for k=1:5
                [data_count,rdline]=data_read(data,fseg,data_count);
                tec_num=16;
                if k==5
                    tec_num=9;
                end
                for j=1:tec_num
                    index=1+(j-1)*5;
                    TEC(count,count2)=str2num(rdline(index:(index+4)))/10;
                    count2=count2+1;
                end
            end
            count=count+1;
        end
    end
end

TEC2=TEC(1:length(gim_lat),:);
[a,b]=contourf(gim_lon,gim_lat,TEC2,100);
set(b,'edgecolor','none')
caxis([0 50])
colorbar
axis equal
axis([-180 180 -90 90])
set(gca,'xtick',-180:30:180,'ytick',-90:30:90,'fontsize',12)
xlabel('Longitude (\circE)')
ylabel('Latitude (\circN)')
title('2018/03/26 0000UT')
%----- Plot the Global coast map -----
load coast
hold on
plot(long,lat,'k')

print -dpng GIM_01.jpg

結果

程式碼改寫：能繪製1天的 GIM 資料 (25張)

副程式：data_read.m

參考 程式碼範例：副程式：data_read.m

主程式

• 改寫說明：請參考程式碼註解的部分
• 大部分都只是程式碼邏輯的改寫
• 建議改寫的方法：
  • caxis([cmin cmax])：改變colorbar的上下邊界，理應按照需求調整邊界
• 額外加入的方法：

  • colormap(Colormap Name)：改變colorbar的顏色對應表

    Colormap Name	Color Scale
    parula (預設)
    jet (推薦使用)
    hsv
    hot

• 程式碼：

%=======================
% Read the GIM TEC data
%=======================
clear all
close all

% 輸入檔名
fn='CODG0010.15I';

fid=fopen(fn,'r');
data=fread(fid,inf);
fclose(fid);
fseg=find(data==10);
fseg=[0 fseg'];
data=char(data');
data_count=0;
headeryn=1;
count=1;

%----- Setting GIM parameters -----
gim_lat=87.5:-2.5:-87.5;
gim_lon=-180:5:180;

while data_count<(length(fseg)-1)
    [data_count,rdline]=data_read(data,fseg,data_count);
    %----- Read the GIM file header -----
    if headeryn==1
        if rdline(61:77)=='# OF MAPS IN FILE'
            map_num=str2num(rdline(1:6));
        end
        if rdline(61:73)=='END OF HEADER'
            headeryn=0;
        end
    else
        %----- Read GIM TEC data -----
        if rdline(61:71)=='END OF FILE'
            break
        end
        if rdline(61:76)=='START OF RMS MAP'
            break
        end
        if rdline(61:80)=='LAT/LON1/LON2/DLON/H'
            count2=1;
            for k=1:5
                [data_count,rdline]=data_read(data,fseg,data_count);
                tec_num=16;
                if k==5
                    tec_num=9;
                end
                for j=1:tec_num
                    index=1+(j-1)*5;
                    TEC(count,count2)=str2num(rdline(index:(index+4)))/10;
                    count2=count2+1;
                end
            end
            count=count+1;
        end
    end
end

% 加入迴圈，設定範圍從 0~24，共25張
for i=0:length(TEC)/length(gim_lat)-1
    % 一整天的資料筆數：length(TEC)
    % 每小時的資料筆數：length(gim_lat)
    TEC2=TEC(1+length(gim_lat)*i:length(gim_lat)*(i+1),:);
    [a,b]=contourf(gim_lon,gim_lat,TEC2,100);
    set(b,'edgecolor','none')
    caxis([0 50])
    colorbar
    % 改變 colorbar 的顏色樣式
    colormap('jet')
    axis equal
    axis([-180 180 -90 90])
    set(gca,'xtick',-180:30:180,'ytick',-90:30:90,'fontsize',12)
    xlabel('Longitude (\circE)')
    ylabel('Latitude (\circN)')
    
    % 改變每小時的圖片標題
    if i<10
        hr=['0' num2str(i)];
    else
        hr=num2str(i);
    end
    
    title(['2015/01/01 ' hr '00UT'])
    %----- Plot the Global coast map -----
    load coast
    hold on
    plot(long,lat,'k')
    
    % 改變輸出每小時的圖片檔名
    print('-dpng', ['GIM_' hr '.png'])
end

結果

• 因為25張圖過多，所以便將其製作成gif
  

Reference 參考資料

1. 高層大氣科學 上課講義
2. colormap 官方文件：https://www.mathworks.com/help/matlab/ref/colormap.html