Архив

master
L_DelOff 2021-06-12 23:28:40 +03:00
parent a9f20b053a
commit db195a3307
1 changed files with 410 additions and 0 deletions

410
loaddata_edr_old.m Normal file
View File

@ -0,0 +1,410 @@
function loaddata_edr_old %DATA=loaddata_edr(filename)
%% By L_DelOff
% Загрузка данных из файла *.EDR
% Описание файла приведено в:
% !docs/!AFRL ASCII and Binary File Format Descriptions.pdf
%% Поехали
warning off
sattelite_range=15:18; % Выбор с какого спутника требуются данные
year_range=2015; %
month_range=9; %
%% Скачать данные в формате архива
if 1
for sattelite=sattelite_range %6:20
for year=year_range % года
for month=month_range % месяца
download_fun(sattelite, year, month); %скачивание и обработка нужных данных
end
end
end
end
%% Распаковать архив
if 1
for sattelite=sattelite_range %6:20
for year=year_range % года
for month=month_range % месяца
extract_data_fun(sattelite, year, month); %скачивание и обработка нужных данных
end
end
end
end
%% Преобраховать EDR в mat
if 1
for sattelite=sattelite_range %6:20
for year=year_range % года
for month=month_range % месяца
convert_data_fun(sattelite, year, month); %скачивание и обработка нужных данных
end
end
end
end
end
function convert_data_fun(sattelite, year, month)
type_sensor='ssies';
foldername=['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)];
target=ls(foldername);
[str row]=size(target);
for i=1:str
filename=target(i,:);%end-3:end);
counter=0;
for j=1:length(filename)
if filename(j)==' '
counter=counter+1;
end
end
filename=filename(1:end-counter);
if length(filename)>4
if filename(end-3:end)=='.EDR'
message=['Расшифровка файла:',filename,'\n'];
fprintf(message);
fid = fopen([foldername,'/',filename], 'rt'); % открыть для чтения
counter=1;
while ~feof(fid) % пока не конец файла
frame=[]; % здесь будет 1 кадр (1 минута) в текстовом формате
for i=1:114
if ~feof(fid)
frame(i).str = fgets(fid); % Последовательно читаем из файла строки
end
end
if counter==1
DATA=getframe_edr(frame); % Считываем кадр
else
DATA(counter)=getframe_edr(frame); % Считываем кадр
end
counter=counter+1;
if mod(counter,100)==0
disp(num2str(counter));
end
end
save([foldername,'/',filename(1:end-4),'.mat'],'DATA');
message=['Расшифровка файла закончена\n'];
fprintf(message);
end
end
end
end
function extract_data_fun(sattelite, year, month)
type_sensor='ssies';
foldername=['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)];
target=ls(foldername);
[str row]=size(target);
for i=1:str
if target(i,end-2:end)=='.gz'
message=['Распаковка:',target(i,:),'\n'];
fprintf(message);
gunzip([foldername,'/',target(i,:)],foldername);
message=['Распакован\n'];
fprintf(message);
end
end
end
function download_fun(sattelite, year, month)
% https://satdat.ngdc.noaa.gov/dmsp/data/f18/ssies/2020/05/
% https://satdat.ngdc.noaa.gov/dmsp/data/f18/ssies/2020/05/PS.CKGWC_SC.U_DI.A_GP.SIES3-F18-R99990-B9999090-APGA_AR.GLOBAL_DD.20200501_TP.000001-235959_DF.EDR.gz
type_sensor='ssies';
message=['Ищем данные спутника F' num2str(sattelite) ' по дате:' num2str_new(year,4) '/' num2str_new(month,2) '...'];
fprintf(message)
url_site=['https://satdat.ngdc.noaa.gov/dmsp/data/f' num2str_new(sattelite,2) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2) '/'];
fullList = webread(url_site);
message='готово\n';
fprintf(message)
url=extract_url(fullList);
for i=1:length(url)
message=['число ' num2str(i) '...'];
fprintf(message);
url_file = [url_site url(i).file];
filename = ['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2) '/' url(i).file];
mkdir(['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)])
outfilename = websave(filename,url_file);
message='скачано\n';
fprintf(message);
end
a=1;
end
function url=extract_url(text)
counter=0;
i=1;
while i<length(text)
if text(i)=='"'
a=i;
i=i+1;
while text(i)~='"'
i=i+1;
end
b=i;
if b-a>10
c=text(b-6:b-1);
if c=='EDR.gz'
counter=counter+1;
url(counter).file=text(a+1:b-1);
end
end
end
i=i+1;
end
end
% %% Инициализация
% %load DATA
% %% Открытие файла если уже нет DATA
% if exist('DATA')~=1
% fid = fopen('DATA/PS.CKGWC_SC.U_DI.A_GP.SIES3-F18-R99990-B9999090-APGA_AR.GLOBAL_DD.20140831_TP.000001-235959_DF.EDR', 'rt'); % открыть для чтения
% counter=1;
% while ~feof(fid) % пока не конец файла
% frame=[]; % здесь будет 1 кадр (1 минута) в текстовом формате
% for i=1:114
% if ~feof(fid)
% frame(i).str = fgets(fid); % Последовательно читаем из файла строки
% end
% end
%
% if counter==1
% DATA=getframe_edr(frame); % Считываем кадр
% else
% DATA(counter)=getframe_edr(frame); % Считываем кадр
% end
% counter=counter+1;
% if mod(counter,100)==0
% disp(num2str(counter));
% end
% end
% end
%
% for i=1:length(DATA)
% x(i)=(DATA(i).GeoLong(1));
% y(i)=(DATA(i).GeoLat(1));
% end
% plot(x,y)
% xlabel('Geographic longitude (degrees, east)')
% ylabel('Geographic latitude (degrees, north)')
%end
function FRAME=getframe_edr(frame)
HEALTH=0; % "здоровье кадра"
%% строка 1: Пустая строка (Blank line)
if isspace(frame(1).str) HEALTH=HEALTH+1;
end
%% строка 2: Шапка + данные о ПО(если это первый кадр в файле)
if frame(2).str~="" HEALTH=HEALTH+1; end
%% строка 3: 5 значений
z = textscan(frame(3).str,'%f','EmptyValue',-Inf);
FRAME.N_record_bin = z{1}(1); %1. Номер записи в двоичном файле, откуда был взят этот EDR.(Number of the record in the binary file from where this EDR was taken.)
FRAME.N_EDR_bin = z{1}(2); %2. Номер EDR в записи в двоичном файле (1-3).(Number of the EDR within the record in the binary file (1-3).)
FRAME.SatID = z{1}(3); %3. ID спутника.(Satellite Flight ID (two digit integer))
FRAME.year = fix(z{1}(4)/1e4); %4. Год, месяц, день, час, минута
FRAME.month = mod(fix(z{1}(4)/1e2),1e2);
FRAME.day = mod(z{1}(4),1e2);
FRAME.hour = fix(z{1}(5)/1e2);
FRAME.minute = mod(z{1}(5),1e2);
if length(z{1})==5 HEALTH=HEALTH+1; end
%% строка 4: ЭФЕМЕРИДА (EPHEMERIS)
if frame(4).str(1:end-1)=="EPHEMERIS" HEALTH=HEALTH+1; end
%% строки 5-7: по 6 значений
for i=0:2
a = textscan(frame(5+i).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
% данные на секунде 00 20 40 (i=0 1 2)
FRAME.GeoLat(i+1)=a{1}(1); % 1. Географическая широта (градусы, север)(Geographic latitude (degrees, north))
FRAME.GeoLong(i+1)=a{1}(2); % 2. Географическая долгота (градусы, восток)(Geographic longitude (degrees, east))
FRAME.ApLat(i+1)=a{1}(3); % 3. Вершинная широта (градусы, север)(Apex latitude (degrees, north))
FRAME.ApLong(i+1)=a{1}(4); % 4. Вершинная долгота (градусы, восток)(Apex longitude (degrees, east))
FRAME.ApLocalTime(i+1)=a{1}(5);% 5. Местное время Apex (часы)(градусы, восток)(Apex local time (hours))
FRAME.SatAlt(i+1)=a{1}(6); % 6. Высота (км)(Satellite altitude (km))
end
%% строка 8: Потенциал спутника (“SATELLITE POTENTIAL, LAST = SOURCE”)
if frame(8).str(1:end-1)=="SATELLITE POTENTIAL, LAST = SOURCE" HEALTH=HEALTH+1; end
%% строки 9-10: 16 значений (2х8)
a = textscan(frame(9).str,'%f','EmptyValue',-Inf);
b = textscan(frame(10).str,'%f','EmptyValue',-Inf);
if length(a{1})==8 HEALTH=HEALTH+1; end
if length(b{1})==8 HEALTH=HEALTH+1; end
FRAME.SatPot(1:8)=a{1}(1:8); % значения для потенциала спутника (Vbias + VIP) в вольтах каждые 4 секунды
FRAME.SatPot(9:15)=b{1}(1:7); % values for the satellite potential (Vbias+VIP) in volts every 4 seconds
FRAME.SatPotSrc=b{1}(8); % Измеритель потенциала (an integer (1-2) indicating the satellite potential source where:
% 1 - as set by the on-board microprocessor
% 2 - as set by the SENPOT sensor)
%% строка 11: Потенциал спутника (“PRIMARY PLASMA DENSITY, THEN SOURCE”)
if frame(11).str(1:end-1)=="PRIMARY PLASMA DENSITY, THEN SOURCE" HEALTH=HEALTH+1; end
%% строки 12-21: 60 значений (6х10)
% Односекундные средние значения плотности первичной плазмы (# / см3)
% One-second averages of the primary plasma density (#/cm3)
z=[];
for i=12:21
a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
z=[z a{1}.'];
end
FRAME.PrimPlasmaDen=z;
%% строка 22: Источник данных по плазме
% An integer (1-3) indicating the plasma density source where:
% 1 - Ion density from SM sensor
% 2 - Ion density from DM sensor
% 3 - Electron density from EP sensor (DC Mode)
a = textscan(frame(22).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.PrimPlasmaDenSrc=a{1};
%% строка 23: Горизонтальный ионный дрейф (“HORIZONTAL ION DRIFT VELOCS”)
if frame(23).str(1:end-1)=="HORIZONTAL ION DRIFT VELOCS" HEALTH=HEALTH+1; end
%% строки 24-33: 60 значений (6х10)
% Односекундные значения скорости горизонтального дрейфа (м/с) для времен от HHMM00 до HHMM59 с 6 значениями в строке и всего 60 значениями.
% One-second values for the horizontal drift speed (m/s) for times HHMM00 through HHMM59 with 6 values per line and 60 values total.
z=[];
for i=24:33
a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
z=[z a{1}.'];
end
FRAME.HorIonDrift=z;
%% строка 34: Вертикальный ионный дрейф (“VERTICAL ION DRIFT VELOCS”)
if frame(34).str(1:end-1)=="VERTICAL ION DRIFT VELOCS" HEALTH=HEALTH+1; end
%% строки 35-44: 60 значений (6х10)
% Односекундные значения скорости вертикального дрейфа (м/с) для времен от HHMM00 до HHMM59 с 6 значениями в строке и всего 60 значениями.
% One-second values for the vertical drift speed (m/s) for times HHMM00 through HHMM59 with 6 values per line and 60 values total.
z=[];
for i=35:44
a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
z=[z a{1}.'];
end
FRAME.VerIonDrift=z;
%% строка 45: (“CKL ANALYSES, THEN SOURCE”)
if frame(45).str(1:end-1)=="CKL ANALYSES, THEN SOURCE" HEALTH=HEALTH+1; end
%% строка 46-63:
for i=1:6 % для секунд: 5 15 25 35 45 55
z=[];
a = textscan(frame(46+(3*(i-1))).str,'%f','EmptyValue',-Inf);
if length(a{1})==4 HEALTH=HEALTH+1; end
FRAME.CKL_RMS(i)=a{1}(1); % CKL Analysis: (RMS ΔN)/N (%)
FRAME.CKL_T1(i)=a{1}(2); % CKL Analysis: T1
FRAME.CKL_p1(i)=a{1}(3); % CKL Analysis: p1
FRAME.CKL_CKL(i)=a{1}(4); % CKL Analysis: CKL
a = textscan(frame(47+(3*(i-1))).str,'%f','EmptyValue',-Inf);
if length(a{1})==8 HEALTH=HEALTH+1; end
z=[z a{1}.'];
a = textscan(frame(48+(3*(i-1))).str,'%f','EmptyValue',-Inf);
if length(a{1})==8 HEALTH=HEALTH+1; end
z=[z a{1}.'];
FRAME.CKL_PDS(i,:)=z(1:end-1); %CKL Analysis: Decimated power density spectrum (PDS) for time period.
end
%% строка 64: Источник данных по CLK
% An integer (1-3) indicating data source used for CKL calculation where:
% 1 - SM density data only
% 2 - SM density and filter data
% 3 - EP DC mode density data
a = textscan(frame(64).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.CLK_Src=a{1};
%% строка 112 (потом понадобится)
a = textscan(frame(112).str,'%f','EmptyValue',-Inf);
if length(a{1})==7 HEALTH=HEALTH+1; end
FRAME.ADC_temp=a{1}(2); %2: ADC temperature (degrees C)
FRAME.SEP_temp=a{1}(3); %3: SEP temperature
FRAME.DM_offsetV=a{1}(4); %4: DM offset voltage (volts) if IES-2 or RPA plasma plate potential (volts) if SSIES-3
FRAME.DM_mode=a{1}(5); %5: DM mode (integer)
% SSIES-2:
% -1: Undefined
% 0: Normal
% 1-8: H+
% 9: FIBA
% SSIES-3:
% 0: Slow
% 1: Normal
FRAME.EP_mode=a{1}(6); % 6: EP mode (0-6 : A,B,BS,C,D,DS,E) (integer)
FRAME.EP_VIP=a{1}(7); % 7: VIP at EDR start (volts)
%% строки 65-80
if (FRAME.EP_mode==0)||(FRAME.EP_mode==1)||(FRAME.EP_mode==2)||(FRAME.EP_mode==6)
% строка 65 (“EP SWEEP ANALYSES SETS”)
if frame(65).str(1:end-1)=="EP SWEEP ANALYSES SETS" HEALTH=HEALTH+1; end
% EP Sweep analyses (every 4 seconds). There are 15 EP sweep analysis
% sets. Each is valid for either 4 (modes A, B and BS) or 2 (mode E)
% seconds centered on the time specified in the set. Each analysis set
% contains the following parameters:
for i=1:15
a = textscan(frame(66+(i-1)).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
FRAME.EP_SCT(i)=a{1}(1); % 1: Sweep center time (UT, seconds) (integer)
FRAME.EP_Eden(i)=a{1}(2); % 2: Electron density (el/cm3)
FRAME.EP_Etemp(i)=a{1}(3); % 3: Electron temperature (degrees K)
FRAME.EP_SatPot(i)=a{1}(4); % 4: Satellite potential (volts)
FRAME.EP_Qualif(i)=a{1}(5); % 5: Analysis qualifier (integer). Set to zero if the on-board microprocessor did not perform the analysis, per the flag in element 415. If the on-board microprocessor was in use, then it is set to the MP EP flags word from Word 11 of Cycle 1 in the telemetry. This word can also be zero.
FRAME.EP_PES(i)=a{1}(6); % 6: EP photo-electron surrogate value.
end
end
%% строка 81 (“EP ANALYSES SOURCE”)
if frame(81).str(1:end-1)=="EP ANALYSES SOURCE" HEALTH=HEALTH+1; end
%% строка 82
% An Integer (1-2) indicating EP analysis source where
% 1 Ground processing analysis
% 2 On-board microprocessor analysis
a = textscan(frame(82).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.EP_src=a{1};
%% строка 83 (“RPA SWEEP ANALYSES SETS, THEN SOURCE”)
if frame(83).str(1:end-1)=="RPA SWEEP ANALYSES SETS, THEN SOURCE" HEALTH=HEALTH+1; end
%% строка 84-98 (“RPA SWEEP ANALYSES SETS, THEN SOURCE”)
% RPA Sweep analyses (every 4 seconds). There are 15 RPA sweep analysis
% sets. Each is valid for the 4 seconds centered on the time specified in
% the set. Each analysis set contains the following parameters:
for i=1:15
a = textscan(frame(84+(i-1)).str,'%f','EmptyValue',-Inf);
if length(a{1})==8 HEALTH=HEALTH+1; end
FRAME.RPA_SCT(i)=a{1}(1); % 1: Sweep center time (UT, seconds) (integer)
FRAME.RPA_O2den(i)=a{1}(2); % 2: O+ density (ion/cm3)
FRAME.RPA_HHeden(i)=a{1}(3); % 3: Total (H+ + He+) density (ion/cm3)
FRAME.RPA_LIF(i)=a{1}(4); % 4: Light ion flag (integer)
% 0 - No light ion
% 1 - Light ion is H+
% 2 - Light ion is He+
% 3± = 3 + 10000 x (H+ fraction)
FRAME.RPA_IonNemp(i)=a{1}(5); % 5: Ion temperature (degrees K)
FRAME.RPA_RamIonDrift(i)=a{1}(6);% 6: Ram ion drift velocity (m/s)
FRAME.RPA_qualif(i)=a{1}(7); % 7: Analysis qualifier (integer)
% 0 - Analysis terminated unsuccessfully
% 1 - Successful analysis
FRAME.RPA_totIonDens(i)=a{1}(8);% 8: RPA-derived total ion density
end
%% строка 99
% An integer (1-2) indicating RPA analysis source where:
% 1 - Ground processing analysis
% 2 - On-board microprocessor analysis
a = textscan(frame(99).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.RPA_src=a{1};
%% строка 100 (“DM ION DENSITY”)
if frame(100).str(1:end-1)=="DM ION DENSITY" HEALTH=HEALTH+1; end
%% строки 101-110: 60 значений (6х10)
%DM one-second average ion density (ion/cm3) with 6 values per line and 60
% values total
z=[];
for i=101:110
a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
if length(a{1})==6 HEALTH=HEALTH+1; end
z=[z a{1}.'];
end
FRAME.DM_IonDen=z;
%% строка 111 (“ENGINEERING DATA”)
if frame(111).str(1:end-1)=="ENGINEERING DATA" HEALTH=HEALTH+1; end
%% строка 113 (“FILLER”)
if frame(113).str(1:end-1)=="FILLER" HEALTH=HEALTH+1; end
%% строка 114 (“FILLER”)
a = textscan(frame(114).str,'%f','EmptyValue',-Inf);
if length(a{1})>=1 HEALTH=HEALTH+1; end
FRAME.Filler=a{1};
FRAME.HEALTH=HEALTH;
end