Wind Analysis (with code)
Danielle Boccelli
12/5/2018
This analysis considers the data from a PurpleAir sensor in Point Richmond from 2017-09-30 to 2018-12-02.
#Preparing data for visualization
#Generating Date from Epoch information
c[,"Date"] = as.Date(as.POSIXct(c$Epoch, origin='1970-01-01', tz="Etc/GMT+5"))
#Reordering Data
c = c[order(c$Date),]
c = c[,c(25,2:9,11,13:17,20,23:24,10,12,18:19,21:22)]
#Remove rows that containing missing information
c = c[complete.cases(c),]#Generate standard directions based on wind direction angle
#Creating binary variables
c[,"S"] = ifelse(c[,"WindDir"] <= 11.25 | c[,"WindDir"] > 348.75, 1, 0)
c[,"SSW"] = ifelse(c[,"WindDir"] <= 33.75 & c[,"WindDir"] > 11.25, 1, 0)
c[,"SW"] = ifelse(c[,"WindDir"] <= 56.25 & c[,"WindDir"] > 33.75, 1, 0)
c[,"WSW"] = ifelse(c[,"WindDir"] <= 78.75 & c[,"WindDir"] > 56.25, 1, 0)
c[,"W"] = ifelse(c[,"WindDir"] <= 101.25 & c[,"WindDir"] > 78.75, 1, 0)
c[,"WNW"] = ifelse(c[,"WindDir"] <= 123.75 & c[,"WindDir"] > 101.25, 1, 0)
c[,"NW"] = ifelse(c[,"WindDir"] <= 146.25 & c[,"WindDir"] > 123.75, 1, 0)
c[,"NNW"] = ifelse(c[,"WindDir"] <= 168.75 & c[,"WindDir"] > 146.25, 1, 0)
c[,"N"] = ifelse(c[,"WindDir"] <= 191.25 & c[,"WindDir"] > 168.75, 1, 0)
c[,"NNE"] = ifelse(c[,"WindDir"] <= 213.75 & c[,"WindDir"] > 191.25, 1, 0)
c[,"NE"] = ifelse(c[,"WindDir"] <= 236.25 & c[,"WindDir"] > 213.75, 1, 0)
c[,"ENE"] = ifelse(c[,"WindDir"] <= 258.75 & c[,"WindDir"] > 236.25, 1, 0)
c[,"E"] = ifelse(c[,"WindDir"] <= 281.25 & c[,"WindDir"] > 258.75, 1, 0)
c[,"ESE"] = ifelse(c[,"WindDir"] <= 303.75 & c[,"WindDir"] > 281.25, 1, 0)
c[,"SE"] = ifelse(c[,"WindDir"] <= 326.25 & c[,"WindDir"] > 303.75, 1, 0)
c[,"SSE"] = ifelse(c[,"WindDir"] <= 348.75 & c[,"WindDir"] > 326.25, 1, 0)
#Creating a Wind Direction Category factor
c[,"WindDirCat"] = factor(colnames(c[26:40])[max.col(c[,26:40])],
levels = c("N","NNE","NE","ENE","E","ESE","SE","SSE",
"S","SSW","SW","WSW","W","WNW","NW","NNW"),
labels = c("N","NNE","NE","ENE","E","ESE","SE","SSE",
"S","SSW","SW","WSW","W","WNW","NW","NNW"))Wind Direction
A look at the direction of the wind in the data set. The graphs show that the wind usually blows from the bay toward the community.
#Create bar plot
g1 = ggplot(data = c,
aes(x = WindDirCat)) +
labs(y = "Number of measurements",
x = "") +
geom_bar(stat = "count") +
geom_rect(aes(xmin = 12, xmax = 15, ymin = 61000, ymax = 71500), alpha = 0.1, fill = "lightgray") +
annotate("text", label = "Refinery", x = 12.1, y = 70000, color = "red", hjust = 0) +
annotate("text", label = "Community", x = 12.1, y = 67500, color = "black", hjust = 0) +
annotate("text", label = "Park", x = 12.1, y = 65000, color = "#4A9216", hjust = 0) +
annotate("text", label = "Bay", x = 12.1, y = 62500, color = "#1577BB", hjust = 0) +
theme(axis.title = element_text(size = 8),
title = element_text(size=9.5),
axis.text.x = element_text(size = 7,
color = c("red","red",
"black","black","black","black",
"#4A9216","#4A9216","#4A9216",
"#1577BB","#1577BB","#1577BB","#1577BB","#1577BB",
"red","red")))
#Create radial plot
g2 = ggplot(data = c,
aes(x = WindDirCat)) +
geom_bar(stat = "count", col = "black", width = 0.5) +
coord_polar(theta = "x") +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
title = element_text(size=9.5),
axis.text.x = element_text(color = c("red","red",
"black","black","black","black",
"#4A9216","#4A9216","#4A9216",
"#1577BB","#1577BB","#1577BB","#1577BB","#1577BB",
"red","red")))
#Print plots
grid.arrange(g1, g2, ncol = 2)
All Chemicals
These plots show the average measurement for all chemical species by wind direction. These graphs show that the average measurements of the combination of all chemical species are highest when the wind is blowing from the refinery.
#Aggregate chemical data by wind direction
c_agg = aggregate(c[,2:18],
by = list(WindDirCat = c$WindDirCat),
FUN = mean)
c_agg = gather(c_agg, key = "Chemical",
value = "Mean",
-WindDirCat)
#Create bar plot
g1 = ggplot(data = c_agg,
aes(x = WindDirCat, y = Mean, fill = Chemical)) +
geom_bar(stat = "identity", col = "black", width = 0.5) +
geom_rect(aes(xmin = 12, xmax = 15, ymin = 24, ymax = 28), alpha = 0.1, fill = "lightgray") +
annotate("text", label = "Refinery", x = 12.1, y = 27.5, color = "red", hjust = 0) +
annotate("text", label = "Community", x = 12.1, y = 26.5, color = "black", hjust = 0) +
annotate("text", label = "Park", x = 12.1, y = 25.5, color = "#4A9216", hjust = 0) +
annotate("text", label = "Bay", x = 12.1, y = 24.5, color = "#1577BB", hjust = 0) +
labs(y = "Mean detection",
x = "") +
theme(axis.title = element_text(size = 8),
title = element_text(size=9.5),
legend.position = "none",
axis.text.x = element_text(size = 7,
color = c("red","red",
"black","black","black","black",
"#4A9216","#4A9216","#4A9216",
"#1577BB","#1577BB","#1577BB","#1577BB","#1577BB",
"red","red"))) +
scale_fill_manual(values=c("#F8766D", "#E7851E", "#D09400", "#B2A100", "#89AC00", "#45B500",
"#00BC51", "#00C087", "#00C0B2", "#00BCD6", "#00B3F2", "#29A3FF",
"#9C8DFF", "#D277FF", "#F166E8", "#FF61C7", "#FF689E"), drop = FALSE)
#Create radial plot
g2 = ggplot(data = c_agg,
aes(x = WindDirCat, y = Mean, fill = Chemical)) +
geom_bar(stat = "identity", col = "black", width = 0.5) +
coord_polar(theta = "x") +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
title = element_text(size=9.5),
legend.text=element_text(size=7),
legend.title=element_blank(),
legend.position = "bottom",
legend.key.width = unit(0.3, "cm"),
legend.key.height = unit(0.3, "cm"),
axis.text.x = element_text(color = c("red","red",
"black","black","black","black",
"#4A9216","#4A9216","#4A9216",
"#1577BB","#1577BB","#1577BB","#1577BB","#1577BB",
"red","red"))) +
scale_fill_manual(labels = c("Ammonia", "Benzene", "Carbon Black", "Ethylbenzene", "n-Heptane",
"n-Hexane","Hydrogen Sulfide", "Methylpentane", "n-Octane", "PM 2.5",
"Toluene", "1,2,3-Trimethylbenz", "1,2,4-Trimethylbenz",
"1,3,5-Trimethylbenz", "2,2,4-Trimethylpent",
"p-Xylene", "o-Xylene"),
values=c("#F8766D", "#E7851E", "#D09400", "#B2A100", "#89AC00", "#45B500",
"#00BC51", "#00C087", "#00C0B2", "#00BCD6", "#00B3F2", "#29A3FF",
"#9C8DFF", "#D277FF", "#F166E8", "#FF61C7", "#FF689E"), drop = FALSE)
#Print plots
grid.arrange(g1, g2, ncol = 2, widths = c(3,3))
Each Chemical
These plots show the average measurement of each chemical species by wind direction.
#A function that builds a radial plot for a given chemical
chemical_plot = function(chemical, chemical_name){
ggplot(data = subset(c_agg, Chemical == chemical),
aes(x = WindDirCat,
y = Mean,
fill = Chemical)) +
labs(title = chemical_name) +
geom_bar(stat = "identity",
col = "black",
width = 0.5) +
coord_polar(theta = "x") +
theme(axis.title = element_blank(),
legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
title = element_text(size=9.5),
axis.text.x = element_text(color = c("red","red",
"black","black","black","black",
"#4A9216","#4A9216","#4A9216",
"#1577BB","#1577BB","#1577BB","#1577BB","#1577BB",
"red","red")))+
scale_fill_manual(values=c("#F8766D", "#E7851E", "#D09400", "#B2A100", "#89AC00", "#45B500",
"#00BC51", "#00C087", "#00C0B2", "#00BCD6", "#00B3F2", "#29A3FF",
"#9C8DFF", "#D277FF", "#F166E8", "#FF61C7", "#FF689E"), drop = FALSE)}#Aggregating the measurement data by wind direction
c_agg = aggregate(c[,2:18],
by = list(WindDirCat = c$WindDirCat),
function(x) mean(x, na.rm = TRUE))
c_agg = gather(c_agg,
key = "Chemical",
value = "Mean",
-WindDirCat)
c_agg$Chemical = as.factor(c_agg$Chemical)#Plots for each chemical species using the previously defined function
g1 = chemical_plot("Ammonia", "Ammonia")
g2 = chemical_plot("Benzene", "Benzene")
g3 = chemical_plot("Carbon", "Carbon Black")
g4 = chemical_plot("Ebenz", "Ethylbenzene")
g5 = chemical_plot("Hept", "n-Heptane")
g6 = chemical_plot("Hex", "n-Hexane")
g7 = chemical_plot("Hsulf", "Hydrogen Sulfide")
g8 = chemical_plot("Mpent", "Methylpentane")
g9 = chemical_plot("Oct", "n-Octane")
g10 = chemical_plot("PM25", "PM 2.5")
g11 = chemical_plot("Tol", "Toluene")
g12 = chemical_plot("Trimeth123", "1,2,3-Trimethylbenzene")
g13 = chemical_plot("Trimeth124", "1,2,4-Trimethylbenzene")
g14 = chemical_plot("Trimeth135", "1,3,5-Trimethylbenzene")
g15 = chemical_plot("Trimeth224", "2,2,4-Trimethylpentane")
g16 = chemical_plot("Xymp", "p-Xylene")
g17 = chemical_plot("Xyo", "o-Xylene")1
#Page 1
grid.arrange(g1, g2, g3, g4, g5, g6, ncol = 3)
2
#Page 2
grid.arrange(g7, g8, g9, g10, g11, g12, ncol = 3)
3
#Page 3
grid.arrange(g13, g14, g15, g16, g17, ncol = 3)