R/Sensors_Processing.R

Defines functions rain_PolySPL5 wind_PolySPL8 Process_RawIMU IMU_processAcc IMU_processHeading Process_Ramses ra_single Process_pH_SBE Process_DO_AADI_SVU Process_DO_Bittig

#**************************************************

# Data Processing

#**************************************************

#**************************************************
#*
# Compute do
#*
#**************************************************

# Inputs : C1phase,C2phase,temp, Pres : data from the Optode. tempCTD,salCTD, PRESCTD : Data from the CTD

Process_DO_Bittig <- function(C1phase,C2phase,temp, Pres,tempCTD,salCTD, PRESCTD,COEF=NULL) {
  
  #tempCTD <- approx(PRESCTD, tempCTD, Pres, rule=2)$y #Approx tempCTD on DO pressure
  tempCTD <- temp #We use the temperature from the Optode
  
  #Approx tempCTD on DO pressure
  salCTD <- approx(PRESCTD,salCTD, Pres, rule=2,ties = "mean")$y 
  
  
  
  #COEF if NULL (for plotting)
  if (is.null(COEF)){
    COEF <- c(5.6725661e-03,8.2915275e-05,1.0033795e-06,6.2236942e-02,-9.3470722e-05,-1.4554620e-02,1.2110645e-03) # From Henry
  }
  
  
  TCPhase <-  (C1phase - C2phase) 
  
  KSV <-  COEF[1] + (COEF[2]  * temp) + (COEF[3] * temp * temp)
  
  #print(KSV[20:100])
  
  P0 <- COEF[4]  + COEF[5] * temp
  
  #print(P0[20:100])
  
  PC <- COEF[6]  + COEF[7] * TCPhase
  
  #print(PC[20:100])
  pO2 <- ((P0/PC)-1) / KSV
  
  
  rhumid=1
  atm_press=1.01325
  
  pH2O = rhumid * (exp(24.4543-(67.4509*(100/(tempCTD+273.15)))-(4.8489*log(((273.15+tempCTD)/100)))-0.000544*salCTD))
  
  th0=1-(0.999025+0.00001426*tempCTD-0.00000006436*tempCTD^2)
  sca_T=log((298.15-tempCTD)/(273.15+tempCTD))
  
  
  oxy_sol=((exp(2.00856+3.224*sca_T+3.99063*sca_T^2+4.80299*sca_T^3+0.978188*sca_T^4+1.71069*sca_T^5+salCTD*(-0.00624097-0.00693498*sca_T-0.00690358*sca_T^2-0.00429155*sca_T^3)-0.00000031168*salCTD^2))/0.022391903)
  oxy_sol1 <- oxy_sol 
  
  #oxy_sol=oxy_sol*P_atm*(((1-pH2O/P_atm)*(1-th0*P_atm))/((1-pH2O)*(1-th0)));
  # pressure corrected O2 solubility / umol atm / l
  oxy_sol  = oxy_sol *  atm_press*(((1-pH2O/atm_press)*(1-th0*atm_press))/((1-pH2O)*(1-th0)))
  oxy_sol2 <- oxy_sol 
  
  # oxygen concentration in umol/l (salinity corrected)
  O2conc_sal=(pO2*oxy_sol)/((atm_press-pH2O)*0.20946*1013.25)
  O2conc_sal3 <- O2conc_sal 
  
  # correcion de pression
  O2conc_sal = O2conc_sal * (1+ ((3.2)/100 *(Pres/1000)))
  
  return(O2conc_sal)
  
}

#**************************************************
#*
# Compute do
#*
#**************************************************

# Inputs : C1phase,C2phase,temp, Pres : data from the Optode. tempCTD,salCTD, PRESCTD : Data from the CTD

Process_DO_AADI_SVU <- function(C1phase,C2phase,temp, Pres,tempCTD,salCTD, PRESCTD,COEF=NULL, PHASECOEF0=NULL) {
  
  tempCTD <- approx(PRESCTD, tempCTD, Pres, rule=2, ties = "mean")$y #Approx tempCTD on DO pressure
  #tempCTD <- temp #We use the temperature from the Optode
  
  #Approx tempCTD on DO pressure
  salCTD <- approx(PRESCTD,salCTD, Pres, rule=2,ties = "mean")$y 
  
  #COEF if NULL (for plotting)
  if (is.null(COEF)){
    COEF <- c(5.6725661e-03,8.2915275e-05,1.0033795e-06,6.2236942e-02,-9.3470722e-05,-1.4554620e-02,1.2110645e-03) # From Henry
  }
  if (is.null(PHASECOEF0)){
    PHASECOEF0 <- 0 # From Henry
  }
  
  TCPhase <-  (C1phase - C2phase) + PHASECOEF0
  
  # first part of pressure effect
  TCPhase = TCPhase + 0.1 * Pres/1000 # do O2-independent phase adjustment
  
  KSV <-  COEF[1] + (COEF[2]  * temp) + (COEF[3] * temp * temp)
  
  #print(KSV[20:100])
  
  P0 <- COEF[4]  + COEF[5] * temp
  
  #print(P0[20:100])
  
  PC <- COEF[6]  + COEF[7] * TCPhase
  
  #print(PC[20:100])
  # oxygen concentration in umol/l (in freshwater)
  O2_molar_fresh <- ((P0/PC)-1) / KSV
  
  rhumid = 1
  atm_press=1.01325
  
  pH2Ofresh = rhumid * (exp(24.4543-(67.4509*(100/(tempCTD+273.15)))-(4.8489*log(((273.15+tempCTD)/100)))-0.000544*0))
  pH2O = rhumid * (exp(24.4543-(67.4509*(100/(tempCTD+273.15)))-(4.8489*log(((273.15+tempCTD)/100)))-0.000544*salCTD))
  
  #th0=1-(0.999025+0.00001426*tempCTD-0.00000006436*tempCTD^2)
  sca_T=log((298.15-tempCTD)/(273.15+tempCTD))
  Scorr=((exp(salCTD*(-0.00624097-0.00693498*sca_T-0.00690358*sca_T^2-0.00429155*sca_T^3)-0.00000031168*salCTD^2)))
  # oxygen concentration in umol/l (salinity corrected)
  O2_molar_salt=O2_molar_fresh*(atm_press-pH2Ofresh)/(atm_press-pH2O)*Scorr

  # correcion de pression
  #O2conc_sal = O2conc_sal * (1+ ((3.2)/100 *(Pres/1000)))
  pcfactor=4.19+0.022*tempCTD
  O2conc_sal = O2_molar_salt * (1+ ((pcfactor)/100 *(Pres/1000))) # in umol/L, pressure corrected, salinity compensated
  
  # divide by density to get umol/kg for DOXY
  
  return(O2conc_sal)
  
}




#**************************************************
#*
# Compute pH
#*
#**************************************************

# Inputs : 

Process_pH_SBE<-function(data,NumberPhase="ASC",k0=-1.392151,k2=-1.0798E-03,coefsp=c(2.5064E-05,-4.4107E-08,4.7311E-11,-2.8822E-14,9.2132E-18,-1.1965E-21)){
  phtot=NULL
  
  dataCTD<-data$sbe41[data$sbe41$PhaseName==NumberPhase,]
  datapH<-data$sbeph[data$sbeph$PhaseName==NumberPhase,]
  
  if ((length(datapH$pH_mV)>5) & (length(dataCTD$Temperature_degC)>5)){
    #il y a assez de data
    
    ## choix method pour approxfun
    if (length(unique(dataCTD$Pressure_dbar))>2){
      #il y a au moins 2 points differents
      approx_method="linear"
    } else {
      approx_method = "constant"
    }
    
    #interpolation des donnees CTD
    t<-approxfun(dataCTD$Pressure_dbar,dataCTD$Temperature_degC,rule = 2,ties="mean",method=approx_method)
    S<-approxfun(dataCTD$Pressure_dbar,dataCTD$Salinity_PSU,rule = 2,ties="mean",method=approx_method)
    
    #calcul
    Press<-datapH$Pressure_dbar
    Vrs<-datapH$pH_mV
    Temp<-t(Press)
    Tk<-273.15+Temp #degrees Kelvin
    Salt<-S(Press)
    
    #P<-10*p #passage en bar
    
    # ************************************************************************
    #  SET SOME CONSTANTS
    # ************************************************************************
    #Universal gas constant, (R) , http://physics.nist.gov/cgi-bin/cuu/Value?r
    R    = 8.31446 # J/(mol K) 
    Fa    = 96485  #Faraday constant Coulomb / mol
    
    ln10 = log(10) #natural log of 10
    
    # ************************************************************************
    #  CALCULATE PHYSICAL AND THERMODYNAMIC DATA
    # Dickson, A. G., Sabine, C. L., & Christian, J. R. (2007). Guide to best
    # practices for ocean CO2 measurements.
    # ************************************************************************
    
    # IONIC STRENGTH OF SEAWATER (mol / kg H2O)
    # Varified units by comparing to Dickson et al. 2007: Chap 5, p10 Table 2
    # Dickson et al. 2007: Chap 5, p13 Eq 34
    IonS = 19.924 * Salt / (1000 - 1.005 * Salt)
    
    # MEAN SEAWATER SULFATE CONCENTRATION (mol / kg solution)
    # This wants to be mol/kg sw  as KHSO4 is on that scale
    # Dickson et al. 2007: Chap 5, p10 Table 2
    Stotal = (0.14 / 96.062) * (Salt / 1.80655)
    
    # MEAN SEAWATER CHLORIDE CONCENTRATION  (mol / kg H20)
    # this wants to be mol/kg H2O as activity is on mol/kg H2O scale
    # Dickson et al. 2007: Chap 5, p10 Table 2
    Cltotal = 0.99889 / 35.453 * Salt / 1.80655 # %(mol / kg solution)
    Cltotal = Cltotal /(1 - 0.001005 * Salt) #  % (mol / kg H20)
    
    # BISULFIDE DISSCIATION CONSTANT AT T,S AND IONIC STRENGTH(mol/kg solution)
    # Dickson et al. 2007: Chap 5, p12 Eq 33
    Khso4 = exp(-4276.1 / Tk + 141.328 - 23.093 * log(Tk) + 
                  (-13856 / Tk + 324.57 - 47.986 * log(Tk)) * IonS^ 0.5 + 
                  (35474 / Tk - 771.54 + 114.723 * log(Tk)) * IonS - 
                  2698 / Tk * IonS^ 1.5 + 1776 / Tk * IonS ^ 2 + 
                  log(1 - 0.001005 * Salt))
    
    # Millero 1983 Chemical Oceanography vol 8
    # partial molar volume and compressibility of HSO4 in seawater. 
    deltaVHSO4 = -18.03 + 0.0466 * Temp + 0.000316 * Temp^ 2
    KappaHSO4 = (-4.53 + 0.09 * Temp) / 1000
    
    #  Press changed from dbar to bar here by / 10
    lnKhso4fac = (-deltaVHSO4 + 0.5 * KappaHSO4 * (Press / 10)) * (Press / 10) / (R * 10 * Tk)
    
    # bisulfate association constant at T, S, P
    Khso4TPS = Khso4 * exp(lnKhso4fac)
    
    # GAMMA +/- HCl, activity coefficient of HCl at T/S, P=1
    # ADH is the Debye Huckel constant, calcualted as a polynomial 
    # fit to data in Khoo et al. 1977, doi:10.1021/ac50009a016
    # See Martz et al. 2010, DOI 10.4319/lom.2010.8.172, p175
    # Typo in paper 2nd term should be e-4 not e-6
    
    ADH = (3.4286e-6 * Temp^ 2 + 6.7524e-4 * Temp + 0.49172143)
    
    log10gammaHCl = -ADH * sqrt(IonS) / (1 + 1.394 * sqrt(IonS)) +  (0.08885 - 0.000111 * Temp) * IonS;
    # Millero 1983 partial molar volume of HCl in seawater
    deltaVHCl = 17.85 + 0.1044 * Temp - 0.001316 * Temp^ 2
    
    # effect of pressure on activity coefficient of HCl, divide by 2 because
    # its a mean activity coefficient, divide by 10 for units in the cm3 to F
    # conversion.
    
    log10gammaHCLtP = log10gammaHCl + deltaVHCl*(Press/10)/(R*Tk*ln10)/2/10
    
    #  Sensor reference potential
    
    # ************************************************************************
    k0T = k0 + k2 * Temp # % Temp  in deg C
    
    # CALCULATE PRESSURE CORRECTION (POLYNOMIAL FUNCTION OF PRESSURE)
    # ALL SENSORS HAVE A PRESSURE RESPONSE WHICH IS DETERMINED IN THE LAB
    # AND CONTAINED IN THE POLYNOMIAL Pcoefs
    f<-function(p){coefsp[1]*p*1+coefsp[2]*(p*1)^2+coefsp[3]*(p*1)^3+coefsp[4]*(p*1)^4+coefsp[5]*(p*1)^5+coefsp[6]*(p*1)^6}
    pcorr = f(Press)
    k0TP  = k0T + pcorr
    
    
    #  pH on free scale then corrected to get to pH total on mol/kg sw scale
    #    pHinsituFree = (Vrs - k0TP) / (R * Tk / F * ln10) + ...
    #                   log(Cltotal) / ln10 + 2 * log10gammaHCLtP
    #  this will be mol kg H2O  need to convert to mol/kg sw
    phfree = (Vrs - k0TP) / (R * Tk / Fa * ln10) + log(Cltotal) / ln10 + 2 * log10gammaHCLtP # %mol/kg-H2O scale
    
    # CONVERT TO mol/kg-sw scale - JP 2/4/16
    phfree = phfree - log10(1 - 0.001005 * Salt)  #mol/kg-sw scale
    
    # convert to total proton scale
    phtot = phfree - log10(1 + Stotal / Khso4TPS)
    
    
    
  }
  else {
    warning("Not enough data for pH processing in phase :",NumberPhase)
    phtot<-rep(NA,length(datapH$pH_mV))
  }
  
  return(phtot)
  
}

#**************************************************
#*
# Compute Ramses single spectra ######
#* Input : X = c(ramses_int_time,ramses_dark_count, I) 
#* where I is the spectra in count to be processed
#* 
#* Output : Physical units uW/cm2/nm vector with the same length than I
#*
#**************************************************

ra_single<-function(x,B0,B1,S,B0_Dark,B1_Dark){
t<-as.numeric(x[1]) #Integration time
offset<-as.numeric(x[2]) #Offset = Dark_average
I<-x[-(1:2)] #count spectra

#Etape1 Normalisation
M<-I/65535
  
#Etape2 Background Substraction
B<-B0 + t*B1/8192
C<-M-B

offset<-offset/65535
offset<-offset-B0_Dark-t*B1_Dark/8192

D<-C-offset

#Etape3 Integration time normalisation
E<-D*8192/t

#Change for uW/cm2/nm
E<-E/10
  
#Etape4 Sensitivity
return(E/S)
  
}

#**************************************************

Process_Ramses<-function(data,PixelStart=1,PixelStop=200,
                         PixelBinning=2,
                         calib_file="SAM.*AllCal.txt",
                         ramses_cal=NULL,
                         InWater="auto"){
  
  if (is.null(ramses_cal)){
    if (!file.exists(calib_file)){  
      
      #test 1 : avec pattern
      calib_file_pattern<-calib_file
      calib_file<-list.files(pattern = calib_file)[1]
      
      if (!file.exists(calib_file)){
        cat("!! No Ramses calibration file for: ",calib_file_pattern,"\n")
        
        #test 2 : generic
        calib_file<-list.files(pattern = "SAM.*AllCal.txt")[1]
        
        cat("!! Default Ramses calibration is used \n")
      }
      
    }
    
    if (file.exists(calib_file)){
      cat("Open RAMSES cal file: ",calib_file,"\n")
      ramses_cal<-read.table(calib_file,header = T,sep="\t")
    }
  }
  
if (!is.null(ramses_cal)){
  sq<-seq(PixelStart,PixelStop,by=PixelBinning)
  
  wave<-sapply(1:(length(sq)),function(i){mean(ramses_cal$Wave[c(sq[i],sq[i]+PixelBinning-1)])})
  B0<-sapply(1:(length(sq)),function(i){mean(ramses_cal$B0[c(sq[i],sq[i]+PixelBinning-1)])})
  B1<-sapply(1:(length(sq)),function(i){mean(ramses_cal$B1[c(sq[i],sq[i]+PixelBinning-1)])})
  
  indDark<-ramses_cal$Wave == -1
  B0_Dark=mean(ramses_cal$B0[indDark],na.rm = T)
  B1_Dark=mean(ramses_cal$B1[indDark],na.rm = T)
  
  ind<-c(grep("ramses_int_time",colnames(data)),grep("ramses_dark_count",colnames(data)),grep("ramses_raw_count",colnames(data)))
  
  dataCal<-NULL
  
  if (InWater=="yes"){
    S<-sapply(1:(length(sq)),function(i){mean(ramses_cal$S[c(sq[i],sq[i]+PixelBinning-1)])})
    cat("Ramses : apply in water calibration \n")
    #Process calibration
    dataCal<-t(apply(data[,ind],1,ra_single,B0=B0,B1=B1,S=S,B0_Dark,B1_Dark))
    
    sum_inAir=0
  } 
  
  if (InWater=="no"){
    S<-sapply(1:(length(sq)),function(i){mean(ramses_cal$Sair[c(sq[i],sq[i]+PixelBinning-1)])})
    cat("Ramses : apply in Air calibration \n")
    
    #Process calibration
    dataCal<-t(apply(data[,ind],1,ra_single,B0=B0,B1=B1,S=S,B0_Dark,B1_Dark))
    
    sum_inAir=nrow(dataCal)
  }
  
  if (InWater=="auto"){
    S<-sapply(1:(length(sq)),function(i){mean(ramses_cal$S[c(sq[i],sq[i]+PixelBinning-1)])})
    Sair<-sapply(1:(length(sq)),function(i){mean(ramses_cal$Sair[c(sq[i],sq[i]+PixelBinning-1)])})
    cat("Ramses : apply auto in air calibration \n")
    
    #Process calibration in water
    inAir<-data$PhaseName=="SUR"
    dataCal<-t(apply(data[!inAir,ind],1,ra_single,B0=B0,B1=B1,S=S,B0_Dark,B1_Dark))
    
    #Process calibration in air
    sum_inAir<-sum(inAir)
    cat("Data in air:",sum_inAir,"\n")
    dataCal_air<-t(apply(data[inAir,ind],1,ra_single,B0=B0,B1=B1,S=Sair,B0_Dark,B1_Dark))
    
    dataCal<-rbind(dataCal,dataCal_air)
  }
  

  
  wave<-round(wave*100)/100
  colnames(dataCal)<-paste("ramses_sig",wave,sep="_")
  
  #listOut
  listOut<-list(dataCal=dataCal,calib_file=calib_file,InWater=InWater,sum_inAir=sum_inAir)
}
else {
  warning("Ramses, no calibration file found:",calib_file,"\n")
  dataCal<-NULL
  #listOut
  listOut<-list(dataCal=dataCal,calib_file=NULL,InWater=InWater,sum_inAir=NULL)
}
  
return(listOut) 
  
}

#**************************************************
#*
# Compute IMU ######
#*
#**************************************************

# Calcul le Heading a partir des donnees raw Mag
# via Process_RawIMU
# Inputs RawMag et calibration
# Outputs Heading en degres

IMU_processHeading<-function(RawMag,imu_cal){
  
  compass_cal<-imu_cal$COMPASS
  mag_cal<-imu_cal$MAGNETOMETER
  
  # Compensation simple et Orientation
  PhyMagx=RawMag[1]+mag_cal$mx0
  PhyMagy=RawMag[3]+mag_cal$mz0
  PhyMagz=RawMag[2]+mag_cal$my0
  
  # Compensation compas
  PhyMagx=PhyMagx+compass_cal$hi1
  PhyMagy=PhyMagy+compass_cal$hi2
  
  PhyMagx=PhyMagx*compass_cal$si11 + PhyMagy*compass_cal$si12
  PhyMagy=PhyMagx*compass_cal$si21 + PhyMagy*compass_cal$si22
  
  # Calcul de l'angle
  fHead = atan2(PhyMagy,PhyMagx)
  
  # On retourne le résultat
  return( fHead *180.0 / pi )
  
}

#**************************************************
# Calcul le tilt et l'acceleration totale en fonction des 3 accelerations
# via Process_RawIMU
# Inputs RawAcc et calibration
# Outputs tilt en degres et acceleration en g

IMU_processAcc<-function(RawAcc,acc_cal){
  
  #Calibration et orientation
  PhyAccx=4*acc_cal$axg*(RawAcc[1]+acc_cal$ax0)/65536
  PhyAccy=4*acc_cal$azg*(RawAcc[3]+acc_cal$az0)/65536
  PhyAccz=4*acc_cal$ayg*(RawAcc[2]+acc_cal$ay0)/65536
  
  # Calcul du Tilt
  fTilt = atan2(sqrt(PhyAccx*PhyAccx + PhyAccy*PhyAccy),PhyAccz)
  fTilt = fTilt*180.0 / pi
  
  # Calcul du module de l'acceleration
  fAccTot = sqrt(PhyAccx*PhyAccx + PhyAccy*PhyAccy + PhyAccz*PhyAccz)
  
  
  # On retourne le résultat
  return( c(fTilt,fAccTot,PhyAccx,PhyAccy,PhyAccz) )
  
}

#**************************************************
# Calcul le heading, tilt et l'acceleration totale en fonction des 3 accelerations
# appelee par cts5_ProcessData
# Inputs : data=dataprofile$data$imu
#         imu_cal=Meta$SENSORS$SENSOR_IMU 
#
# Outputs data avec heading,tilt,acceleration

# data<-dataprofile$data$imu
# imu_cal<-Meta$SENSORS$SENSOR_IMU
Process_RawIMU<-function(data,imu_cal){
  heading<-NULL
  tilt<-NULL
  acceleration<-NULL
  PhyAccx<-NULL
  PhyAccy<-NULL
  PhyAccz<-NULL
  
  if (("ACCELEROMETER" %in% names(imu_cal)) & ("COMPASS" %in% names(imu_cal)) & ("MAGNETOMETER" %in% names(imu_cal))){
  
    indAcc<-grep("RawA",colnames(data))
    indMag<-grep("RawM",colnames(data))
    
    for (i in 1:nrow(data)){
      fheading<-IMU_processHeading(as.numeric(data[i,indMag]),imu_cal)
      
      tempacc<-IMU_processAcc(as.numeric(data[i,indAcc]),imu_cal$ACCELEROMETER)
      
      heading<-c(heading,fheading)
      tilt<-c(tilt,tempacc[1])
      acceleration<-c(acceleration,tempacc[2])
      PhyAccx<-c(PhyAccx,tempacc[3])
      PhyAccy<-c(PhyAccy,tempacc[4])
      PhyAccz<-c(PhyAccz,tempacc[5])
    }
    
    data<-cbind(data,heading,tilt,acceleration,PhyAccx,PhyAccy,PhyAccz)
  }
  else {
    warning("NO IMU calibration \n")
    data<-NULL
  }
    
  return(data)
}


#**************************************************
#*
# Compute PAL ######
#*
#**************************************************
# data<-dataprofile$data$pal

# MaNystuen2005_CS windcoef=c(-45.534,3.7163,-0.0984,0.0009)
# Nystuen2015_appendix windcoef=c(2.0871,0.4904,-0.031,0.0005)
# Med windcoef=c(107.93,-3.683,0.03143,0) ## Rapport Pochi
wind_PolySPL8<-function(data,windcoef=c(-45.534,3.7163,-0.0984,0.0009)){
  wind<-NULL
  if (!is.null(data$f_8000Hz)){
    SPL8<-data$f_8000Hz
    wind<-windcoef[1]+windcoef[2]*SPL8+windcoef[3]*(SPL8)^2+windcoef[4]*(SPL8)^3
  }
  
  return(wind)
}

rain_PolySPL5<-function(data,raincoef=c(-255.6,15.839,-0.328,0.0023)){
  rain<-NULL
  if (!is.null(data$f_5000Hz)){
    SPL5<-data$f_5000Hz
    rain<-raincoef[1]+raincoef[2]*SPL5+raincoef[3]*(SPL5)^2+raincoef[4]*(SPL5)^3
  }
  
  return(rain)
}
EdLeymarie/USEA_R documentation built on July 16, 2025, 1 p.m.