ODIN_CONA.Rmd

---
title: "ODIN pre-CONA deployment"
author: "Gustavo Olivares"
output: html_document
---
```{r load_libraries,echo=FALSE}
require('openair')
require('reshape2')
```
## Introduction
The purpose of this analysis is to document the performance of the ODIN during their pre-CONA deployment in Christchurch.  

## ODIN data
The units **ODIN_02-07** were deployed at Rangiora (**ODIN_02** failed to record data)
```{r load_odin}

## ODIN_03
odin_03 <- read.table("/home/gustavo/data/CONA/ODIN/deployment/odin_03.data",
                      header=T, quote="")
#force GMT as the time zone to avoid openair issues with daylight saving switches
#The actual time zone is 'NZST'
odin_03$date=as.POSIXct(paste(odin_03$Date,odin_03$Time),tz='GMT')
odin_03$Time<-NULL
odin_03$Date<-NULL
odin_03$Battery<-5*odin_03$Battery/1024

## ODIN_04
odin_04 <- read.table("/home/gustavo/data/CONA/ODIN/deployment/odin_04.data",
                      header=T, quote="")
#force GMT as the time zone to avoid openair issues with daylight saving switches
#The actual time zone is 'NZST'
odin_04$date=as.POSIXct(paste(odin_04$Date,odin_04$Time),tz='GMT')
odin_04$Time<-NULL
odin_04$Date<-NULL
odin_04$Battery<-5*odin_04$Battery/1024

## ODIN_05
odin_05 <- read.table("/home/gustavo/data/CONA/ODIN/deployment/odin_05.data",
                      header=T, quote="")
#force GMT as the time zone to avoid openair issues with daylight saving switches
#The actual time zone is 'NZST'
odin_05$date=as.POSIXct(paste(odin_05$Date,odin_05$Time),tz='GMT')
odin_05$Time<-NULL
odin_05$Date<-NULL
odin_05$Battery<-5*odin_05$Battery/1024

## ODIN_06
odin_06 <- read.table("/home/gustavo/data/CONA/ODIN/deployment/odin_06.data",
                      header=T, quote="")
#force GMT as the time zone to avoid openair issues with daylight saving switches
#The actual time zone is 'NZST'
odin_06$date=as.POSIXct(paste(odin_06$Date,odin_06$Time),tz='GMT')
odin_06$Time<-NULL
odin_06$Date<-NULL
odin_06$Battery<-5*odin_06$Battery/1024

## ODIN_07
odin_07 <- read.table("/home/gustavo/data/CONA/ODIN/deployment/odin_07.data",
                      header=T, quote="")
#force GMT as the time zone to avoid openair issues with daylight saving switches
#The actual time zone is 'NZST'
odin_07$date=as.POSIXct(paste(odin_07$Date,odin_07$Time),tz='GMT')
odin_07$Time<-NULL
odin_07$Date<-NULL
odin_07$Battery<-5*odin_07$Battery/1024

```

## ECan data
The data from the Rangiora site was obtained from Environment Canterbury's data catalogue and then corrected for proper date handling.

```{r load_ecan}

download.file(url = "http://data.ecan.govt.nz/data/29/Air/Air%20quality%20data%20for%20a%20monitored%20site%20(Hourly)/CSV?SiteId=5&StartDate=11%2F08%2F2015&EndDate=16%2F08%2F2015",destfile = "ecan_data.csv",method = "curl")
system("sed -i 's/a.m./AM/g' ecan_data.csv")
system("sed -i 's/p.m./PM/g' ecan_data.csv")
ecan_data_raw <- read.csv("ecan_data.csv",stringsAsFactors=FALSE)
ecan_data_raw$date<-as.POSIXct(ecan_data_raw$DateTime,format = "%d/%m/%Y %I:%M:%S %p",tz='GMT')
ecan_data<-as.data.frame(ecan_data_raw[,c('date','PM10.FDMS','Temperature..2m')])
names(ecan_data)<-c('date','PM10.FDMS','Temperature..1m')
```

## Merging the data
ECan's data was provided as 10 minute values while ODIN reports every 1 minute
so before merging the data, the timebase must be homogenized
```{r data_merge}
names(odin_03)<-c('Dust.03','Humidity.03','Temperature.03','Battery.03','date')
names(odin_04)<-c('Dust.04','Humidity.04','Temperature.04','Battery.04','date')
names(odin_05)<-c('Dust.05','Humidity.05','Temperature.05','Battery.05','date')
names(odin_06)<-c('Dust.06','Humidity.06','Temperature.06','Battery.06','date')
names(odin_07)<-c('Dust.07','Humidity.07','Temperature.07','Battery.07','date')

odin <- merge(odin_03,odin_04,by='date',all=TRUE)
odin <- merge(odin,odin_05,by='date',all=TRUE)
odin <- merge(odin,odin_06,by='date',all=TRUE)
odin <- merge(odin,odin_07,by='date',all=TRUE)
odin.10min<-timeAverage(selectByDate(odin,start = '2015-08-10',end = '2015-08-16'),avg.time='10 min')
all_merged.10min<-merge(odin.10min,ecan_data,by='date',all=TRUE)
```
## Time sync
```{r time_sync}
lag_test=ccf(all_merged.10min$Temperature.07,
             all_merged.10min$Temperature..1m,
             na.action=na.pass,
             lag.max=100,
             type='covariance',
             ylab='Correlation',
             main='Temperature correlation as function of clock lag')
odin_lag=lag_test$lag[which.max(lag_test$acf)]
```
ECan's record is behind by `r sprintf('%d',10*odin_lag)` minutes with respect
to ODIN data.  
The correction was applied to the ODIN data as follows:  
```{r correct_timing}
odin$date=odin$date-odin_lag*10*60
odin.10min<-timeAverage(selectByDate(odin,start = '2015-08-10',end = '2015-08-16'),avg.time='10 min')
all_merged.10min<-merge(odin.10min,ecan_data,by='date',all=TRUE)
lag_test=ccf(all_merged.10min$Temperature.03,
             all_merged.10min$Temperature..1m,
             na.action=na.pass,
             lag.max=100,
             type='covariance',
             ylab='Correlation',
             main='Temperature correlation as function of clock lag')
```

## Remove drift from ODIN raw data
It has been documented that the dust sensors suffer from significant drift, therefore a linear fit of the baseline drift is estimated.
```{r}
# Estimate the baseline from a simple linear regression
all_merged.10min$ODIN_drift.03<-predict(lm(all_merged.10min$Dust.03~seq(all_merged.10min$Dust.03)),newdata = all_merged.10min)
all_merged.10min$ODIN_drift.04<-predict(lm(all_merged.10min$Dust.04~seq(all_merged.10min$Dust.04)),newdata = all_merged.10min)
all_merged.10min$ODIN_drift.05<-predict(lm(all_merged.10min$Dust.05~seq(all_merged.10min$Dust.05)),newdata = all_merged.10min)
all_merged.10min$ODIN_drift.06<-predict(lm(all_merged.10min$Dust.06~seq(all_merged.10min$Dust.06)),newdata = all_merged.10min)
all_merged.10min$ODIN_drift.07<-predict(lm(all_merged.10min$Dust.07~seq(all_merged.10min$Dust.07)),newdata = all_merged.10min)

# Remove the baseline drift from the raw ODIN data

all_merged.10min$Dust.03.raw <- all_merged.10min$Dust.03
all_merged.10min$Dust.03.detrend<-all_merged.10min$Dust.03.raw - all_merged.10min$ODIN_drift.03

all_merged.10min$Dust.04.raw <- all_merged.10min$Dust.04
all_merged.10min$Dust.04.detrend<-all_merged.10min$Dust.04.raw - all_merged.10min$ODIN_drift.04

all_merged.10min$Dust.05.raw <- all_merged.10min$Dust.05
all_merged.10min$Dust.05.detrend<-all_merged.10min$Dust.05.raw - all_merged.10min$ODIN_drift.05

all_merged.10min$Dust.06.raw <- all_merged.10min$Dust.06
all_merged.10min$Dust.06.detrend<-all_merged.10min$Dust.06.raw - all_merged.10min$ODIN_drift.06

all_merged.10min$Dust.07.raw <- all_merged.10min$Dust.07
all_merged.10min$Dust.07.detrend<-all_merged.10min$Dust.07.raw - all_merged.10min$ODIN_drift.07

## Testing not correcting drift
all_merged.10min$Dust.03.detrend<-all_merged.10min$Dust.03.raw
all_merged.10min$Dust.04.detrend<-all_merged.10min$Dust.04.raw
all_merged.10min$Dust.05.detrend<-all_merged.10min$Dust.05.raw
all_merged.10min$Dust.06.detrend<-all_merged.10min$Dust.06.raw
all_merged.10min$Dust.07.detrend<-all_merged.10min$Dust.07.raw

```

## Calculate the temperature interference
First we divide the temperature range for each unit.
```{r}
all_merged.10min$Temperature.03.bin<-cut(all_merged.10min$Temperature.03,breaks = c(0,5,10,15,20,25),labels = c('2.5','7.5','12.5','17.5','22.5'))
all_merged.10min$Temperature.04.bin<-cut(all_merged.10min$Temperature.04,breaks = c(0,5,10,15,20,25),labels = c('2.5','7.5','12.5','17.5','22.5'))
all_merged.10min$Temperature.05.bin<-cut(all_merged.10min$Temperature.05,breaks = c(0,5,10,15,20,25),labels = c('2.5','7.5','12.5','17.5','22.5'))
all_merged.10min$Temperature.06.bin<-cut(all_merged.10min$Temperature.06,breaks = c(0,5,10,15,20,25),labels = c('2.5','7.5','12.5','17.5','22.5'))
all_merged.10min$Temperature.07.bin<-cut(all_merged.10min$Temperature.07,breaks = c(0,5,10,15,20,25),labels = c('2.5','7.5','12.5','17.5','22.5'))
Temp <- c(2.5,7.5,12.5,17.5,22.5)

Dust.03<-tapply(all_merged.10min$Dust.03.detrend,all_merged.10min$Temperature.03.bin,quantile,0.25)
Dust.04<-tapply(all_merged.10min$Dust.04.detrend,all_merged.10min$Temperature.04.bin,quantile,0.25)
Dust.05<-tapply(all_merged.10min$Dust.05.detrend,all_merged.10min$Temperature.05.bin,quantile,0.25)
Dust.06<-tapply(all_merged.10min$Dust.06.detrend,all_merged.10min$Temperature.06.bin,quantile,0.25)
Dust.07<-tapply(all_merged.10min$Dust.07.detrend,all_merged.10min$Temperature.07.bin,quantile,0.25)

TC_Dust.03 <- data.frame(Dust.03.detrend = Dust.03,Temperature.03 = Temp)
TC_Dust.04 <- data.frame(Dust.04.detrend = Dust.04,Temperature.04 = Temp)
TC_Dust.05 <- data.frame(Dust.05.detrend = Dust.05,Temperature.05 = Temp)
TC_Dust.06 <- data.frame(Dust.06.detrend = Dust.06,Temperature.06 = Temp)
TC_Dust.07 <- data.frame(Dust.07.detrend = Dust.07,Temperature.07 = Temp)

```
Now we calculate the linear regression for the minimum dust response in each temperature bin and subtract it from the detrended data
```{r}
summary(odin.03_T<-lm(data = TC_Dust.03,Dust.03.detrend~Temperature.03))
summary(odin.04_T<-lm(data = TC_Dust.04,Dust.04.detrend~Temperature.04))
summary(odin.05_T<-lm(data = TC_Dust.05,Dust.05.detrend~Temperature.05))
summary(odin.06_T<-lm(data = TC_Dust.06,Dust.06.detrend~Temperature.06))
summary(odin.07_T<-lm(data = TC_Dust.07,Dust.07.detrend~Temperature.07))

all_merged.10min$Dust.03.corr <- all_merged.10min$Dust.03.detrend - predict(odin.03_T,newdata = all_merged.10min)
all_merged.10min$Dust.04.corr <- all_merged.10min$Dust.04.detrend - predict(odin.04_T,newdata = all_merged.10min)
all_merged.10min$Dust.05.corr <- all_merged.10min$Dust.05.detrend - predict(odin.05_T,newdata = all_merged.10min)
all_merged.10min$Dust.06.corr <- all_merged.10min$Dust.06.detrend - predict(odin.06_T,newdata = all_merged.10min)
all_merged.10min$Dust.07.corr <- all_merged.10min$Dust.07.detrend - predict(odin.07_T,newdata = all_merged.10min)

```



## Dust performance using ECan data for calibration
With ECan's PM data available, a more accurate calibration can be applied
to the corrected Dust signal from ODIN.

$Dust_{calibrated}=A*Dust_{corrected}+B$


### Full dataset 1 hour  PM$_{2.5}$ fdms
```{r full_1_pm25}
all_merged.1hr<-timeAverage(all_merged.10min,avg.time='1 hour')

summary(odin3.lm.full.1hr.pm2.5<-
          lm(data=all_merged.1hr,PM10.FDMS~
               Dust.03.corr))
summary(odin4.lm.full.1hr.pm2.5<-
          lm(data=all_merged.1hr,PM10.FDMS~
               Dust.04.corr))
summary(odin5.lm.full.1hr.pm2.5<-
          lm(data=all_merged.1hr,PM10.FDMS~
               Dust.05.corr))
summary(odin6.lm.full.1hr.pm2.5<-
          lm(data=all_merged.1hr,PM10.FDMS~
               Dust.06.corr))
summary(odin7.lm.full.1hr.pm2.5<-
          lm(data=all_merged.1hr,PM10.FDMS~
               Dust.07.corr))
```
### Calibrated Dust
```{r}

all_merged.1hr$Dust.03.cal<-predict(odin3.lm.full.1hr.pm2.5,newdata = all_merged.1hr)
all_merged.1hr$Dust.04.cal<-predict(odin4.lm.full.1hr.pm2.5,newdata = all_merged.1hr)
all_merged.1hr$Dust.05.cal<-predict(odin5.lm.full.1hr.pm2.5,newdata = all_merged.1hr)
all_merged.1hr$Dust.06.cal<-predict(odin6.lm.full.1hr.pm2.5,newdata = all_merged.1hr)
all_merged.1hr$Dust.07.cal<-predict(odin7.lm.full.1hr.pm2.5,newdata = all_merged.1hr)

timePlot(all_merged.1hr,pollutant = c('PM10.FDMS',
                                      'Dust.03.corr',
                                      'Dust.04.corr',
                                      'Dust.05.corr',
                                      'Dust.06.corr',
                                      'Dust.07.corr')
         ,group = TRUE)
timePlot(all_merged.1hr,pollutant = c('PM10.FDMS',
                                      'Dust.03.cal',
                                      'Dust.04.cal',
                                      'Dust.05.cal',
                                      'Dust.06.cal',
                                      'Dust.07.cal')
         ,group = TRUE)
timeVariation(all_merged.1hr,pollutant = c('PM10.FDMS',
                                           'Dust.03.cal',
                                           'Dust.04.cal',
                                           'Dust.05.cal',
                                           'Dust.06.cal',
                                           'Dust.07.cal'))
timeVariation(all_merged.1hr,pollutant = c('PM10.FDMS',
                                      'Dust.03.corr',
                                      'Dust.04.corr',
                                      'Dust.05.corr',
                                      'Dust.06.corr',
                                      'Dust.07.corr'))
timeVariation(all_merged.1hr,pollutant = c('Dust.03.corr',
                                      'Dust.04.corr',
                                      'Dust.05.corr',
                                      'Dust.06.corr',
                                      'Dust.07.corr'))
```