-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathair_qual_analysis.R
86 lines (68 loc) · 3.96 KB
/
air_qual_analysis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#code for looking at air quality over parks
#variables:
#AOD (aerosols): no unit
#Methane: ppbv
#Nitrogen Dioxide (NO2): 1/square cm
setwd("/Users/hannah/Desktop/SpaceApps2020/")
#read in data from 2020 - obtained from Giovanni
f2020 <- cbind(data.frame(lapply(c("aerosols_2020.csv", "aerosol2_2020.csv", "methane_2020.csv", "no2_2020.csv"), read.csv, header=TRUE, skip = 8)))
f2020 <- f2020[,c(1,2,4,6,8)]
colnames(f2020) <- c("Date", "Aerosols1", "Aerosols2", "Methane", "NO2")
f2020 <- f2020[-which(f2020[,1]=="2/29/20"),] #ignores leap day
rownames(f2020) <- f2020$Date
f2020 <- f2020[,-1]
#read in data from 2019 - obtained from Giovanni
f2019 <- cbind(data.frame(lapply(c("aerosols_2019.csv", "aerosol2_2019.csv", "methane_2019.csv", "no2_2019.csv"), read.csv, header=TRUE, skip = 8)))
f2019 <- f2019[1:146,c(2,4,6,8)]
colnames(f2019) <- c("Aerosols1", "Aerosols2", "Methane", "NO2")
rownames(f2019) <- rownames(f2020)
f2019$Aerosols1[f2019$Aerosols1==-9.999000e+03] <- NA
f2019$Aerosols2[f2019$Aerosols2==-9.999000e+03] <- NA
f2019$Methane[f2019$Methane==-9999.000] <- NA
f2019$NO2[f2019$NO2==-1.267651e+30] <- NA
f2020$Aerosols1[f2020$Aerosols1==-9.999000e+03] <- NA
f2020$Aerosols2[f2020$Aerosols2==-9.999000e+03] <- NA
f2020$Methane[f2020$Methane==-9999.000] <- NA
f2020$NO2[f2020$NO2==-1.27e+30] <- NA
#COMPARISON IN SCATTERPLOTS - ugly!!
#compare aerosols
plot(f2019[,1], type="p", col="#566456", xlab="Date", ylab="Aerosols (unitless)", main="Aerosols over Great Smoky Park", pch=16, cex=0.8)
points(f2019[,2], col="#566456", pch=16, cex=0.8)
points(f2020[,1], col="#9CAEA9", pch=16, cex=0.8)
legend("topright", legend=c("2019", "2020"), col=c("#566456", "#9CAEA9"), lty=1, cex=0.8, lwd=2)
axis(side=1, labels=rownames(f2019))
#compare methane
plot(f2019[,3], type="p", col="#566456", xlab="Date", ylab="Methane (ppbv)", main="Methane over Great Smoky Park", pch=16, cex=0.8)
points(f2020[,3], col="#9CAEA9", pch=16, cex=0.8)
legend("topleft", legend=c("2019", "2020"), col=c("#566456", "#9CAEA9"), lty=1, cex=0.8)
#compare no2
plot(f2019[,4], type="p", col="#566456", xlab="Date", ylab="NO2 (1/square cm)", main="NO2 over Great Smoky Park", pch=16, cex=0.8)
points(f2020[,4], col="#9CAEA9", pch=16, cex=0.8)
legend("topright", legend=c("2019", "2020"), col=c("#566456", "#9CAEA9"), lty=1, cex=0.8)
###########################################################
##Using bar graphs to compare averages in 2019/2020 by month
###########################################################
month <- data.frame(rownames(f2019))
month <- lapply(month, function(y) gsub("(.*)/.*/.*", "\\1", y))
f2019 <- cbind(f2019, month)
f2020 <- cbind(f2020, month)
colnames(f2019) <- c("Aerosols1", "Aerosols2", "Methane", "NO2", "month")
colnames(f2020) <- c("Aerosols1", "Aerosols2", "Methane", "NO2", "month")
aero <- data.frame(matrix(ncol=5, nrow=2))
colnames(aero) <- c("Jan", "Feb", "Mar", "Apr", "May")
rownames(aero) <- c("2019", "2020")
no2 <- data.frame(matrix(ncol=5, nrow=2))
colnames(no2) <- c("Jan", "Feb", "Mar", "Apr", "May")
rownames(no2) <- c("2019", "2020")
for (i in c(1:5)){
temp19 <- f2019[which(f2019$month==i),]
temp20 <- f2020[which(f2020$month==i),]
aero[1,i] <- sum(rbind(temp19[,1], temp19[,2]), na.rm=TRUE)/(nrow(temp19)*2)
aero[2,i] <- sum(rbind(temp20[,1], temp20[,2]), na.rm=TRUE)/(nrow(temp20)*2)
no2[1,i] <- sum(temp19[,4], na.rm=TRUE)/nrow(temp19)
no2[2,i] <- sum(temp20[,4], na.rm=TRUE)/nrow(temp20)
}
barplot(as.matrix(aero), names.arg=colnames(aero), beside=T, main="Average Aerosol Levels in Great Smoky Park", xlab="Month", ylab="Aerosols (no units)", col=c("#566456", "#9CAEA9"))
legend("topleft", inset=0.05,legend=c("2019", "2020"), col=c("#566456", "#9CAEA9"), pch=15, cex=0.8, bty="n", pt.cex=1.8)
barplot(as.matrix(no2), names.arg=colnames(no2), beside=T, main="Average NO2 Levels in Great Smoky Park", xlab="Month", ylab="NO2 (1/square cm)", col=c("#566456", "#9CAEA9"))
legend("topright",legend=c("2019", "2020"), col=c("#566456", "#9CAEA9"), pch=15, cex=0.8, bty="n", pt.cex=1.8)