volcano_2011_ECan.R

# Did the volcano hit Canterbury in 2011?
library('ggplot2')
library('openair')

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

library('lubridate')

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

library('corrgram')
library('corrplot')

# Get ECan Data

# SiteNo,SiteName,LongName,D_S_Sites_Key,AirShed
# 1,St Albans,Christchurch - St Albans,1,Christchurch
# 2,Riccarton Road,Christchurch - Riccarton Road,2,Christchurch
# 3,Woolston,Christchurch - Woolston,3,Christchurch
# 4,Kaiapoi,Kaiapoi - Kaiapoi,4,Kaiapoi
# 5,Rangiora,Rangiora - Rangiora,5,Rangiora
# 6,Lincoln,Lincoln - Lincoln,6,Lincoln
# 7,Ashburton,Ashburton - Ashburton,7,Ashburton
# 9,Geraldine,Geraldine - Geraldine,9,Geraldine
# 10,Timaru Anzac Square,Timaru - Timaru Anzac Square,10,Timaru
# 11,Washdyke,Timaru - Washdyke,11,Timaru
# 12,Waimate Stadium,Waimate - Waimate Stadium,12,Waimate
# 54,Timaru Grey Rd,Timaru - Timaru Grey Rd,54,Timaru

site_id <- c('1','2','3','4','5','6','7','9','10','11','12','54') # ECan sites
xx <- 0
for (site in site_id){
  download.file(url = paste0("http://data.ecan.govt.nz/data/29/Air/Air%20quality%20data%20for%20a%20monitored%20site%20(Hourly)/CSV?SiteId=",site,"&StartDate=01%2F05%2F2011&EndDate=30%2F09%2F2011"),destfile = paste0("ecan_data_",site,".csv"),method = "curl")
  system(paste0("sed -i 's/a.m./AM/g' ecan_data_",site,".csv"))
  system(paste0("sed -i 's/p.m./PM/g' ecan_data_",site,".csv"))
  if (file.info(paste0("ecan_data_",site,".csv"))$size>100000) {
    ecan_data_raw <- read.csv(paste0("./ecan_data_",site,".csv"),stringsAsFactors=FALSE)
    ecan_data_raw$date<-as.POSIXct(ecan_data_raw$DateTime,format = "%d/%m/%Y %H:%M:%S",tz='NZST')
    ecan_data_raw <- timeAverage(ecan_data_raw,avg.time = '10 min')
    idx_pm10 <- grep('PM10',names(ecan_data_raw))
    idx_pm2.5 <- grep('PM2',names(ecan_data_raw))
    idx_so2 <- grep('SO2',names(ecan_data_raw))
    idx_date <- grep('date',names(ecan_data_raw))
    idx <- c(idx_date,idx_pm10,idx_pm2.5,idx_so2)
    if (length(idx)<=1){
      next
    }
    eval(parse(text=paste0("ecan_data_raw_",site," <- ecan_data_raw[,names(ecan_data_raw)[idx]]")))
    eval(parse(text=paste0("names(ecan_data_raw_",site,") <- c('date',paste0(names(ecan_data_raw_",site,")[2:length(idx)],'_",site,"'))")))
    if (xx==0){
      eval(parse(text=paste0("merged_data <- ecan_data_raw_",site)))
      xx <- xx +1
    }
    else {
      eval(parse(text=paste0("merged_data <- merge(merged_data,ecan_data_raw_",site,",by='date', all=TRUE)")))
    }
  }
}

# Plots for the extended period (May to September)

timePlot(merged_data,pollutant = names(merged_data)[grep('SO',names(merged_data))],group = TRUE,avg.time = '1 day')

timePlot(merged_data,pollutant = names(merged_data)[grep('PM10',names(merged_data))],group = TRUE,avg.time = '1 day')

timePlot(merged_data,pollutant = names(merged_data)[grep('PM2',names(merged_data))],group = TRUE,avg.time = '1 day')

# Plots for the period around the 14th of June
short_data <- selectByDate(merged_data,start = '2011-06-10', end = '2011-06-22')

timePlot(short_data,pollutant = names(short_data)[grep('SO',names(short_data))],group = TRUE, avg.time = '1 hour')

timePlot(short_data,pollutant = names(short_data)[grep('PM10',names(short_data))],group = TRUE, avg.time = '1 hour')

timePlot(short_data,pollutant = names(short_data)[grep('PM2',names(short_data))],group = TRUE, avg.time = '1 hour')