【Matlab】地磁觀測資料分析 - subplot, xticklabel, yticklabel的應用

tags: Matlab 高層大氣科學

目標

1. 下載任意觀測站任意一個月的地磁觀測資料

2. 繪製H-component一個月的變化情形，格式如下圖所示 (範例圖)

3. 於圖上標上測站名稱、觀測年份與月份。

   

地磁觀測資料 載點

1. 全球：http://www.intermagnet.org/index-eng.php
2. 日本：http://www.kakioka-jma.go.jp/obsdata/metadata/?locale=en

地磁資料格式

程式碼範例：只能繪製1天的地磁資料

副程式：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 Magnetometer data
%============================
clear all
close all

fn='kak20180301dmin.min';

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;

while data_count<(length(fseg)-1)
    [data_count,rdline]=data_read(data,fseg,data_count);
    %----- Read the file header -----
    if headeryn==1
        if rdline(1:4)=='DATE'
            headeryn=0;
        end
    else
       %----- Read Magnetometer data -----
       [data_count,rdline]=data_read(data,fseg,data_count);
       obs_time(count)=str2num(rdline(12:13))+str2num(rdline(15:16))/60;
       obs_D(count)=str2num(rdline(31:40));
       obs_H(count)=str2num(rdline(41:50));
       obs_Z(count)=str2num(rdline(51:60));
       obs_F(count)=str2num(rdline(61:70));
       count=count+1;
    end
end
plot(obs_time,obs_H,'r.-')
set(gca,'linewidth',2,'fontsize',12)
set(gca,'xtick',0:2:24,'xlim',[0 24])
xlabel('Time (UT)')
ylabel('H component')
title('Kakioka 2018/03/01')

結果

程式碼改寫：能繪製1個月的地磁資料

副程式：data_read.m

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

主程式

• 改寫說明：請參考程式碼註解的部分
• 大部分都只是程式碼邏輯的改寫
• 額外加入的方法：
  • subplot(m,n,p)：設定子圖
  • set(gca,'YTickLabel',[])：設定Y軸標籤上的資料
  • set(gca,'XTickLabel', 1+i*10:1:12+i*10)：設定X軸標籤上的資料
  • set(gca,'XGrid','on')：設定X軸上的網格是否開啟
• 程式碼：

%============================
% Read the Magnetometer data
%============================
clear all
close all

% 加入 subplot，讓3張圖能同時繪製在1張圖上
for i=0:2
    subplot(3,1,i+1)
    % 每張 subplot 有11筆資料 (也就是要讀檔案11次)
    for j=1+i*10:11+i*10
        
        % 由於檔案名稱中的日期是 01~09，所以要先將1~9轉換成01~09
        if j<10
            date=['0' num2str(j)];
        else
            date=num2str(j);
        end
        fn=['kak201707' date 'dmin.min'];
        
        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;
        
        while data_count<(length(fseg)-1)
            [data_count,rdline]=data_read(data,fseg,data_count);
            %----- Read the file header -----
            if headeryn==1
                if rdline(1:4)=='DATE'
                    headeryn=0;
                end
            else
                %----- Read Magnetometer data -----
                [data_count,rdline]=data_read(data,fseg,data_count);
                % rdline(12:13) 是獲得時
                % rdline(15:16))/60 也是將剩餘的分和秒轉換成時
                % 從上述可知每個檔案(1天)，obs_time(count)都一樣
                % 所以必須每天+24hr，累加時間上去，才能將11筆資料畫在同一個subplot上
                obs_time(count)=str2num(rdline(12:13))+str2num(rdline(15:16))/60+(j-1)*24;
                obs_D(count)=str2num(rdline(31:40));
                obs_H(count)=str2num(rdline(41:50));
                obs_Z(count)=str2num(rdline(51:60));
                obs_F(count)=str2num(rdline(61:70));
                count=count+1;
            end
        end
        
        % 每讀1次檔案，也就是每1天，繪製1次圖
        plot(obs_time,obs_H,"r.-")
        hold on
    end
    % 設定 title 在整張圖的最上方
    if i==0
        title('Kakioka 2017/07')
    end
    % 設定 x label 在整張圖的最下方
    if i==2
        xlabel('Date')
    end
    % 每張 subplot，設定 linewidth, fontsize, y label, YTickLabel, xtick, xlim, XTickLabel, XGrid
    set(gca,'linewidth',2,'fontsize',14)
    ylabel('H component')
    set(gca,'YTickLabel',[]);  
    set(gca,'xtick',i*240:24:264+i*240,'xlim',[i*240 264+i*240],'XTickLabel', 1+i*10:1:12+i*10,'XGrid','on')
end

結果

Reference 參考資料

1. 高層大氣科學 上課講義