Source Code for 1-min MAGDAS data

版本 1.0.0 (9.6 KB) 作者: Godfrey Ojerheghan
This code can reconstruct the H, D, F, Z magnetogram from the .mgd files
0.0
(0)
18.0 次下载
更新时间 2021/8/19

查看许可证

function [MagData,Header,STATDATA,CORRECT_INF]=read_1M(FileName)
%============================================================================
%
% [MagData,Header,STATDATA,CORRECT_INF]=read_1M(FileName)
%
% Read 1-Minute data of MAGDAS staraged data
%
% Input: FileName
% Outout: MagData : Magnetic Data [H,D,Z,F,Inc_EW,Inc_NS,Temp] (1440x7/1day)
% Header : Header Information (Structure array)
% STATDATA : Status Data of the Magnetometer (Structure Array)
% CORRECT_INF : GPS Time Correct Information (Structure Array)
%
% 2003/10/09 Ver. 1.0 by K. Kitamura
% 2004/07/22 Ver. 1.1 (時???Z?ｳ回?狽窿wッダ?[の繰り返し０に対応?j
%============================================================================
%clear all
FileName=('ABU_MIN_201102060000.mgd')
%================= Open file ===================
[Fid,Message]=fopen(FileName,'r','l');
if Fid<0,
disp(Message)
return
end
% ============== Find magnetic data ======================
ReadUnit=1000; FindMag=0; NumRead=0;
Status=fseek(Fid,0,'bof');
while FindMag==0,
[Buff,Count]=fread(Fid,ReadUnit,'uchar');
[Message,Errnum]=ferror(Fid);
if Errnum~=0,
error(Message)
end
NumRead=NumRead+1;
Find1A=find((Buff-26)==0);
if ~isempty(Find1A),
for i=1:length(Find1A),
if Find1A(i)~=ReadUnit,
if Buff(Find1A(i)+1)==0,
MagPos=(NumRead-1)*ReadUnit+Find1A(i)+1;
FindMag=1; break;
end
else
[Buff,Count]=fread(Fid,ReadUnit,'uchar');
[Message,Errnum]=ferror(Fid);
if Errnum~=0,
error(Message)
end
NumRead=NumRead+1;
if Buff(1)==0,
MagPos=(NumRead-1)*ReadUnit+1;
FindMag=1; break;
end
end
end
end
end
fclose(Fid);
% ============== Read header =============================
[Fid,Message]=fopen(FileName,'rt');
if Fid<0,
disp(Message)
return
end
Status=fseek(Fid,0,'bof');
head=fgetl(Fid);
Header.DATA_TYPE=sscanf(head(10:end),'%s');
head=fgetl(Fid);
Header.SERIAL_NUM=str2num(head(11:end));
head=fgetl(Fid);
Header.STATION_NAME=sscanf(head(13:end),'%s');
head=fgetl(Fid);
Header.GEODETIC_LATITUDE=str2num(head(18:end));
head=fgetl(Fid);
Header.GEODETIC_LONGITUDE=str2num(head(19:end));
head=fgetl(Fid);
Header.RECORD_TIME=sscanf(head(12:end),'%s');
head=fgetl(Fid);
Header.REPORTED=sscanf(head(9:end),'%s');
head=fgetl(Fid);
Header.SAMPLING_TIME=str2num(head(14:end));
head=fgetl(Fid);
Header.MAG_RANGE=str2num(head(10:end));
head=fgetl(Fid);
Header.HEADER_DATA_NUM=str2num(head(18:end));
if Header.HEADER_DATA_NUM~=0
for i=1:Header.HEADER_DATA_NUM
head=fgetl(Fid);
STATDATA(i).HEADER_DATA_TIME=sscanf(head(17:end),'%s');
head=fgetl(Fid);
STATDATA(i).BAIAS_H=str2num(head(8:end));
head=fgetl(Fid);
STATDATA(i).BAIAS_D=str2num(head(8:end));
head=fgetl(Fid);
STATDATA(i).BAIAS_Z=str2num(head(8:end));
head=fgetl(Fid);
STATDATA(i).BAIAS_F=str2num(head(8:end));
end
else
STATDATA=[];
end
head=fgetl(Fid);
Header.CORRECT_NUM=str2num(head(12:end));
if Header.CORRECT_NUM~=0
for i=1:Header.CORRECT_NUM
head=fgetl(Fid);
CORRECT_tmp=sscanf(head(12:end),'%s');
CORRECT_INF(i).TIME=CORRECT_tmp(1:8);
CORRECT_INF(i).VALUE=str2num(CORRECT_tmp(10:end));
end
else
CORRECT_INF=[];
end
fclose(Fid);
% HEADER{1}=Header;
% HEADER{2}=STATDATA;
% HEADER{3}=CORRECT_INF;
%================== Read Data ====================
[Fid,Message]=fopen(FileName,'r','l');
% ============== Calculate size of magnetic data =========
Status=fseek(Fid,0,'eof');
EOFid=ftell(Fid);
DataSize=EOFid-MagPos;
%========================================
Status=fseek(Fid,MagPos,'bof');
[MagData,Count]=fread(Fid,[7,DataSize/28],'float32');
%========== convert error data (NaN) to 99999.99=================
[errorrow,errorcol]=find(isnan(MagData)==1);
MagData(errorrow,errorcol)=999999.9;
MagData=MagData';
fclose(Fid);
thr1= [0:60:1440];
thr1= [ 60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440];
subplot(4,1,1)
plot(MagData(:,1))
set(gca,'XTick',[thr1])
set(gca,'XLim',[0 1440])
set (gca,'FontSize',10)
xlabel('Local time (Hrs)','FontSize',10)
ylabel('H (nT)','FontSize',10)
title('ABU\it{H(nT) 06 Feb 2011}','FontSize',14)
%title('\it{Diurnal variation of Day to day variability of H values 061106/07}','FontSize',14)
%
subplot(4,1,2)
plot(MagData(:,2))
set(gca,'XTick',[thr1])
set(gca,'XLim',[0 1440])
set (gca,'FontSize',10)
xlabel('Local time (Hrs)','FontSize',10)
ylabel('D(nT)','FontSize',10)
title('ABU\it{D(nT) 06 Feb 2011}','FontSize',14)
%title('\it{Diurnal variation of Day to day variability of D values 061106/07}','FontSize',14)
%
subplot(4,1,3)
plot(MagData(:,3))
set(gca,'XTick',[thr1])
set(gca,'XLim',[0 1440])
set (gca,'FontSize',10)
xlabel('Local time (Hrs)','FontSize',10)
ylabel('Z(nT)','FontSize',10)
title('ABU\it{Z(nT) 06 Feb 2011}','FontSize',14)
subplot(4,1,4)
plot(MagData(:,4))
set(gca,'XTick',[thr1])
set(gca,'XLim',[0 1440])
set (gca,'FontSize',10)
xlabel('Local time (Hrs)','FontSize',10)
ylabel('Z(nT)','FontSize',10)
title('ABU\it{Z(nT) 06 Feb 2011}','FontSize',14)
dummy=GenerateMin(MagData);
Hmin=dummy(:,1);
Dmin=dummy(:,2);
Zmin=dummy(:,3);
Fmin=dummy(:,4);
Mag_dat=[Hmin;Dmin;Zmin;Fmin];
% [nrow,ncols]= size(Mag_dat);
Mag_data=(reshape(Mag_dat, length(Mag_dat)/4, 4));
% ddata1=[ddata1;(reshape(data1, ncols, length(data1)/ncols))']
Mag_data=[thr1', Mag_data];
save H06Feb2011ABU.mat;
%title('\it{Diurnal variation of Z values 061106/07}','FontSize',14)
%save dummy.dat Hmin Dmin Zmin Fmin -ascii
%----------------------------------------
% %efforts to evaluate hourly means
% %-----------------------------------
function y=GenerateMin(x)
% Reduce data to every minute average
deltat=60;
Hrpday=24;
for i=1:Hrpday,
for j=1:4,
aux=x((i-1).*deltat+1:i.*deltat,j); % Take every 60 seconds data
I=find(not(isnan(aux)));
if length(I)>=10,
y(i,j)=mean(aux(I));
end
end
end

引用格式

Godfrey Ojerheghan (2024). Source Code for 1-min MAGDAS data (https://www.mathworks.com/matlabcentral/fileexchange/97884-source-code-for-1-min-magdas-data), MATLAB Central File Exchange. 检索时间: .

MATLAB 版本兼容性
创建方式 R2021a
兼容任何版本
平台兼容性
Windows macOS Linux
标签 添加标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

版本 已发布 发行说明
1.0.0