Mapping PM 2.5 Emissions (Part 2)
Ian Maddaus
This project will map annual PM 2.5 emissions divided by square miles and population over four years, 1999, 2002, 2005, and 2008.
Part 1: here.
All the data files that I use are available here:
This data originally came from here: https://www.census.gov/popest/data/historical/index.html
Loading Data
Picking up from the earlier project, start by reading two new data files into R using fread. Because the census is conducted only once every 10 years the data for these years are estimates. The data is split into two separate files, the 98-99 data was in a text file which I processed into a csv before bringing it into R and the 2000-2008 data was in a separate csv file available from the census.gov website.
CensusEstimate9899 <- fread('9899CensusEstimate.csv', header = TRUE)
CensusEstimate0008 <- fread('CO-EST2008-ALLDATA.csv', header = TRUE)
head(CensusEstimate9899)## Fips State / County 7/1/98 Estimate 7/1/99 Estimate
## 1: NA United States 270,248,003 272,690,813
## 2: 1 Alabama 4,351,037 4,369,862
## 3: 1001 Autauga County, AL 42,193 43,140
## 4: 1003 Baldwin County, AL 132,857 135,820
## 5: 1005 Barbour County, AL 26,936 26,726
## 6: 1007 Bibb County, AL 18,987 19,601
## Numeric Population Change 98-99 Percent Pop Change
## 1: 2,442,810 0.9
## 2: 18,825 0.4
## 3: 947 2.2
## 4: 2,963 2.2
## 5: -210 -0.8
## 6: 614 3.2head(CensusEstimate0008)## SUMLEV REGION DIVISION STATE COUNTY STNAME CTYNAME
## 1: 40 3 6 1 0 Alabama Alabama
## 2: 50 3 6 1 1 Alabama Autauga County
## 3: 50 3 6 1 3 Alabama Baldwin County
## 4: 50 3 6 1 5 Alabama Barbour County
## 5: 50 3 6 1 7 Alabama Bibb County
## 6: 50 3 6 1 9 Alabama Blount County
## CENSUS2000POP ESTIMATESBASE2000 POPESTIMATE2000 POPESTIMATE2001
## 1: 4447100 4447355 4451687 4462832
## 2: 43671 43671 43874 44437
## 3: 140415 140415 141354 144936
## 4: 29038 29038 29032 29210
## 5: 20826 19889 19937 20942
## 6: 51024 51023 51177 51979
## POPESTIMATE2002 POPESTIMATE2003 POPESTIMATE2004 POPESTIMATE2005
## 1: 4469906 4486598 4506574 4537299
## 2: 45160 45766 46941 47882
## 3: 148031 151502 156276 162149
## 4: 29266 29441 29402 29372
## 5: 20983 21038 21209 21375
## 6: 52737 53660 54381 54915
## POPESTIMATE2006 POPESTIMATE2007 POPESTIMATE2008 NPOPCHG_2000
## 1: 4587564 4626595 4661900 4332
## 2: 49039 49830 50364 203
## 3: 168154 171748 174439 939
## 4: 29420 29460 29309 -6
## 5: 21307 21483 21629 48
## 6: 55708 56502 57441 154
## NPOPCHG_2001 NPOPCHG_2002 NPOPCHG_2003 NPOPCHG_2004 NPOPCHG_2005
## 1: 11145 7074 16692 19976 30725
## 2: 563 723 606 1175 941
## 3: 3582 3095 3471 4774 5873
## 4: 178 56 175 -39 -30
## 5: 1005 41 55 171 166
## 6: 802 758 923 721 534
## NPOPCHG_2006 NPOPCHG_2007 NPOPCHG_2008 BIRTHS2000 BIRTHS2001 BIRTHS2002
## 1: 50265 39031 35305 15058 62358 59529
## 2: 1157 791 534 142 635 593
## 3: 6005 3594 2691 428 1890 1753
## 4: 48 40 -151 84 429 377
## 5: -68 176 146 67 272 296
## 6: 793 794 939 165 687 677
## BIRTHS2003 BIRTHS2004 BIRTHS2005 BIRTHS2006 BIRTHS2007 BIRTHS2008
## 1: 58571 59547 60076 61332 63046 62940
## 2: 613 635 674 646 661 718
## 3: 1795 1846 2025 2199 2199 2267
## 4: 360 313 369 363 396 329
## 5: 247 275 251 268 252 257
## 6: 661 689 685 675 724 729
## DEATHS2000 DEATHS2001 DEATHS2002 DEATHS2003 DEATHS2004 DEATHS2005
## 1: 10564 45163 45754 45883 46913 46549
## 2: 79 354 367 401 396 407
## 3: 342 1340 1446 1472 1553 1529
## 4: 70 323 334 296 284 301
## 5: 59 193 214 221 235 219
## 6: 87 480 511 556 514 542
## DEATHS2006 DEATHS2007 DEATHS2008 NATURALINC2000 NATURALINC2001
## 1: 47004 46516 47057 4494 17195
## 2: 455 445 448 63 281
## 3: 1630 1653 1672 86 550
## 4: 308 299 299 14 106
## 5: 233 217 218 8 79
## 6: 546 534 554 78 207
## NATURALINC2002 NATURALINC2003 NATURALINC2004 NATURALINC2005
## 1: 13775 12688 12634 13527
## 2: 226 212 239 267
## 3: 307 323 293 496
## 4: 43 64 29 68
## 5: 82 26 40 32
## 6: 166 105 175 143
## NATURALINC2006 NATURALINC2007 NATURALINC2008 INTERNATIONALMIG2000
## 1: 14328 16530 15883 1270
## 2: 191 216 270 6
## 3: 569 546 595 41
## 4: 55 97 30 8
## 5: 35 35 39 1
## 6: 129 190 175 19
## INTERNATIONALMIG2001 INTERNATIONALMIG2002 INTERNATIONALMIG2003
## 1: 4806 3871 2080
## 2: 5 -24 -72
## 3: 179 163 109
## 4: 44 41 33
## 5: 6 4 2
## 6: 97 94 75
## INTERNATIONALMIG2004 INTERNATIONALMIG2005 INTERNATIONALMIG2006
## 1: 4483 3558 4286
## 2: 39 -5 22
## 3: 138 129 139
## 4: 35 33 35
## 5: 7 5 7
## 6: 79 80 83
## INTERNATIONALMIG2007 INTERNATIONALMIG2008 DOMESTICMIG2000
## 1: 3366 3460 -1474
## 2: 0 5 133
## 3: 117 119 837
## 4: 30 30 -29
## 5: 5 5 38
## 6: 74 74 60
## DOMESTICMIG2001 DOMESTICMIG2002 DOMESTICMIG2003 DOMESTICMIG2004
## 1: -8850 -8616 4821 5123
## 2: 318 549 499 956
## 3: 2818 2595 2931 4392
## 4: -205 -212 -92 -290
## 5: 918 -24 54 141
## 6: 536 535 767 517
## DOMESTICMIG2005 DOMESTICMIG2006 DOMESTICMIG2007 DOMESTICMIG2008
## 1: 16248 32945 18496 15118
## 2: 723 970 561 245
## 3: 5185 5304 2932 2055
## 4: -318 -229 -266 -203
## 5: 148 -101 133 92
## 6: 362 609 509 668
## NETMIG2000 NETMIG2001 NETMIG2002 NETMIG2003 NETMIG2004 NETMIG2005
## 1: -204 -4044 -4745 6901 9606 19806
## 2: 139 323 525 427 995 718
## 3: 878 2997 2758 3040 4530 5314
## 4: -21 -161 -171 -59 -255 -285
## 5: 39 924 -20 56 148 153
## 6: 79 633 629 842 596 442
## NETMIG2006 NETMIG2007 NETMIG2008 RESIDUAL2000 RESIDUAL2001 RESIDUAL2002
## 1: 37231 21862 18578 42 -2006 -1956
## 2: 992 561 250 1 -41 -28
## 3: 5443 3049 2174 -25 35 30
## 4: -194 -236 -173 1 233 184
## 5: -94 138 97 1 2 -21
## 6: 692 583 742 -3 -38 -37
## RESIDUAL2003 RESIDUAL2004 RESIDUAL2005 RESIDUAL2006 RESIDUAL2007
## 1: -2897 -2264 -2608 -1294 639
## 2: -33 -59 -44 -26 14
## 3: 108 -49 63 -7 -1
## 4: 170 187 187 187 179
## 5: -27 -17 -19 -9 3
## 6: -24 -50 -51 -28 21
## RESIDUAL2008 GQESTIMATESBASE2000 GQESTIMATES2000 GQESTIMATES2001
## 1: 844 114720 114720 116446
## 2: 14 260 260 260
## 3: -78 2274 2274 2331
## 4: -8 2718 2718 2701
## 5: 10 302 302 1198
## 6: 22 610 610 610
## GQESTIMATES2002 GQESTIMATES2003 GQESTIMATES2004 GQESTIMATES2005
## 1: 115651 118109 116744 117725
## 2: 260 260 260 260
## 3: 2335 2359 2198 2186
## 4: 2893 3149 3162 3162
## 5: 1234 1217 1297 1299
## 6: 610 610 610 610
## GQESTIMATES2006 GQESTIMATES2007 GQESTIMATES2008 RBIRTH2001 RBIRTH2002
## 1: 119152 120244 120244 13.99021 13.32828
## 2: 260 260 260 14.38099 13.23705
## 3: 2166 2333 2333 13.20340 11.96722
## 4: 3167 3148 3148 14.73164 12.89418
## 5: 1300 1297 1297 13.30757 14.12045
## 6: 610 610 610 13.31963 12.93021
## RBIRTH2003 RBIRTH2004 RBIRTH2005 RBIRTH2006 RBIRTH2007 RBIRTH2008
## 1: 13.07899 13.24271 13.28546 13.44283 13.68459 13.55225
## 2: 13.48349 13.69907 14.21596 13.33044 13.37123 14.33220
## 3: 11.98532 11.99566 12.71885 13.31505 12.93902 13.09697
## 4: 12.26430 10.63848 12.55657 12.34862 13.45109 11.19638
## 5: 11.75603 13.01868 11.78847 12.55799 11.77845 11.92243
## 6: 12.42516 12.75442 12.53477 12.20361 12.90438 12.79587
## RDEATH2001 RDEATH2002 RDEATH2003 RDEATH2004 RDEATH2005 RDEATH2006
## 1: 10.132459 10.244116 10.245739 10.433026 10.294041 10.302401
## 2: 8.017121 8.192239 8.820359 8.543044 8.584415 9.389090
## 3: 9.361137 9.871419 9.828633 10.091689 9.603517 9.869726
## 4: 11.091652 11.423490 10.083976 9.652805 10.242624 10.477616
## 5: 9.442501 10.208706 10.518550 11.125050 10.285553 10.917951
## 6: 9.306293 9.759731 10.451423 9.514906 9.918021 9.871365
## RDEATH2007 RDEATH2008 RNATURALINC2001 RNATURALINC2002 RNATURALINC2003
## 1: 10.096635 10.132320 3.857752 3.084161 2.833248
## 2: 9.001810 8.942651 6.363873 5.044812 4.663133
## 3: 9.726333 9.659519 3.842258 2.095799 2.156691
## 4: 10.156250 10.175433 3.639985 1.470689 2.180319
## 5: 10.142557 10.113194 3.865065 3.911747 1.237476
## 6: 9.517868 9.724160 4.013339 3.170480 1.973740
## RNATURALINC2004 RNATURALINC2005 RNATURALINC2006 RNATURALINC2007
## 1: 2.8096872 2.991418 3.140431 3.587956
## 2: 5.1560292 5.631545 3.941354 4.369418
## 3: 1.9039697 3.115333 3.445321 3.212691
## 4: 0.9856737 2.313948 1.871003 3.294837
## 5: 1.8936256 1.502912 1.640036 1.635896
## 6: 3.2395109 2.616747 2.332246 3.386507
## RNATURALINC2008 RINTERNATIONALMIG2001 RINTERNATIONALMIG2002
## 1: 3.419930 1.0782410 0.8666995
## 2: 5.389544 0.1132362 -0.5357322
## 3: 3.437449 1.2504803 1.1127533
## 4: 1.020946 1.5109371 1.4022847
## 5: 1.809241 0.2935493 0.1908169
## 6: 3.071711 1.8806468 1.7953321
## RINTERNATIONALMIG2003 RINTERNATIONALMIG2004 RINTERNATIONALMIG2005
## 1: 0.4644669 0.9969786 0.7868310
## 2: -1.5837054 0.8413604 -0.1054596
## 3: 0.7277996 0.8967503 0.8102379
## 4: 1.1242271 1.1896062 1.1229455
## 5: 0.0951905 0.3313845 0.2348300
## 6: 1.4098142 1.4624078 1.4639145
## RINTERNATIONALMIG2006 RINTERNATIONALMIG2007 RINTERNATIONALMIG2008
## 1: 0.9394114 0.7306147 0.74500767
## 2: 0.4539780 0.0000000 0.09980638
## 3: 0.8416515 0.6884337 0.68748971
## 4: 1.1906382 1.0190217 1.02094642
## 5: 0.3280071 0.2336995 0.23195398
## 6: 1.5005921 1.3189555 1.29889506
## RDOMESTICMIG2001 RDOMESTICMIG2002 RDOMESTICMIG2003 RDOMESTICMIG2004
## 1: -1.985525 -1.929084 1.076536 1.139309
## 2: 7.201821 12.254875 10.975958 20.624117
## 3: 19.686332 17.715306 19.570465 28.540052
## 4: -7.039593 -7.250838 -3.134209 -9.856737
## 5: 44.913036 -1.144902 2.570143 6.675030
## 6: 10.392028 10.218114 14.417700 9.570441
## RDOMESTICMIG2005 RDOMESTICMIG2006 RDOMESTICMIG2007 RDOMESTICMIG2008
## 1: 3.593151 7.220930 4.014691 3.255210
## 2: 15.249465 20.016302 11.348350 4.890512
## 3: 32.566538 32.115966 17.252031 11.872196
## 4: -10.821111 -7.790176 -9.035326 -6.908404
## 5: 6.950967 -4.732674 6.216406 4.267953
## 6: 6.624213 11.010369 9.072275 11.725161
## RNETMIG2001 RNETMIG2002 RNETMIG2003 RNETMIG2004 RNETMIG2005 RNETMIG2006
## 1: -0.9072839 -1.0623842 1.541003 2.136287 4.379982 8.160342
## 2: 7.3150570 11.7191424 9.392253 21.465477 15.144005 20.470280
## 3: 20.9368123 18.8280591 20.298264 29.436802 33.376776 32.957618
## 4: -5.5286563 -5.8485533 -2.009982 -8.667131 -9.698166 -6.599537
## 5: 45.2065853 -0.9540847 2.665334 7.006415 7.185797 -4.404667
## 6: 12.2726744 12.0134459 15.827514 11.032849 8.088128 12.510961
## RNETMIG2007 RNETMIG2008
## 1: 4.745306 4.000217
## 2: 11.348350 4.990319
## 3: 17.940465 12.559686
## 4: -8.016304 -5.887458
## 5: 6.450105 4.499907
## 6: 10.391231 13.024056Format the data
The population data in the 98 - 99 data set is character data and has commas as a thousands separator. Just using the as.numeric function will return NA’s because it can’t handle the commas. Remove the commas using gsub then convert to numeric data.
CensusEstimate9899$`7/1/99 Estimate`<- as.numeric(gsub(",","",CensusEstimate9899$`7/1/99 Estimate`))The FIPS codes in the 98-99 population data are numeric data which sometimes is only 4 digits for codes that start with a 0. Convert to five digit integers using sprintf and as.integer.
CensusEstimate9899$Fips <- sprintf("%05d", as.integer(CensusEstimate9899$Fips))
head(CensusEstimate9899)## Fips State / County 7/1/98 Estimate 7/1/99 Estimate
## 1: 000NA United States 270,248,003 272690813
## 2: 00001 Alabama 4,351,037 4369862
## 3: 01001 Autauga County, AL 42,193 43140
## 4: 01003 Baldwin County, AL 132,857 135820
## 5: 01005 Barbour County, AL 26,936 26726
## 6: 01007 Bibb County, AL 18,987 19601
## Numeric Population Change 98-99 Percent Pop Change
## 1: 2,442,810 0.9
## 2: 18,825 0.4
## 3: 947 2.2
## 4: 2,963 2.2
## 5: -210 -0.8
## 6: 614 3.2Now we take the Emissions / SqMile data and divide it by the population of each county. Similar to the divideEmissionsBySqMiles function, this for loop finds the emissions / SqMile data and divides it by the population of each county.
vec <- numeric()
for (i in NinetyNine$fips){
vec <- append(vec, NinetyNine[NinetyNine$fips == i, emissionsSqMiles] / as.numeric(CensusEstimate9899[CensusEstimate9899$Fips == i, '7/1/99 Estimate', with=FALSE]))
}Cbind the results of the function to the NinetyNine data.table and multiply by 100,000 to get tons per 100,000 people.
NinetyNine <- cbind(NinetyNine, emissionsSqMilePerCapita = vec * 100000)
rm(vec,i)
head(NinetyNine)## fips year Emissions emissionsSqMiles emissionsSqMilePerCapita
## 1: 09001 1999 5205.503 8.330844 0.9901946
## 2: 09003 1999 5791.588 7.878115 0.9495468
## 3: 09005 1999 2168.551 2.355625 1.2914683
## 4: 09007 1999 1905.879 5.161123 3.4075593
## 5: 09009 1999 4377.132 7.239105 0.9126364
## 6: 09011 1999 4123.281 6.201430 2.5204044Get FIPS code in 2000 - 2008 data. The fips code in this data set is split into two columns, State and County. States codes that are lower than 10 are just one digit, and county codes that are less than 100 are two or one digit. Use sprintf and as.integer to make them the correct two and three digit lengths.
CensusEstimate0008$STATE <- sprintf("%02d", as.integer(CensusEstimate0008$STATE))
CensusEstimate0008$COUNTY <- sprintf("%03d", as.integer(CensusEstimate0008$COUNTY)) Then concatenate them to create the full five digit Fips code.
CensusEstimate0008$FIPS <- paste(CensusEstimate0008$STATE, CensusEstimate0008$COUNTY, sep="")Create a function similar to the earlier divideEmissionsBySqMiles function that divides the Emissions per square mile by the population of each county.
divideEmissionsPerSqMileByPopulation <- function (emissionsData, populationData, columnNumber) {
vec <- numeric()
for (i in emissionsData$fips){
vec <- append(vec, as.numeric(emissionsData[emissionsData$fips == i, 4, with=FALSE]) / as.numeric(populationData[populationData$FIPS == i, columnNumber, with=FALSE]))
}
return(vec)
}Run the function and multiply the result by 100,000 to get tons of emissions per 100,000 people.
ZeroTwo <- cbind(ZeroTwo, emissionsSqMilePerCapita = divideEmissionsPerSqMileByPopulation(ZeroTwo,CensusEstimate0008, 12) * 100000)
ZeroFive <- cbind(ZeroFive, emissionsSqMilePerCapita = divideEmissionsPerSqMileByPopulation(ZeroFive,CensusEstimate0008, 15) * 100000)
ZeroEight <- cbind(ZeroEight, emissionsSqMilePerCapita = divideEmissionsPerSqMileByPopulation(ZeroEight,CensusEstimate0008, 18) * 100000)Split the emissions into deciles
DecileEmissionsCapitaTons <- quantile(c(NinetyNine$emissionsSqMilePerCapita, ZeroEight$emissionsSqMilePerCapita, ZeroFive$emissionsSqMilePerCapita, ZeroTwo$emissionsSqMilePerCapita), probs=seq(0,1, by=0.1), na.rm = TRUE)
DecileEmissionsCapitaTons## 0% 10% 20% 30% 40%
## 3.698330e-02 1.344701e+00 2.236449e+00 3.117958e+00 4.166124e+00
## 50% 60% 70% 80% 90%
## 5.485297e+00 7.275185e+00 9.730774e+00 1.382390e+01 2.327746e+01
## 100%
## 2.199207e+04Mapping the Data
Create color buckets
NinetyNine$colorbucketsTonsEmissions <- as.numeric(cut(NinetyNine$emissionsSqMilePerCapita, DecileEmissionsCapitaTons))
ZeroTwo$colorbucketsTonsEmissions <- as.numeric(cut(ZeroTwo$emissionsSqMilePerCapita, DecileEmissionsCapitaTons))
ZeroFive$colorbucketsTonsEmissions <- as.numeric(cut(ZeroFive$emissionsSqMilePerCapita, DecileEmissionsCapitaTons))
ZeroEight$colorbucketsTonsEmissions <- as.numeric(cut(ZeroEight$emissionsSqMilePerCapita, DecileEmissionsCapitaTons))Then create vectors for each year that will tell R which color bucket to use for each county.
NinetyNineColorsMatched <- NinetyNine$colorbucketsTonsEmissions[match(cnty.fips, as.numeric(NinetyNine$fips))]
ZeroTwoColorsMatched <- ZeroTwo$colorbucketsTonsEmissions[match(cnty.fips, as.numeric(ZeroTwo$fips))]
ZeroFiveColorsMatched <- ZeroFive$colorbucketsTonsEmissions[match(cnty.fips, as.numeric(ZeroFive$fips))]
ZeroEightColorsMatched <- ZeroEight$colorbucketsTonsEmissions[match(cnty.fips, as.numeric(ZeroEight$fips))]## Warning in match(cnty.fips, as.numeric(ZeroEight$fips)): NAs introduced by
## coercionCreate the legend for each map.
EmissionsPerCapitaLegend.txt <- c("0-1.34", "1.34-2.34", "2.34-3.12", "3.12-4.17", "4.17-5.49", "5.49-7.28", "7.28-9.73", "9.73-13.82", "13.82-23.28", ">23.28")Create a color palette. I ended up using one I found in http://colorbrewer2.org https://gka.github.io/palettes/
Emissions_palette4 <- c("#ffe4e1", "#ffc9bd", "#ffad98", "#ff8f73", "#ff6d4e", "#ff3d21", "#ee0001", "#cb0002", "#ab0002", "#8b0000")Print the emissions maps.
1999 Emissions Map
#1999 Emissions Map
png (file = 'PM25Emissions1999PerCapita.png', width = 800, height = 800, pointsize = 12)
map("county", col = Emissions_palette4[NinetyNineColorsMatched], fill = TRUE, resolution = 0, lty = 0, projection = "polyconic")
map("state", col = "white", fill = FALSE, add = TRUE, lty = 1, lwd = 0.2, projection="polyconic")
title("PM2.5 Emissions in Pounds / SQ Mile / Per Capita / Year - 1999")
legend("top", fill = (Emissions_palette4), legend = (EmissionsPerCapitaLegend.txt), ncol = 2)
dev.off() ## quartz_off_screen
## 2
2002 Emissions Map
#2002 Emissions Map
png (file = 'PM25Emissions2002PerCapita.png', width = 800, height = 800, pointsize = 12)
map("county", col = Emissions_palette4[ZeroTwoColorsMatched], fill = TRUE, resolution = 0, lty = 0, projection = "polyconic")
map("state", col = "white", fill = FALSE, add = TRUE, lty = 1, lwd = 0.2, projection="polyconic")
title("PM2.5 Emissions in Pounds / SQ Mile / Per Capita / Year - 2002")
legend("top", fill = (Emissions_palette4), legend = (EmissionsPerCapitaLegend.txt), ncol = 2)
dev.off() ## quartz_off_screen
## 2
2005 Emissions Map
#2005 Emissions Map
png (file = 'PM25Emissions2005PerCapita.png', width = 800, height = 800, pointsize = 12)
map("county", col = Emissions_palette4[ZeroFiveColorsMatched], fill = TRUE, resolution = 0, lty = 0, projection = "polyconic")
map("state", col = "white", fill = FALSE, add = TRUE, lty = 1, lwd = 0.2, projection="polyconic")
title("PM2.5 Emissions in Pounds / SQ Mile / Per Capita / Year - 2005")
legend("top", fill = (Emissions_palette4), legend = (EmissionsPerCapitaLegend.txt), ncol = 2)
dev.off() ## quartz_off_screen
## 2
2008 Emissions Map
#2008 Emissions Map
png (file = 'PM25Emissions2008PerCapita.png', width = 800, height = 800, pointsize = 12)
map("county", col = Emissions_palette4[ZeroEightColorsMatched], fill = TRUE, resolution = 0, lty = 0, projection = "polyconic")
map("state", col = "white", fill = FALSE, add = TRUE, lty = 1, lwd = 0.2, projection="polyconic")
title("PM2.5 Emissions in Pounds / SQ Mile / Per Capita / Year - 2008")
legend("top", fill = (Emissions_palette4), legend = (EmissionsPerCapitaLegend.txt), ncol = 2)
dev.off() ## quartz_off_screen
## 2