ExploratoryDataAssignment2/plot4.R at master · ibanmd/ExploratoryDataAssignment2 · GitHub

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
83
84
85
86
87
88
89
90
91
92
## Instructions for plot4.R
## ------------------------------------------------------------------------------------
## Across the United States, how have emissions from coal combustion-related sources
## changed from 1999–2008?

## Download and read in data
download.file(url="https://d396qusza40orc.cloudfront.net/exdata%2Fdata%2FNEI_data.zip",
              destfile="Exploratory2.zip",
              method="curl")
unzip(zipfile="Exploratory2.zip")
NEI <- readRDS("summarySCC_PM25.rds")
SCC <- readRDS("Source_Classification_Code.rds")

### --------------------------------------------------------------------------###
### Prep work
### --------------------------------------------------------------------------###

## Load ggplot2 and data.table
library(ggplot2)
library(data.table)

## Subset for coal combustion
## These were found visually by looking at grep("coal|Coal", names(table(SCC$EI.Sector)), value=TRUE)
## Find the SCC codes which refer to coal combustion, and make a character vector since they are characters in NEI
coalSCC <- as.character(SCC[grep("Coal", SCC$EI.Sector),]$SCC)

## Subset NEI according to the SCC values for coal combustion, and make it a data.table
coal <- NEI[NEI$SCC %in% coalSCC,]
coal <- data.table(coal)

## Sum Emissions by year
plot4 <- coal[,sum(Emissions), by=year]

## Add new column for Emissions measured in tens of thousands of tons (instead of tons as in V1)
plot4[,V2 := round(V1/10000, 2)]

## Find out how many observations there were per year, and make into numeric vector
yearcount <- as.numeric(table(coal$year))

## New column with variable representing "average emissions per observation, in tons"
plot4[,V3 := round(V1/yearcount,2)]


### --------------------------------------------------------------------------###
### Now we make the FOURTH plot and save it
### --------------------------------------------------------------------------###

png("/Users/Mario/Desktop/ExploratoryDataAssignment2/plot4.png", width=600, height=600)

## Set up the first plot, Total Emissions ~ Year
plot(V2 ~ year,
     data=plot4,
     xlab="Year",
     ylab="",
     type="b",
     col="red",
     lwd=2,
     pch=19,
     yaxt="n",
     xaxt="n",
     ylim=c(30,120),
     main="Total and Average Emissions from Coal Combustion by Year (Whole Country)"
)

## Custom axis so that both plots look better and are easier to read
axis(side=2, at=seq(30,120,15))
axis(side=1, at=c(1999,2002,2005,2008))

## Label the first set of points
text(plot4[,year], plot4[,V2], labels=round(plot4[,V2],1), pos=3)

## Add the points for the second plot, Average Coal Emissions per Observation
points(V3 ~ year,
       data=plot4,
       type="b",
       col="blue",
       lwd=2,
       pch=19)

## Add labels to the second set of points
text(plot4[,year], plot4[,V3], labels=round(plot4[,V3],1), pos=3)

## Add a legend to the top right
legend("topright",
       legend=c("Total Coal Combustion Emissions, in Ten Thousands of Tons", "Average Coal Combustion Emissions per Observation, in Tons"),
       col=c("red","blue"),
       lty=1,
       lwd=2)

dev.off()