【Matlab】GPS-TEC 的資料取得、轉換，並繪製變化圖

tags: Matlab 高層大氣科學

Data Source 資料來源

• Scripps Orbit and Permanent Array Center (SOPAC)：http://sopac.ucsd.edu/

• 進入官網後，可選擇以下兩個方式下載資料

  • Data Browser
  • Download Data：建議使用此種方式

  

Download Data 下載資料

1. 選擇使用 ftp

2. 登入方式

   • 帳號：anonymous
   • 密碼：你的 email 信箱

3. 進入以下資料夾：pub >> rinex >> 選擇年份 >> 選擇天數

   • 或者直接點選此連結：ftp://garner.ucsd.edu/pub/rinex/

4. 如果要大量下載data，請改用以下方式

   • 使用此軟體 FileZilla 進行下載：https://filezilla-project.org/download.php
   • 主機：填你要下載的資料夾網址、路徑，ex：ftp://garner.ucsd.edu/pub/rinex/2003/
   • 使用者名稱：也就是登入的帳號，anonymous
   • 密碼：也就是登入的密碼，你的 email 信箱
   • 填完後，按快速連線
   • 找到要下載的資料夾，按右鍵選擇下載，就能夠將所有資料夾中的檔案一次全部下載

   

   

Transform Data 資料轉換

• 轉換程式：下載
• 轉換步驟：
  1. 將剛才從 SOPAC 取得的資料和轉換程式放在同一層資料夾
  2. 開啟 Matlab
  3. 先執行 d2o.m，並輸入
     • 年分，ex：2003
     • 開始天數，ex：301
     • 結束天數，ex：301
  4. 再執行 gpstec2.m，並輸入
     • 年分，ex：2003
     • 開始天數，ex：301
     • 結束天數，ex：301
  5. 會得到很多檔案，主要會用到的有
     • a 開頭的 .dat檔：緯度
     • o 開頭的 .dat檔：經度
     • v 開頭的 .dat檔：垂直TEC
     • e 開頭的 .dat檔：天頂角
     • s 開頭的 .dat檔：斜向TEC
       

Data Format 資料格式

畫觀測資料位置

• 1筆資料
• a 開頭的 .dat檔：緯度
• o 開頭的 .dat檔：經度

Code 程式碼

% 0 means no data, so change 0 to nan
fn = 'ag101269.dat';
lats = load(fn);
f = find(lats == 0);
lats(f) = nan;

fn = 'og101269.dat';
lons = load(fn);
f = find(lons == 0);
lons(f) = nan;

plot(lons, lats, 'k.')
xlabel('Longitude (\circE)')
ylabel('Latitude (\circN)')

Result 結果

畫垂直 TEC 隨時間變化

• 1筆資料
• v 開頭的 .dat檔：垂直TEC

Code 程式碼

fn = 'vg101269.dat';
vtec = load(fn);
f = find(vtec == 0);
vtec(f) = nan;

plot((0:2879)/120, vtec)
xlabel('Time (UT)')
ylabel('VTEC (TECU)')

Result 結果

繪製全球垂直TEC分布

• 多筆資料：通常是同天數，所有測站的資料
• a 開頭的 .dat檔：緯度
• o 開頭的 .dat檔：經度
• v 開頭的 .dat檔：垂直TEC
• 小提醒：由於測站數量龐大，且每一個測站的資料量也很多，故跑一次(一張圖)，需要的時間約10分鐘上下，時間端看電腦效能而定。

Code 程式碼

clear all

% 找到所有v 開頭的 .dat檔
list = dir('v*.dat');

% 每兩小時繪製一張圖，共13張
for t=0:12
    close all
    lons=[];
    lats=[];
    vtec=[];
    % 每30秒一筆資料，故2小時就是相差2*120筆資料
    time_index = 1+t*2*120;
    % 最後一張圖，資料會不足2*120筆資料
    if time_index > 2880
        time_index = 2880;
    end
    for i=1:length(list)
        fn=list(i).name;
        tmp = load(['o' fn(2:end)]);
        if ~isempty(tmp)
            lons = [lons tmp(time_index, :)];
            tmp = load(['a' fn(2:end)]);
            lats = [lats tmp(time_index, :)];
            tmp = load(['v' fn(2:end)]);
            vtec = [vtec tmp(time_index, :)];
        end
    end
    
    f = find(lons==0 | lats==0 | vtec==0);
    lons(f)=[];
    lats(f)=[];
    vtec(f)=[];
    gmap
    caxis([0 70])
    colorbar
    colormap('jet')
    hold on
    % 5 means dot scale
    scatter(lons, lats, 5, vtec, 'filled');
    shading flat
    
    if t*2<10
        tt=['0' num2str(t*2)];
    else
        tt=num2str(t*2);
    end
    
    title(['2011/09/26 ' tt '00UT'])
    print('-dpng', ['GPS TEC_' tt '00UT.png'])
end

Result 結果

Reference 參考資料

1. 高層大氣科學 上課講義
2. Scripps Orbit and Permanent Array Center (SOPAC)：http://sopac.ucsd.edu/
3. FileZilla：https://filezilla-project.org/download.php