IEC 2018 plots

Script to generate plots for the IEC2018 presentation about ATCHU

Libraries

library(librarian)
shelf(openair,
      reshape2,
      ggplot2,
      readr,
      latex2exp)

Warning in shelf(openair, reshape2, ggplot2, readr, latex2exp): cran_repo = '@CRAN@' is not a valid URL. 
                    Defaulting to cran_repo = 'https://cran.r-project.org'.

Load data

data.wide <- read_delim("~/data/ATCHU/ATCHU_dataset.csv",delim = ',',
                        col_types = cols(col_datetime(format = ""),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double(),
                                      col_double()))

Warning in rbind(names(probs), probs_f): number of columns of result is not
a multiple of vector length (arg 1)

Warning: 16635 parsing failures.
row # A tibble: 5 x 5 col     row col   expected   actual     file                             expected   <int> <chr> <chr>      <chr>      <chr>                            actual 1     1 <NA>  16 columns 17 columns '~/data/ATCHU/ATCHU_dataset.csv' file 2     2 <NA>  16 columns 17 columns '~/data/ATCHU/ATCHU_dataset.csv' row 3     3 <NA>  16 columns 17 columns '~/data/ATCHU/ATCHU_dataset.csv' col 4     4 <NA>  16 columns 17 columns '~/data/ATCHU/ATCHU_dataset.csv' expected 5     5 <NA>  16 columns 17 columns '~/data/ATCHU/ATCHU_dataset.csv'
... ................. ... .................................................................... ........ .................................................................... ...... .................................................................... .... .................................................................... ... .................................................................... ... .................................................................... ........ ....................................................................
See problems(...) for more details.

attr(data.wide$date,"tzone") <- "NZST"

Calculate extra vars

data.wide$NOxHS <- data.wide$NOHS + data.wide$NO2HS
data.wide$NOxMP <- data.wide$NOMP + data.wide$NO2MP
data.wide$BC.NOxHS <- data.wide$BCHS / data.wide$NOxHS

Load Grimm data for HS

load("~/data/ATCHU/grimm_HS.RData")
Dp.HS <- Dp
names(grimm_dndlogdp.HS) <- c(names(grimm_dndlogdp.HS)[1:2],paste0("HS",names(grimm_dndlogdp.HS)[3:18]))
data.wide <- merge(data.wide,grimm_dndlogdp.HS,by = 'date',all = TRUE)

Load Grimm data for MP

load("~/data/ATCHU/grimm_MP.RData")
names(grimm_dndlogdp.MP) <- c(names(grimm_dndlogdp.MP)[1:2],paste0("MP",names(grimm_dndlogdp.MP)[3:33]))
data.wide <- merge(data.wide,grimm_dndlogdp.MP,by = 'date',all = TRUE)

Load SMPS data from UoA

load("~/data/ATCHU/smps_UA.RData")
Dp.smps.UA <- Dp.UA
data.wide <- merge(data.wide,smps_10min_wide,by = 'date',all = TRUE)

Load 2B NOx data from UoA

parsed2B <- read_csv("~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.txt",
                     col_types = cols(date = col_datetime(format = "%Y/%m/%d %H:%M:%S GMT")))

Warning in rbind(names(probs), probs_f): number of columns of result is not
a multiple of vector length (arg 1)

Warning: 1435 parsing failures.
row # A tibble: 5 x 5 col     row col   expected   actual     file                                   expected   <int> <chr> <chr>      <chr>      <chr>                                  actual 1 11460 O3    an integer 4294966000 '~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.… file 2 11461 O3    an integer 4294966000 '~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.… row 3 11462 O3    an integer 4294966000 '~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.… col 4 11463 O3    an integer 4294966000 '~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.… expected 5 11464 O3    an integer 4294966000 '~/data/ATCHU/RAW/UoA/2BNOX/parsed2B.…
... ................. ... .......................................................................... ........ .......................................................................... ...... .......................................................................... .... .......................................................................... ... .......................................................................... ... .......................................................................... ........ ..........................................................................
See problems(...) for more details.

nox.UA <- subset(parsed2B,(date < as.POSIXct("2015-04-29 20:00:00")))[,1:2]
names(nox.UA) <- c('date','NOx.UA')
data.wide <- merge(data.wide,nox.UA,by='date',all = TRUE)

Load Grimm data for UoA

load("~/data/ATCHU/grimm_UA.RData")
data.wide <- merge(data.wide,grimm_dndlogdp.UA, by = 'date', all = TRUE)

Average up to 1 hour

data.wide.1hr <- timeAverage(data.wide,avg.time = '1 hour')

Correcting the scales for the Grimm data

data.wide.1hr[,21:35] <- data.wide.1hr[,21:35] / 10 # Halsley Drive
data.wide.1hr[,38:68] <- data.wide.1hr[,38:68] / 100 # Musick Point
data.wide.1hr[,172:202] <- data.wide.1hr[,172:202] / 100 #Auckland University

Plots

All pollution roses are normalised so that the fraction of measurements of certain level is shown.

Describe Halsley Drive

HSn300 = dN/dlogDp at 300nm from the Grimm

timePlot(data.wide.1hr,pollutant = c('BCHS','NOHS','NO2HS','O3HS','wsHS','HSn300') , y.relation = 'free')

timeVariation(data.wide.1hr,
              pollutant = c('BCHS','NOxHS'),
              normalise = TRUE)

pollutionRose(data.wide.1hr,
              pollutant = 'BCHS',
              ws='wsHS',
              wd='wdHS',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,500,1500,3000),
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'NOHS',
              ws='wsHS',
              wd='wdHS',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,2,4,6),
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'HSn300',
              ws='wsHS',
              wd='wdHS',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,2,4,6),
              statistic = 'prop.count')

Describe UoA

50 = dN/dLogDp at 50nm from the SMPS n265 = dN/dLogDp at 265nm from the Grimm

timePlot(data.wide.1hr,pollutant = c('NOx.UA','BCUoA','50','n265'),y.relation = 'free')

pollutionRose(data.wide.1hr,
              pollutant = 'BCUoA',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,500,1500,3000),
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'NOx.UA',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,2,5,10,20),
              statistic = 'prop.count')

Warning in windRose(mydata, pollutant = pollutant, paddle = paddle, seg =
seg, : Some values are below minimum break.

pollutionRose(data.wide.1hr,
              pollutant = '50',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'n265',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              statistic = 'prop.count')

Describe Musick Point

n265 = dN/dLogDp at 265nm from the Grimm

timePlot(data.wide.1hr,pollutant = c('NOMP','NO2MP','N6MP','O3MP','wsMP','MPn265'),y.relation = 'free')

pollutionRose(data.wide.1hr,
              pollutant = 'NOMP',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,2,4,6),
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'N6MP',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              breaks = c(0,2000,5000,10000,20000),
              statistic = 'prop.count')

pollutionRose(data.wide.1hr,
              pollutant = 'MPn265',
              ws='wsMP',
              wd='wdMP',
              normalise = TRUE,
              cols = 'heat',
              statistic = 'prop.count')

Compare Size Distributions by wind sector

sw_sector.HS <- (data.wide.1hr$wdHS < 255) & (data.wide.1hr$wdHS > 195)
ne_sector.HS <- (data.wide.1hr$wdHS < 75) & (data.wide.1hr$wdHS > 15)
sizedist_sw.HS <- colMeans(data.wide.1hr[sw_sector.HS,21:35],na.rm = TRUE)
sizedist_ne.HS <- colMeans(data.wide.1hr[ne_sector.HS,21:35],na.rm = TRUE)

sw_sector.MP <- (data.wide.1hr$wdMP < 255) & (data.wide.1hr$wdMP > 195)
ne_sector.MP <- (data.wide.1hr$wdMP < 75) & (data.wide.1hr$wdMP > 15)
sizedist_sw.MP <- colMeans(data.wide.1hr[sw_sector.MP,38:68],na.rm = TRUE)
sizedist_ne.MP <- colMeans(data.wide.1hr[ne_sector.MP,38:68],na.rm = TRUE)

sizedist_sw.UA <- colMeans(data.wide.1hr[sw_sector.MP,70:169], na.rm = TRUE)
sizedist_ne.UA <- colMeans(data.wide.1hr[ne_sector.MP,70:169], na.rm = TRUE)

sizedist_sw.grimm.UA <- colMeans(data.wide.1hr[sw_sector.MP,172:202], na.rm = TRUE)
sizedist_ne.grimm.UA <- colMeans(data.wide.1hr[ne_sector.MP,172:202], na.rm = TRUE)

NOxsizedist_sw.UA <- colMeans(data.wide.1hr[sw_sector.MP,70:169]/data.wide.1hr$NOx.UA[sw_sector.MP], na.rm = TRUE)
NOxsizedist_ne.UA <- colMeans(data.wide.1hr[ne_sector.MP,70:169]/data.wide.1hr$NOx.UA[ne_sector.MP], na.rm = TRUE)

cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

Comparing size distribution between Grimm and SMPS at the university (Sanity check)

sanity1 <- timeAverage(smps_10min_wide,avg.time = '1 hour')
sanity2 <- timeAverage(grimm_dndlogdp.UA,avg.time = '1 hour')
sanity.check <- merge(sanity1,sanity2,by = 'date')

sanity.sizedist_sw.UA <- colMeans(sanity.check[,3:102], na.rm = TRUE)
sanity.sizedist_sw.grimm.UA <- colMeans(sanity.check[,104:134], na.rm = TRUE) / 100
ggplot() +
  geom_line(aes(x=Dp.smps.UA,y=sanity.sizedist_sw.UA,colour = 'UoA SMPS')) +
  geom_line(aes(x=Dp.UA,y=sanity.sizedist_sw.grimm.UA,colour = 'UoA Grimm')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette) +
  ylab(TeX("$(\\frac{dN}{dlogD_p})\\[pt/cc\\]$")) + 
  xlab(TeX("$D_p\\[nm\\]"))

University roof

ggplot() +
  geom_line(aes(x=Dp.smps.UA,y=sizedist_ne.UA,colour = 'UoA North-East (smps)')) +
  geom_line(aes(x=Dp.smps.UA,y=sizedist_sw.UA,colour = 'UoA South-West (smps)')) +
  geom_line(aes(x=Dp.UA,y=sizedist_ne.grimm.UA,colour = 'UoA North-East (grimm)')) +
  geom_line(aes(x=Dp.UA,y=sizedist_sw.grimm.UA,colour = 'UoA South-West (grimm)')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_sw.HS,colour = 'HS South-West')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_ne.HS,colour = 'HS North-East')) +
#  geom_line(aes(x=Dp.MP,y=sizedist_sw.MP,colour = 'MP South-West (grimm)')) +
#  geom_line(aes(x=Dp.MP,y=sizedist_ne.MP,colour = 'MP North-East (grimm)')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette) +
  ylab(TeX("$(\\frac{dN}{dlogD_p})\\[pt/cc\\]$")) + 
  xlab(TeX("$D_p\\[nm\\]"))

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 1 rows containing missing values (geom_path).

Musick Point

ggplot() +
  #geom_line(aes(x=Dp.smps.UA,y=sizedist_ne.UA,colour = 'UoA North-East (smps)')) +
  #geom_line(aes(x=Dp.smps.UA,y=sizedist_sw.UA,colour = 'UoA South-West (smps)')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_sw.HS,colour = 'HS South-West')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_ne.HS,colour = 'HS North-East')) +
  geom_line(aes(x=Dp.MP,y=sizedist_sw.MP,colour = 'MP South-West (grimm)')) +
  geom_line(aes(x=Dp.MP,y=sizedist_ne.MP,colour = 'MP North-East (grimm)')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette) +
  ylab(TeX("$(\\frac{dN}{dlogD_p})\\[pt/cc\\]$")) + 
  xlab(TeX("$D_p\\[nm\\]"))

Halsley School

ggplot() +
  #geom_line(aes(x=Dp.smps.UA,y=sizedist_ne.UA,colour = 'UoA North-East (smps)')) +
  #geom_line(aes(x=Dp.smps.UA,y=sizedist_sw.UA,colour = 'UoA South-West (smps)')) +
  geom_line(aes(x=Dp.HS,y=sizedist_sw.HS,colour = 'HS South-West')) +
  geom_line(aes(x=Dp.HS,y=sizedist_ne.HS,colour = 'HS North-East')) +
  #geom_line(aes(x=Dp.MP,y=sizedist_sw.MP,colour = 'MP South-West')) +
  #geom_line(aes(x=Dp.MP,y=sizedist_ne.MP,colour = 'MP North-East')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette)

Warning: Removed 3 rows containing missing values (geom_path).

Warning: Removed 2 rows containing missing values (geom_path).

Halsley School and Musick Point

ggplot() +
#  geom_line(aes(x=Dp.smps.UA,y=sizedist_ne.UA,colour = 'UoA North-East')) +
#  geom_line(aes(x=Dp.smps.UA,y=sizedist_sw.UA,colour = 'UoA South-West')) +
  geom_line(aes(x=Dp.HS,y=sizedist_sw.HS,colour = 'HS South-West')) +
  geom_line(aes(x=Dp.HS,y=sizedist_ne.HS,colour = 'HS North-East')) +
  geom_line(aes(x=Dp.MP,y=sizedist_sw.MP,colour = 'MP South-West')) +
  geom_line(aes(x=Dp.MP,y=sizedist_ne.MP,colour = 'MP North-East')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette)

Warning: Removed 3 rows containing missing values (geom_path).

Warning: Removed 2 rows containing missing values (geom_path).

Musick Point and Auckland University

ggplot() +
  geom_line(aes(x=Dp.smps.UA,y=sizedist_ne.UA,colour = 'UoA North-East')) +
  geom_line(aes(x=Dp.smps.UA,y=sizedist_sw.UA,colour = 'UoA South-West')) +
  geom_line(aes(x=Dp.UA,y=sizedist_ne.grimm.UA,colour = 'UoA North-East')) +
  geom_line(aes(x=Dp.UA,y=sizedist_sw.grimm.UA,colour = 'UoA South-West')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_sw.HS,colour = 'HS South-West')) +
  #geom_line(aes(x=Dp.HS,y=sizedist_ne.HS,colour = 'HS North-East')) +
  geom_line(aes(x=Dp.MP,y=sizedist_sw.MP,colour = 'MP South-West')) +
  geom_line(aes(x=Dp.MP,y=sizedist_ne.MP,colour = 'MP North-East')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette) +
  ylab(TeX("$(\\frac{dN}{dlogD_p})\\[pt/cc\\]$")) + 
  xlab(TeX("$D_p\\[nm\\]"))

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 1 rows containing missing values (geom_path).

(dN/dlogDp)/NOx at University

ggplot() +
  geom_line(aes(x=Dp.smps.UA,y=NOxsizedist_ne.UA,colour = 'UoA North-East')) +
  geom_line(aes(x=Dp.smps.UA,y=NOxsizedist_sw.UA,colour = 'UoA South-West')) +
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(trans = 'log10') +
  scale_colour_manual(values = cbbPalette) +
  ylab(TeX("$\\frac{(\\frac{dN}{dlogD_p})}{NO_X}\\[pt/cc\\]$")) + 
  xlab(TeX("$D_p\\[nm\\]"))