Skip to content

Commit

Permalink
Merge pull request #11 from cwingard/poc2w
Browse files Browse the repository at this point in the history 
PCO2W Corrections
  • Loading branch information
@danmergens
danmergens authored Sep 12, 2018
2 parents 6b83c06 + 50b5ace commit 4b3d531
Show file tree
Hide file tree
Showing 3 changed files with 194 additions and 63 deletions.
67 changes: 38 additions & 29 deletions ion_functions/data/co2_functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@

import numpy as np
import numexpr as ne
import scipy as sp
from ion_functions.utils import fill_value


Expand Down Expand Up @@ -97,6 +96,8 @@ def pco2_blank(raw_blank):
2014-02-19: Christopher Wingard. Updated comments.
2014-02-28: Christopher Wingard. Updated to except raw blank values
from a sparse array.
2018-03-04: Christopher Wingard. Updated to correctly calculate the
blank based on new code from the vendor.
Usage:
Expand All @@ -115,9 +116,7 @@ def pco2_blank(raw_blank):
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
#blank = -1. * sp.log10(raw_blank / 16384.)
blank = -1. * sp.log10(raw_blank)

blank = raw_blank / 16384.
return blank


Expand Down Expand Up @@ -207,11 +206,11 @@ def pco2_pco2wat(mtype, light, therm, ea434, eb434, ea620, eb620,
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
# reset inputs to arrays
### measurements
# measurements
mtype = np.atleast_1d(mtype)
light = np.atleast_2d(light)
therm = np.atleast_1d(therm)
### calibration coefficients
# calibration coefficients
ea434 = np.atleast_1d(ea434)
eb434 = np.atleast_1d(eb434)
ea620 = np.atleast_1d(ea620)
Expand All @@ -220,7 +219,7 @@ def pco2_pco2wat(mtype, light, therm, ea434, eb434, ea620, eb620,
cala = np.atleast_1d(cala)
calb = np.atleast_1d(calb)
calc = np.atleast_1d(calc)
### blank measurements
# blank measurements
a434blank = np.atleast_1d(a434blank)
a620blank = np.atleast_1d(a620blank)

Expand Down Expand Up @@ -251,6 +250,11 @@ def pco2_calc_pco2(light, therm, ea434, eb434, ea620, eb620,
2013-04-20: Christopher Wingard. Initial python code.
2014-02-19: Christopher Wingard. Updated comments.
2014-03-19: Christopher Wingard. Optimized.
2018-03-04: Christopher Wingard. Updated to correctly calculate pCO2 using
newly formulated code provided by the vendor. Original vendor code
incorrectly calculated the blank correction. Applies additional
corrections to calculations to avoid errors thrown when running a
blank measurement.
Usage:
Expand Down Expand Up @@ -281,43 +285,48 @@ def pco2_calc_pco2(light, therm, ea434, eb434, ea620, eb620,
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
# set constants
ea434 = ea434 - 29.3 * calt
eb620 = eb620 - 70.6 * calt
e1 = ea620 / ea434
e2 = eb620 / ea434
e3 = eb434 / ea434
# set constants -- original vendor formulation, reset below
# ea434 = ea434 - 29.3 * calt
# eb620 = eb620 - 70.6 * calt
# e1 = ea620 / ea434
# e2 = eb620 / ea434
# e3 = eb434 / ea434

# set the e constants, values provided by Sunburst
e1 = 0.0043
e2 = 2.136
e3 = 0.2105

# Extract variables from light array
light = light.astype(np.float)
#DRef1 = light[0] # Dark Reference LED
#DSig1 = light[1] # Dark Signal LED
#R434 = light[2] # 434nm Reference LED intensity
#S434 = light[3] # 434nm Signal Signal LED intensity
#R620 = light[4] # 620nm Reference LED intensity
#S620 = light[5] # 434nm Signal Signal LED intensity
Ratio434 = light[:, 6] # 434nm Ratio
Ratio620 = light[:, 7] # 620nm Ratio

# Convert thermistor counts to degrees C
therm = pco2_thermistor(therm)
# Convert blank light readings to absorbance
a434blank = pco2_blank(a434blank)
a620blank = pco2_blank(a620blank)

# calculate absorbance ratio, correcting for blanks
A434 = -1. * sp.log10(Ratio434 / a434blank) # 434 absorbance
A620 = -1. * sp.log10(Ratio620 / a620blank) # 620 absorbance
Ratio = A620 / A434 # Absorbance ratio
# correct the absorbance ratios using the blanks
AR434 = (Ratio434 / a434blank)
AR620 = (Ratio620 / a620blank)

# map out blank measurements and spoof the ratios to avoid throwing an error
m = np.where(AR434 == AR620)[0]
AR434[m] = 0.99999
AR620[m] = 0.99999

# Calculate the final absorbance ratio
A434 = -1 * np.log10(AR434) # 434 absorbance
A620 = -1 * np.log10(AR620) # 620 absorbance
Ratio = A620 / A434 # Absorbance ratio

# calculate pCO2
V1 = Ratio - e1
V2 = e2 - e3 * Ratio
RCO21 = -1. * sp.log10(V1 / V2)
RCO22 = (therm - calt) * 0.007 + RCO21
RCO21 = -1 * np.log10(V1 / V2)
RCO22 = (therm - calt) * 0.008 + RCO21
Tcoeff = 0.0075778 - 0.0012389 * RCO22 - 0.00048757 * RCO22**2
Tcor_RCO2 = RCO21 + Tcoeff * (therm - calt)
pco2 = 10.**((-1. * calb + (calb**2 - (4. * cala * (calc - Tcor_RCO2)))**0.5) / (2. * cala))
pco2[m] = fill_value # reset the blanks captured earlier to a fill value

return np.real(pco2)

Expand Down
124 changes: 124 additions & 0 deletions ion_functions/data/matlab_scripts/pco2wat/SAMI_pCO2_equations_Sunburst.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
% This m-file reads raw data from Sunburst SAMI pCO2 data file and
% calculates pCO2.
% Code supplied by Reggie Spaulding at Sunburst Sensors in Feb 2018.
%
% Test data file was created from data collected during the deployment of
% sensor SN C0123 in July 2017.
%
% *************************************************************************
clear all

%
% Sensor Specific Calibration coefficients for C0123, 2017
%
CalT=4.6539;
CalA=0.0422;
CalB=0.6761;
CalC=-1.5798;

% e's for calculating pCO2, constants provided by Sunburst based on lab
% determinations with BTB
e1=0.0043;
e2=2.136;
e3=0.2105;

File = 'testFile.txt';
fid = fopen(File);

% Read in SAMI hex data
i=1; j=1; k=1;
while 1
s=fgetl(fid);
if s == -1, break,end % indicates EOF
if ~isempty(s)
if (strcmpi('*',s(1:1))==1 && strcmpi('04',s(6:7))| strcmpi('05',s(6:7)))
if (length(s) == 81)
s = s(2:length(s));
while length(s) < 80
s = [s,fgetl(fid)];
end
AA(i,:)= s;
i=i+1;
end
else
end
end
end
fclose(fid);
[d1,d2]=size(AA);
% *************************************************************************
% Extract data from hex string
for i=1:d1
Type(i)=hex2dec(AA(i,5:6)); % Type 4 - Measurement, Type 5 - Blank
Time(i)=hex2dec(AA(i,7:14)); % Time
t_SAMI(i) = ((Time(i)/86400) +(365.25*4)+(1));
DarkRef(i)=hex2dec(AA(i,15:18)); % Dark Reference LED
DarkSig(i)=hex2dec(AA(i,19:22)); % Dark Signal LED
Ref434(i)=hex2dec(AA(i,23:26)); % 434nm Reference LED intensity
Sig434(i)=hex2dec(AA(i,27:30)); % 434nm Signal Signal LED intensity
Ref620(i)=hex2dec(AA(i,31:34)); % 620nm Reference LED intensity
Sig620(i)=hex2dec(AA(i,35:38)); % 434nm Signal Signal LED intensity
Ratio434(i)=hex2dec(AA(i,39:42)); % 434nm Ratio
Ratio620(i)=hex2dec(AA(i,43:46)); % 620nm Ratio
Battery(i)=hex2dec(AA(i,71:74)); % Raw Battery Value
BattV(i) = Battery(i)*15/4096; % Battery voltage
Therm(i)=hex2dec(AA(i,75:78)); % Thermistor raw value
end

k=1;
for i=1:d1
BFlag(i)=nan;
if Type(i)==5 % Blank measurement
A434Bka(k)=Ratio434(i)/16384;
A620Bka(k)=Ratio620(i)/16384;
blanktime(k)=((Time(i))/(60*60*24)); % Time of blank
datetime1(i)=(Time(i)/(60*60*24)); % Time total
Tc_SAMI(i)=(1./((0.0010183)+(0.000241.*log((Therm(i)./...
(4096-Therm(i))).*17400))+((0.00000015.*log((Therm(i)./...
(4096-Therm(i))).*17400).^3))))-273.15; % Temp
k=k+1;
end
end

% PCO2 Calculations
i2=1;
for i=1:d1
if Type(i)==5
A434Bk=A434Bka(i2);
A620Bk=A620Bka(i2);
calcCO2_1(i)=0;
i2=i2+1;
else if Type(i)==4
k434(i)=A434Bk;
k620(i)=A620Bk;
A434(:,i)=-log10((Ratio434(i)./16384)/k434(i));
A620(:,i)=-log10((Ratio620(i)./16384)/k620(i));
Ratio(:,i)=(A620(:,i))/(A434(:,i));
datetime1(i)=(Time(i)./(60*60*24));
Tc_SAMI(i)=(1./((0.0010183)+(0.000241.*log((Therm(i)./...
(4096-Therm(i))).*17400))+((0.00000015.*log((Therm(i)./...
(4096-Therm(i))).*17400).^3))))-273.15; % Temp
RCO2(:,i) = -log10((Ratio(:,i)-e1)./(e2-e3.*Ratio(:,i)));
Tcor_RCO2i(:,i) = RCO2(:,i)+(0.008.*(Tc_SAMI(i)-CalT));
Tcoeff(:,i) = (0.0075778)-(0.0012389.*Tcor_RCO2i(:,i))-...
(0.00048757.*Tcor_RCO2i(:,i).^2);
Tcor_RCO2(:,i) = RCO2(:,i)+(Tcoeff(:,i)).*(Tc_SAMI(i)-CalT);
calcCO2_1(:,i) = (10.^(((-1*CalB)+((CalB^2)-(4*CalA.*...
(CalC-Tcor_RCO2(:,i)))).^0.5)./(2*CalA)));
end
end
end

% *************************** Write output file ***************************
outfile_CO2=strcat(File(1:end-4),'_out.xls');
col_header={'Time','pCO2','Temp','BattV','DarkRef','DarkSig',...
'Ratio434','Ref434','Sig434','Ratio620','Ref620','Sig620'};
datafinal_CO2=[t_SAMI;calcCO2_1;Tc_SAMI;BattV;DarkRef;DarkSig...
;Ratio434;Ref434;Sig434;Ratio620;Ref620;Sig620];
fid = fopen(outfile_CO2,'w');
fprintf(fid,['Time\t pCO2\t Temp\t BattV\t DarkRef\t DarkSig\t'...
'Ratio434\t Ref434\t Sig434\t Ratio620\t Ref620\t Sig620\r\n']);
fmt = ['%8.4f \t %8.4f \t %8.4f \t %8.2f \t %8f \t %8f \t %8f \t'...
'%8f \t %8f \t %8f \t %8f \t %8f \r\n'];
fprintf(fid,fmt,datafinal_CO2);
fclose(fid);
66 changes: 32 additions & 34 deletions ion_functions/data/test/test_co2_functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,42 +21,41 @@ class Testpco2FunctionsUnit(BaseUnitTestCase):
def setUp(self):
###### Test data for PCO2W ######
raw_strings = np.array([
'*7E2705CBACEE7F007D007D0B2A00BF080500E00187034A008200790B2D00BE080600DE0C1406C98C',
'*7E2704CBACEECB008000880B2900B2080600D300FB0263007F00890B2B00B4080700CE0C5106C884',
'*7E2704CBACEF43007E00890B27014408070189045B0875007E00870B2B0140080201860C5506C601',
'*7E2704CBACEFBB007E00820B2A042B0803051F16182785008000850B2A043C080405390C5506C5A9',
'*7E2704CBACF033007F00840B28054D080606831CCC330A007F00850B290551080406800C5606C556',
'*7E2704CBACF0AB007E00770B2804DE080605F31A672E9F008100790B2F04E0080705F70C5606C5EB',
'*7E2704CBACF1230081007C0B2B00BF080800DB01BF0390007C00790B3000C7080B00EA0C5606C80A',
'*7E2704CBACF19B008000750B2B01D20807023008310E94007E00730B2A01D2080502290C5706C0A2',
'*7E2704CBACF213008000740B2A01D50808042F08501FC6007D00780B3201D8080A04350C5706C0A5',
'*7E2704CBACF28B007F00710B2F0174080203570615189B008100730B2A0174080A03580C5706C125',
'*7E2704CBACF303007E006B0B2C019B080803CC07001CBA007F006F0B300199080803D00C5706C4E3'
'*BC2705D5A7C0E10082005A0CA9090E07CB08E82DCA4B1C0082005A0CA9090E07CD08EC0C3208C38A',
'*BC2704D5A7E2B1007E005A0CB1022F07C40443099F226D007F005A0CAF022F07C404400C3F08BE2E',
'*BC2704D5A7FEC90080005A0CAD028707CC03160B711800008300580CAC028607CC03160C4007389C',
'*BC2704D5A8006F0083005B0CA1028D07DE02F70BA016AF0081005A0CA2028D07E202F70C4007A16D',
'*BC2704D5A802150080005A0C98028D07E902DE0BAE15BF0081005A0C98028E07EB02DF0C40080C6B',
'*BC2704D5A803BB007F00590C98028307EB02EA0B6B161B008100580C94028107EE02EB0C41085372',
'*BC2704D5A80560007F005A0C9A027407E403050B26171F0083005C0C99027607E903060C4008862B',
'*BC2704D5A80707007F005A0CA0027007DE03210AF9182A008000590CA0026D07DE03240C400895E4',
'*BC2704D5A808AD0082005B0CA2026607D203470AC019A00080005C0CA6026607D403490C3F0899FF',
'*BC2704D5A80A54007F00560CAB025407CC03860A701BCF0083005C0CAA025707D003840C3F089BE2',
'*BC2704D5A80BF90080005A0CB3024907C603B60A2D1D99008100580CAF024707C603B90C3F089C58',
'*BC2704D5A80E140080005B0CB3023807C0041009CA20C40082005A0CB3023807C004110C3F08A0D4',
'*BC2704D5A8102F007F00580CB9022D07BA04350999222D0080005A0CB4022E07BC04370C3E08B067',
'*BC2704D5A812E70081005B0CBE022107B00479094A24AD0080005A0CBC022007B8047B0C3E08CEE4',
])

# reagent constants (instrument and reagent bag specific)
self.ea434 = np.ones(11) * 19706.
self.ea620 = np.ones(11) * 34.
self.eb434 = np.ones(11) * 3073.
self.eb620 = np.ones(11) * 44327.
self.calt = np.ones(11) * 16.5
self.cala = np.ones(11) * 0.0459
self.calb = np.ones(11) * 0.6257
self.calc = np.ones(11) * -1.5406
self.calt = np.ones(14) * 4.6539
self.cala = np.ones(14) * 0.0422
self.calb = np.ones(14) * 0.6761
self.calc = np.ones(14) * -1.5798

# expected outputs
self.therm = np.array([18.8526, 18.8765, 18.9245, 18.9485,
18.9485, 18.9485, 18.8765, 19.0686,
19.0686, 19.0446, 18.9725])
self.pco2 = np.array([fill_value, 294.1720, 311.3361, 319.0101,
319.8925, 319.8950, 305.8104, 317.9661,
284.3676, 280.2324, 280.0354])
self.therm = np.array([7.3151, 7.4258, 16.2306, 13.8108, 11.3900,
9.8019, 8.6674, 8.3345, 8.2458, 8.2014,
8.1792, 8.0905, 7.7359, 7.0716])
self.pco2 = np.array([fill_value, 609.8626, 394.3221, 351.6737, 321.4986,
324.0670, 339.4317, 358.3335, 388.1735, 436.8431,
481.1713, 566.8172, 607.5355, 685.8555])

# parse the data strings
self.mtype = np.zeros(11, dtype=np.int)
self.light = np.zeros((11, 14), dtype=np.int)
self.traw = np.zeros(11, dtype=np.int)
for i in range(11):
self.mtype = np.zeros(14)
self.light = np.zeros((14, 14))
self.traw = np.zeros(14)
for i in range(14):
# parse the raw strings into subelements, such as the driver would
# provide.
s = raw_strings[i]
Expand All @@ -72,8 +71,8 @@ def setUp(self):
a434blnk = self.light[i, 6]
a620blnk = self.light[i, 7]

self.a434blnk = np.ones(11) * a434blnk
self.a620blnk = np.ones(11) * a620blnk
self.a434blnk = np.ones(14) * a434blnk
self.a620blnk = np.ones(14) * a620blnk

def test_pco2_pco2wat(self):
"""
Expand All @@ -99,7 +98,7 @@ def test_pco2_pco2wat(self):
### bulk case ###
tout = co2func.pco2_thermistor(self.traw)
pco2out = co2func.pco2_pco2wat(self.mtype, self.light, self.traw,
self.ea434, self.eb434, self.ea620, self.eb620,
fill_value, fill_value, fill_value, fill_value,
self.calt, self.cala, self.calb, self.calc,
self.a434blnk, self.a620blnk)

Expand All @@ -111,8 +110,7 @@ def test_pco2_pco2wat(self):
for mtype in self.mtype:
tout = co2func.pco2_thermistor(self.traw[indx])
pco2out = co2func.pco2_pco2wat(mtype, self.light[indx, :], self.traw[indx],
self.ea434[indx], self.eb434[indx],
self.ea620[indx], self.eb620[indx],
fill_value, fill_value, fill_value, fill_value,
self.calt[indx], self.cala[indx],
self.calb[indx], self.calc[indx],
self.a434blnk[indx], self.a620blnk[indx])
Expand Down

0 comments on commit 4b3d531

Please sign in to comment.