Verfasst am: 28.10.2018, 12:16
Titel: Index Problem, finde par tout nicht das Problem
Hallo zusammen,
ich arbeite an Klimaentwicklungsmodellen und nutze hierfür auch Kohlenstoffdioxid (CO2) und Methan (CH4) Daten. Dieses Jahr war ich wieder mit zwei Spektrometern unterwegs und habe zwei schöne, große Datensätze gesammelt.
Was ich von meiner Wishlist abgehakt habe:
-Die Daten von Analyzer 1 werden erkannt und eingelesen
-Die Daten von Analyzer 2 werden erkannt und eingelesen
-Die Daten von der Masterclock (zum Synchronisieren der Daten) "Autochamber" werden erkannt und eingelesen
-Die voreingestellten Kalibrationen werden auf die Datensätze angewandt.
Ich stecke jetzt an dem Punkt, an dem es darum geht,
-die Daten aller drei Geräte zu synchronisieren
-die Ergebnisse als csv zu exportieren, um anderweitig damit zu arbeiten.
Konkret heißt das, dass ich bei
Code:
% Fehler: Index exceeds array bounds, tritt auf in metPick=metPick(1); % Get the met variables to calculate the flux from
metPick = find(tv_mb >= tv(i));
%grab the closest met perameters which occur after measurment (within 30min)
metPick = metPick(1);
Tair_AC=Tair(metPick);
Pb_AC=Pb(metPick);
WtrLv_Met=WT(metPick);
nicht verstehe, welcher Index hier über welches Array hinausgeht
und zweitens, wie ich am schlauesten die drei Geräte synchronisiere...
Vielen lieben Dank, wenn jemand Lust hat, an diesem umfangreicheren Projekt mitzubasteln.
Da ich mit internationalen Kollegen an diesen Daten arbeiten möchte, ist alles auf englischer Sprache. Falls ihr Fragen habt, immer zu.
Code:
% Mein bisheriger Code: function datT=RwaDataProcessing_PG(fromRaw) % This function will load chamber data based on Elyns process_MB_chamber code. % Inputs: fromRaw = a true or false flag for loading data from files % Outputs: datT = a matlab table which holds the fluxes and data
unzip=false; %set to trueif the files are zipped.
rootdir = 'P:\'; %Set root-directory; gets overwritten below
yr=2018;
tvoffset_a=0; %This offset should be zero
tvoffset_i= 0 %-0.294213; %This is the temporal offset between the isotopes and the autochambers
tvoffset_mci= 0 %0.076; %This is the temporal offset between the mcia3_8 files and the autochambers
if fromRaw
[Tair,par,Kin,RH,Pb,ustar,WT,Tsoil,Tsoil_all,tv_mb]=LoadECData(false);
warning off %supress warnings (there are a few when loading the data)
d = dir([rootdir '1st Run\Incoming\']);
if unzip %should we unzip the data first?
UnzipData(rootdir);
end
fprintf('Isotopes:') for i = 3:length(d) waitbar(i/length(d),f);
filename = d(i).name;
ifstrcmp(filename(1:3),'CO2')% String-compare: If the file is a CO2 file
%load in first column to find start of extra-text at bottom of file and set stopline:
tmp = importfile_MB_chamber_isotope_stopline([rootdir '1st Run\Incoming\' filename]);
stopline = strmatch('-----BEGIN',strvcat(tmp)); clear tmp;
%if no text at end of file, set stopline to infinite
ifisempty(stopline); stopline = inf; end;
try; [Time,CO2_ppm,foo,D13C_VPDB_CO2,foo,foo,foo,H2O_ppm] = ...
importfile_MB_chamber_isotope_2018([rootdir '1st Run\Incoming\' filename],3,stopline);
%add error msg for files with issues
catch; disp(['Could not process ' rootdir '1st Run\Incoming\' filename]); end;
%generate datevecDate Vectors with offset added:
tv_dv = datevec(Time+tvoffset_i);
%round the seconds to the nearest 10 !!!:
tv_dv_rnd = datevec(datenum([tv_dv(:,1:5)round(tv_dv(:,6)/10).*10]));
tv_rnd = datenum(tv_dv_rnd);
%find where the data fits:
[c,ia,ib] = intersect(tv_rnd,tv);
%dump loaded data into the master chonologies:
CO2_i(ib) = CO2_ppm(ia);
D13C(ib) = D13C_VPDB_CO2(ia);
H2O(ib) = H2O_ppm(ia);
% %this is just a progress bar thing so it doesn't look like the computer crashed % if rem(i,30) == 0 % fprintf('%2.2f, ',i/length(d)*100); % end
%this is the strcomp to find the mcia_8 files:
elseifstrcmpi(filename(1:5),'mcia3')
hasData=true;
try tbl = readtable([rootdir '1st Run\Incoming\' filename],'HeaderLines',1);
catch fprintf('\n%s is empty or otherwise problematic',filename);
%print out the file name so you can go see why the program
%didn't load it (probably it was empty).
hasData=false;
end [numOfObs,numOfDat]=size(tbl);
if numOfDat<2
hasData=false;
end
isTV=true;
cnt=0;
if ~isempty(tbl)&hasData
if ~isempty(tbl.Time)
%This loop looks for the end of file:
while isTV
%since there is a bunch of logger junk in the file after the last data point
cnt=cnt+1;
tmp=tbl.Time{cnt};
%Need to dump the data after this line ifstrcmp(tmp(1:5),'-----')
isTV=false;
cut=cnt-1;
else if cnt+1>length(tbl.Time)
isTV=false;
cut=cnt;
end end end
%cutoff everything that follows that '-----' line.
tbl=tbl(1:cut,:);
%get a time vector of the raw data
Time=datenum(tbl.Time,'mm/dd/yy HH:MM:SS.FFF');
tv_dv = datevec(Time+tvoffset_mci);
%rounding Date Vectors to fit into master time vector:
tv_rnd = datevec(datenum([tv_dv(:,1:5)round(tv_dv(:,6)/10).*10]));
%turn back into a single vector:
tv_rnd = datenum(tv_rnd);
%find where the data fits into the master chronology:
[c,ia,ib] = intersect(tv_rnd,tv);
%now put the data in it's place in the chronology:
CH4_mci(ib) = tbl.x_CH4__ppm(ia);
D13C_mci(ib) = tbl.d13C(ia);
H2O_mci(ib) = tbl.x_H2O__ppm(ia);
GasP_torr_mci(ib)= tbl.GasP_torr_se(ia);
RD0(ib) = tbl.RD0_us(ia);
% %if rem(i,30) == 0 %this is just a progress text, you can omit it if you want. % %fprintf('%2.2f, ',i/length(d)*100); % %end end end end end clear tv_dv tv_dv_rnd tv_rnd CO2_ppm D13C_VPDB_CO2 H2O_ppm foo filename c
clear ia ib tbl cnt isTV tmp cut
%do corrections on the D13C for CO2 from Elyn, year2016 / Ignore in 2018
%D13C_recalc = D13C-(0.000008.*CO2_i.^2-0.0307.*CO2_i-0.338);
%%load the autochamber data next, preallocate vectors:
rootdir = 'P:\1st Run\Incoming\';
d = dir(rootdir);
tv = [datenum(2018,6,1,0,0,0):1/8640:datenum(2018,10,1,0,0,0)]';
CH4 = NaN.*ones(length(tv),1);
CO2 = NaN.*ones(length(tv),1);
CalStg = NaN.*ones(length(tv),1);
ChNum = NaN.*ones(length(tv),1);
TC = NaN.*ones(length(tv),25);
fprintf('Autochamber:') for i = 3:length(d)
filename = d(i).name;
ifstrcmpi(filename(1:6),'Auto05')
auto05 = load([rootdir filename]);
Time = datenum(auto05(:,2),1,auto05(:,3),fix(auto05(:,4)./100),rem(auto05(:,4),100),auto05(:,5));
%generate datevec and apply time offset:
tv_dv = datevec(Time+tvoffset_a);
%rounding date vectors just like above:
tv_rnd = datevec(datenum([tv_dv(:,1:5)round(tv_dv(:,6)/10).*10]));
tv_rnd = datenum(tv_rnd);
[c,ia,ib] = intersect(tv_rnd,tv);
CO2(ib) = auto05(ia,18);
CH4(ib) = auto05(ia,44);
CalStg(ib) = auto05(ia,8);
ChNum(ib) = auto05(ia,6);
TC(ib,1:25) = auto05(ia,19:43);
% %if rem(i,30) == 0 % %fprintf('%2.2f, ',i/length(d)*100); % %end end
end clear tv_dv tv_rnd auto05 ia ib c i
fprintf('\n');
save('WorkingData.mat');
%if not loading from raw, just load the work files:
else [Tair,par,Kin,RH,Pb,ustar,WT,Tsoil,Tsoil_all,tv_mb]=LoadECData(false);
load('WorkingData.mat');
end
%perform CH4 calibrations:
calibrationch4logical;
%perform CO2 calibrations:
calibrationCO2logical;
%fprintf('debug');
%useful to have a var that points to the maximum length
%of the vectors so you don't refrence a point outside of it in the
%next loop:
maxVal=length(tv);
%a flag to know when to exit, a for loop is safer then a
%while loop when programing, but since the time the chamber is closed
%varies depending on the time of day, I found a it easier to use:
notEndOfData=true;
%This counts up one for each new flux:
num=0;
%index is needed since I am not using a for loop:
i=1;
while notEndOfData
period=0;
checkForChSwitch=true;
%Need to cycle through and find out when the
% channel switches, everything from this point to that point is a % single chamber measurement period:
while checkForChSwitch
if i+period>maxVal
checkForChSwitch=false;
else if ChNum(i)~=ChNum(i+period)
checkForChSwitch=false;
else
period=period+1;
end end end
endTime=i+period;
%skip calibration, when the channel number is empty or when there is insufficent data to calculate flux
%or when the CO2 data is far too low to be real:
if ~isnan(ChNum(i))&endTime-5>i&CalStg(i)==0&nanmin(CO2(i+2:endTime-2))>100
% Get the met variables to calculate the flux from
metPick = find(tv_mb >= tv(i));
%grab the closest met perameters which occur after measurment (within 30min)
metPick = metPick(1);
Tair_AC=Tair(metPick);
Pb_AC=Pb(metPick);
WtrLv_Met=WT(metPick);
%calculate the fluxes and variables requested:
num=num+1;
%Snip off the edges where flow could be contaminated from channel switching:
[CO2_full(num,1),CO2_fullR2(num,1)] = GetFlux(CO2(i+2:endTime-2,ChNum(i),false)),Tair_AC,Pb_AC;
[CH4_full(num,1),CH4_fullR2(num,1)] = GetFlux(CH4(i+2:endTime-2),ChNum(i),true),Tair_AC,Pb_AC;
[CO2d13Cchangefull(num,1),CO2d13CchangefullR2(num,1)]=GetD13C(d13C_CO2_recalc(i+2:endTime-2),ChNum(i)),Tair_AC,Pb_AC;
[CH4d13Cchangefull(num,1),CH4d13CchangefullR2(num,1)]=GetD13C(d13C_CH4_recalc(i+2:endTime-2),ChNum(i)),Tair_AC,Pb_AC;
%if there is sufficent data, calculate the last 90second flux, omitting the last 20 seconds (when the channel switches):
if period>=11&endTime<maxVal
[CO2_last90sec(num,1),CO2_last90secR2(num,1)]=GetFlux(CO2(endTime-11:endTime-2),ChNum(i),false),Tair_AC,Pb_AC;
[CH4_last90sec(num,1),CH4_last90secR2(num,1)]=GetFlux(CH4(endTime-11:endTime-2),ChNum(i),true),Tair_AC,Pb_AC;
%if there isn't enough data, just use the last 90 seconds (the R2 will reflect if its a problem) elseif endTime-9>1&endTime<maxVal
[CO2_last90sec(num,1),CO2_last90sec(num,1)]=GetFlux(CO2(endTime-9:endTime),ChNum(i),false),Tair_AC,Pb_AC;
[CH4_last90sec(num,1),CH4_last90secR2(num,1)]=GetFlux(CH4(endTime-9:endTime),ChNum(i),true),Tair_AC,Pb_AC;
%if their is less than 90 seconds recorded, then just store blanks:
else
CO2_last90sec(num,1)=nan;
CO2_last90secR2(num,1)=nan;
CH4_last90sec(num,1)=nan;
CH4_last90secR2(num,1)=nan;
end
CH4reci = CH4(i+2:endTime-2).^-1;
p(num,:) = regress(d13C_CH4_recalc(i+2:endTime-2),[ones(size(CH4reci)),CH4reci]);
keelingm(num,1) = p(num,2);
keelingb(num,1) = p(num,1);
%keelingbdiff = keelingb(num,1) - keelingb(num-1,1)
CO2_min(num,1)=nanmin(CO2(i+2:endTime-2));
CO2_max(num,1)=nanmax(CO2(i+2:endTime-2));
CO2_ave(num,1)=nanmean(CO2(i+2:endTime-2));
CH4_min(num,1)=nanmin(CH4(i+2:endTime-2));
CH4_max(num,1)=nanmax(CH4(i+2:endTime-2));
CH4_ave(num,1)=nanmean(CH4(i+2:endTime-2));
CO2_D13C_low(num,1)=nanmin(d13C_CO2_recalc(i+2:endTime-2));
CO2_D13C_mean(num,1)=nanmean(d13C_CO2_recalc(i+2:endTime-2));
CO2_D13C_high(num,1)=nanmax(d13C_CO2_recalc(i+2:endTime-2));
CH4_D13C_low(num,1)=nanmin(d13C_CH4_recalc(i+2:endTime-2));
CH4_D13C_mean(num,1)=nanmean(d13C_CH4_recalc(i+2:endTime-2));
timAvePeriod_sec(num,1)=endTime-i-4*10;
CH4_D13C_high(num,1)=nanmax(d13C_CH4_recalc(i+2:endTime-2));
GasP_torr(num,1)=nanmean(GasP_torr_mci(i+2:endTime-2));
RingDn(num,1)=nanmean(RD0(i+2:endTime-2));
Ch_Num(num,1)=ChNum(i);
TairC(num,1)=Tair_AC;
WtrLv(num,1)=WtrLv_Met;
TV(num,1)=tv(i);
TVend(num,1)=tv(endTime);
Cal_Stg(num,1)=CalStg(i);
end
%add the time monitoring the last chamber to the index:
i=i+period+1;
%check to see if that put us over the end of the file.
if i>maxVal
notEndOfData=false;
end end
%turn the vectors into a data table.
datT=table(CO2_full,CO2_fullR2,CO2_last90sec,CO2_last90secR2,CO2_min,CO2_max,CO2_max,CO2_ave,CO2_D13C_low,CO2_D13C_high,...
CH4_full,CH4_fullR2,CH4_last90sec,CH4_last90secR2,CH4_max,CH4_min,CH4_ave,CH4_D13C_low,CH4_D13C_high,...
GasP_torr,RingDn,TairC,WtrLv,TV,TVend,timAvePeriod_sec,Ch_Num,CO2d13Cchangefull,CO2d13CchangefullR2,CO2_D13C_mean,...
CH4d13Cchangefull,CH4d13CchangefullR2,CH4_D13C_mean,keelingm, keelingb);
function[fl,r2]=GetFlux(conDat,Tair_AC,Pb_AC,ColNum,isMeth)
%calc flux (mean[]fornow);
%Line82 to 85 of working_MB_chamber_Fch4
collar = [14105238-11391213139]./100;
UBC_biomet_constants;
%chamber area from Lai et al. 2012 Biogeosciences
A = 0.21;
h = 0.215; a = 0.52/2;
%volume of spherical cap with height of cap, h = 0.215 m and radius of base of cap, a = 0.52/2 m from Lai et al. 2012 Biogeosci.
cap = pi.*h./6.*(3.*a.^2+h.^2);
%chamber volume
V = A.*collar(ColNum) + cap;
timeStep=10:10:1000;
%Flux = change_in_gas(Dg/Dt)*Pressure(Pa)*VolumeOfChamber(m3)/(SurfaceArea(m2)*Temp(K)*R)
press=Pb_AC*1000;
[dg_dt,r2] = polyfit1(timeStep(1:length(conDat))',conDat,1);%get the slope of the gas change.
fl=dg_dt(1)*press*V/(A*(Tair_AC+ZeroK)*R); %nmol m-2 s-1;
if isMeth
fl*1000;%convert if it's methane data
end
function[change,r2]=GetD13C(conDat,Tair_AC,Pb_AC,ColNum)
%this is the "GetFlux" function by Graham just copy pasted and minor
%changes in the structure to get the slope and R2 of the D13C change
%there are a few useless variables (for this function) in it.
collar = [14105238-11391213139]./100;
UBC_biomet_constants;
A = 0.21; %chamber area from Lai et al. 2012 Biogeosciences
h = 0.215; a = 0.52/2;
cap = pi.*h./6.*(3.*a.^2+h.^2); %volume of spherical cap with height of cap, h = 0.215 m and radius of base of cap, a = 0.52/2 m from Lai et al. 2012 Biogeosci.
V = A.*collar(ColNum) + cap;%chamber volume
timeStep=10:10:1000;
%Flux = change_in_gas(Dg/Dt)*Pressure(Pa)*VolumeOfChamber(m3)/(SurfaceArea(m2)*Temp(K)*R)
press=Pb_AC*1000;
[dg_dt,r2] = polyfit1(timeStep(1:length(conDat))',conDat,1);%get the slope of the gas change.
change=dg_dt(1); %nmol m-2 s-1;
function UnzipData(rootdir)
%This is an unzipping function. Minor changes from Elyn's code to deal
% with changes in file structure.
d=rootdir;
d = dir([rootdir]);
for i = 3:length(d) iflength(d(i).name)>4&~strcmp(d(i).name(1:3),'Inc')&~strcmp(d(i).name(1:3),'arc') ifstrcmp(d(i).name(1:4),'2017')
filename = dir([rootdir '\' d(i).name]);
for ii=3:length(filename)
currentFile=[rootdir '\' d(i).name, '\' filename(ii).name];
fprintf('\n%d w %d:%s',i,ii,currentFile);
ifstrcmp(currentFile(end-2:end),'zip')
unzip(currentFile,[rootdir '\Incoming\']);
fprintf(' File Unzipped');
elseifstrcmp(currentFile(end-2:end),'txt') copyfile(currentFile, [rootdir '\Incoming\']);
fprintf(' File Copied');
else fprintf(' Nothing done with this file');
end end else
filename = dir([rootdir '\' d(i).name]);
%fprintf('\n%d:%s',i,[rootdir '\' d(i).name, '\' filename(4).name]);
fnms= dir([rootdir '\' d(i).name '\' filename(4).name]);
for ii=3:length(fnms) fprintf('\n%d w %d:%s',i,ii,[rootdir '\' d(i).name, '\' filename(4).name '\' fnms(ii).name]);
currentFile=[rootdir '\' d(i).name '\' filename(4).name '\' fnms(ii).name];
ifstrcmp(currentFile(end-2:end),'zip')
unzip(currentFile,[rootdir '\Incoming\']);
fprintf(' File Unzipped');
elseifstrcmp(currentFile(end-2:end),'txt') copyfile(currentFile, [rootdir '\Incoming\']);
fprintf(' File Copied');
else fprintf(' Nothing done with this file');
end end end end
end % for i = 3:58 % filename = dir([rootdir '\' d(i).name]); % %fprintf('\n%d:%s',i,[rootdir '\' d(i).name, '\' filename(4).name]); % fnms= dir([rootdir '\' d(i).name '\' filename(4).name]); % for ii=3:length(fnms) % fprintf('\n%d w %d:%s',i,ii,[rootdir '\' d(i).name, '\' filename(4).name '\' fnms(ii).name]); % currentFile=[rootdir '\' d(i).name '\' filename(4).name '\' fnms(ii).name]; % if strcmp(currentFile(end-2:end),'zip') % unzip(currentFile,[rootdir '\Incoming\']); % fprintf(' File Unzipped'); % elseif strcmp(currentFile(end-2:end),'txt') % copyfile(currentFile, [rootdir '\Incoming\']); % fprintf(' File Copied'); % else % fprintf(' Nothing done with this file'); % end % end % end
%
% for i = 59:166 % filename = dir([rootdir '\' d(i).name]); % for ii=3:length(filename) % currentFile=[rootdir '\' d(i).name, '\' filename(ii).name]; % fprintf('\n%d w %d:%s',i,ii,currentFile); % if strcmp(currentFile(end-2:end),'zip') % unzip(currentFile,[rootdir '\Incoming\']); % fprintf(' File Unzipped'); % elseif strcmp(currentFile(end-2:end),'txt') % copyfile(currentFile, [rootdir '\Incoming\']); % fprintf(' File Copied'); % else % fprintf(' Nothing done with this file'); % end % end % end
%
% for i = 167:247 % filename = dir([rootdir '\' d(i).name]); % %fprintf('\n%d:%s',i,[rootdir '\' d(i).name, '\' filename(4).name]); % fnms= dir([rootdir '\' d(i).name '\' filename(4).name]); % for ii=3:length(fnms) % fprintf('\n%d w %d:%s',i,ii,[rootdir '\' d(i).name, '\' filename(4).name '\' fnms(ii).name]); % currentFile=[rootdir '\' d(i).name '\' filename(4).name '\' fnms(ii).name]; % if strcmp(currentFile(end-2:end),'zip') % unzip(currentFile,[rootdir '\Incoming\']); % fprintf(' File Unzipped'); % elseif strcmp(currentFile(end-2:end),'txt') % copyfile(currentFile, [rootdir '\Incoming\']); % fprintf(' File Copied'); % else % fprintf(' Nothing done with this file'); % end % end % end fprintf('\n');
gibt ein leeres Array zurück, weil keine solche Einträge gefunden wurden. In dem Fall führt die folgende Anweisung zu einem Fehler: wenn es keine Einträge gibt, dann auch keinen ersten.
Ich könnte mir eine Fallunterscheidung vorstellen:
Code:
if ~isempty(metPick)
metPick = metPick(1);
% ... else % Was soll passieren, wenn nichts gefunden wurde? end
1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
gibt ein leeres Array zurück, weil keine solche Einträge gefunden wurden. In dem Fall führt die folgende Anweisung zu einem Fehler: wenn es keine Einträge gibt, dann auch keinen ersten.
Ahh, ich denke, ich sehe das "Problem". Du hast ziemlich sicher Recht und das Array bleibt leer. Ich habe mir die beiden Vektoren tv_mb und tv(i) noch einmal genauer angeschaut und tatsächlich wird die >= Bedingung nicht erfüllt. Jetzt muss ich nur herausfinden, warum in beiden Zeit-Vektoren unterschiedliche Zeitpunkte stehen, wenn sie doch vom gleichen Messzeitraum sind. Hmm, ...danke dir soweit schonmal herzlich!
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS