# load data needed to run all code load("v63i02.RData") ################################################### ### chunk number 1: fig1plot ################################################### # reproduce Figure 1 library(micromap) data("edPov") data("USstates") # bring in education and poverty data for 2010 from American Com. Surey data edpov_ds2 <- merge(edPov, edpov_CDC_ds1, by.x = "state", by.y="State") edpov_ds2 <- edpov_ds2[c(-2, -3)] #dropped ed pov variables from 2000 Census # run create_map_table function to create table for lmplot function # using statistical data and states SpatialPolygonsDataFrame statePolys <- create_map_table(USstates, 'ST') mmplot (stat.data = edpov_ds2, map.data = statePolys, panel.types = c('dot_legend', 'labels', 'dot', 'dot', 'map'), panel.data = list(NA, 'state', 'Poverty_ACS_2006_2010', 'Bachelor_ACS_2006_2010', NA), map.link = c('StateAb','ID'), ord.by='Poverty_ACS_2006_2010', rev.ord = TRUE, grouping = 5, median.row = TRUE, plot.height = 9, plot.width=9, colors=(brewer.pal(6,"Reds"))[2:6], map.color2='lightgray', panel.att = list(list(1, point.type = 20), list(2, header='States', align='left', text.size = .9), list(3, header='Percent Families Living \nBelow Poverty Level', graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size = 1.5, xaxis.ticks = list(5, 7.5, 10, 12.5, 15), xaxis.labels = list(5, 7.5, 10, 12.5, 15), xaxis.title='Percent'), list(4, header='Percent Adults With at\nleast a Bachelors Degree', graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size = 1.5, xaxis.ticks = list(20,30,40,50), xaxis.labels = list(20,30,40, 50), xaxis.title='Percent'), list(5, header='Light Gray Means\nHighlighted Above', inactive.border.color = gray(.7), inactive.border.size = 1))) ################################################### ### chunk number 2: Example code 1: polygon simplification ################################################### # maptools method: library(maptools) Polys_v1 <- thinnedSpatialPoly(txeco,tolerance = 7000, minarea = 1,avoidGEOS = TRUE) plot(Polys_v1) # rgeos method: library(rgeos) Polys_v2 <- gSimplify(txeco, 7000, topologyPreserve = TRUE) Polys_v2 <- SpatialPolygonsDataFrame(Polys_v2, txeco@data) plot(Polys_v2) ################################################### ### chunk number 3: Example code 2 and fig2plot ################################################### wsa.polys <- create_map_table(WSA9, "WSA_9_NM") names(NLA_pH.ds2) NLA_pH.ds2 NLA_pH.ds3 <- NLA_pH.ds2[c(-1:-6,-8,-13)] #remove variables not needed for table 1 str(NLA_pH.ds3) library(plyr) NLA_pH.ds4 <- rename(NLA_pH.ds3, c("Q1.25." = "Q1", "Median.50." = "Median", "Q3.75." = "Q3")) N <- as.integer(NLA_pH.ds4$N) NLA_pH.ds4 <- data.frame(NLA_pH.ds4[c(-6)], N) NLA_pH.ds4 NLA_pH.ds5<-(NLA_pH.ds4[c(-5)]) NLA_pH.ds5 library(xtable) table1 <- xtable(NLA_pH.ds5, label='Table1', caption='National lake assessment pH median estimates and five number summaries of pH for the nine regions.') print(table1, floating.environment='sidewaystable') # example code for draft micromap plot mmplot(stat.data = NLA_pH.ds4, map.data = wsa.polys, panel.types = c("dot_legend", "labels", "box_summary", "dot_cl", "map"), panel.data = list( NA, "Ecoregions", list("Min", "Q1", "Median", "Q3", "Max"), list("Estimate", "LCB95Pct", "UCB95Pct"), NA), ord.by="Estimate", grouping = 3, median.row = FALSE, map.link = c("Ecoregions", "ID")) # reproduce Figure 2 mmplot(stat.data = NLA_pH.ds4, map.data = wsa.polys, panel.types = c('dot_legend', 'labels', 'box_summary', 'dot_cl', 'map'), panel.data = list( NA, 'Ecoregions', list('Min', 'Q1', 'Median', 'Q3', 'Max'), list('Estimate', 'LCB95Pct', 'UCB95Pct'), NA), ord.by='Estimate', grouping = 3, median.row = FALSE, map.link = c("Ecoregions", "ID"), plot.height = 7, plot.width = 7, colors=(brewer.pal(6,"BuPu"))[4:6], map.color2='lightgray', panel.att = list(list(1, point.type = 20, point.size = 2), list(2, header='NARS Reporting\n Regions', align='left', text.size = .8, panel.width = 1.25), list(3, header='Sampled Field\n pH from Lakes', graph.bgcolor=(brewer.pal(5,"Greys"))[1], xaxis.ticks = list(3.0, 5.0, 7.0, 9.0, 11.0), xaxis.labels = list(3.0, 5.0, 7.0, 9.0, 11.0), xaxis.title='pH', graph.bar.size = .4), list(4, header='Estimated Median\n pH with 95% CL', graph.bgcolor=(brewer.pal(5,"Greys"))[1], xaxis.ticks = list(6.5, 7.5, 8.5, 9.5), xaxis.labels = list(6.5, 7.5, 8.5, 9.5), xaxis.title='pH'), list(5, header='Light Gray Means\n Highlighted Above', inactive.border.color = gray(.7), inactive.border.size = 2, panel.width = 2))) ################################################### ### chunk number 4: fig3plot ################################################### # reproduce Figure 3 data1$County.name <- tolower(data1$County.name) #creating myPolyTable containing the simplified Tx counties myPolys<-TXcounty tx.polys <- create_map_table(myPolys, 'ID') tx.polys$ID <- tolower(tx.polys$ID) options(scipen=3) # we'll first log transform our variables before plotting cols = c(3, 4, 5); data1[,cols] = apply(data1[,cols], 2, function(x) log10(x + 1)); # we have to reorder the data first to split it up properly data1 <- data1[rev(order(data1$Total.Leukemia)),] # figure out how to split the data table up nrow(data1)/4 # create 4 objects for each quadrant of the plot with # minor adjustments specifying which lines of data and how to group p1 <- mmplot(stat.data=data1[1:43,], map.data=tx.polys, panel.types=c('map', 'dot_legend', 'labels', 'dot', 'dot','dot'), panel.data=list(NA, NA, 'County.name', 'Total.Leukemia', 'Road.emissions','Airports.emissions'), ord.by='Total.Leukemia', grouping=c(rep(6,6),7), median.row=F, map.link=c('County.name','ID'), map.color2='lightgray', rev.ord=TRUE, panel.att=list(list(2,panel.header.size=.6, point.size=.5), list(3, header='Counties', panel.width=.8,panel.header.size=.6, align='left', text.size=.3), list(4, header='Total Leukemia',panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(5, header='Road Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(6, header='Airport Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(1,inactive.border.color=gray(.7), inactive.border.size=.2, panel.header.size=.6,map.all=T, nodata.fill='white',nodata.border.color=gray(.7),nodata.border.size=.5, active.border.size=.2,active.border.color='black', panel.width=.9))) p1 <- lapply(p1, function(pp) pp + theme(axis.text.x = element_text(size=5))) p2 <- mmplot(stat.data=data1[1:43+43,], map.data=tx.polys, panel.types=c('map', 'dot_legend', 'labels', 'dot', 'dot','dot'), panel.data=list(NA, NA, 'County.name', 'Total.Leukemia', 'Road.emissions','Airports.emissions'), ord.by='Total.Leukemia', grouping=c(rep(6,6),7), median.row=F, map.link=c('County.name','ID'), map.color2='lightgray', rev.ord=TRUE, panel.att=list(list(2,panel.header.size=.6, point.size=.5), list(3, header='Counties', panel.width=.8,panel.header.size=.6, align='left', text.size=.3), list(4, header='Total Leukemia',panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(5, header='Road Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(6, header='Airport Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(1,inactive.border.color=gray(.7), inactive.border.size=.2, panel.header.size=.6,map.all=T, nodata.fill='white',nodata.border.color=gray(.7),nodata.border.size=.5, active.border.size=.2,active.border.color='black', panel.width=.9))) p2 <- lapply(p2, function(pp) pp + theme(axis.text.x = element_text(size=5))) p3 <- mmplot(stat.data=data1[1:42+86,], map.data=tx.polys, panel.types=c('map', 'dot_legend', 'labels', 'dot', 'dot','dot'), panel.data=list(NA, NA, 'County.name', 'Total.Leukemia', 'Road.emissions','Airports.emissions'), ord.by='Total.Leukemia', grouping=c(rep(6,6),7), median.row=F, map.link=c('County.name','ID'), map.color2='lightgray', rev.ord=TRUE, panel.att=list(list(2,panel.header.size=.6, point.size=.5), list(3, header='Counties', panel.width=.8,panel.header.size=.6, align='left', text.size=.3), list(4, header='Total Leukemia',panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(5, header='Road Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(6, header='Airport Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(1,inactive.border.color=gray(.7), inactive.border.size=.2, panel.header.size=.6,map.all=T, nodata.fill='white',nodata.border.color=gray(.7),nodata.border.size=.5, active.border.size=.2,active.border.color='black', panel.width=.9))) p3 <- lapply(p3, function(pp) pp + theme(axis.text.x = element_text(size=5))) p4 <- mmplot(stat.data=data1[1:42+128,], map.data=tx.polys, panel.types=c('map', 'dot_legend', 'labels', 'dot', 'dot','dot'), panel.data=list(NA, NA, 'County.name', 'Total.Leukemia', 'Road.emissions','Airports.emissions'), ord.by='Total.Leukemia', grouping=c(rep(6,6),7), median.row=F, map.link=c('County.name','ID'), map.color2='lightgray', rev.ord=TRUE, panel.att=list(list(2,panel.header.size=.6, point.size=.5), list(3, header='Counties', panel.width=.8,panel.header.size=.6, align='left', text.size=.3), list(4, header='Total Leukemia',panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(5, header='Road Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(6, header='Airport Emissions', panel.header.size=.6, graph.bgcolor=(brewer.pal(5,"Greys"))[1], point.size=.7, xaxis.ticks=c(0,1.5,3),xaxis.labels=c(0,1.5,3)), list(1,inactive.border.color=gray(.7), inactive.border.size=.2, panel.header.size=.6,map.all=T, nodata.fill='white',nodata.border.color=gray(.7),nodata.border.size=.5, active.border.size=.2,active.border.color='black', panel.width=.9))) lwidth <- p4$layout$widths p4 <- lapply(p4, function(pp) pp + theme(axis.text.x = element_text(size=5))) ## 6 panels: map, dot, names, 3 dotplots nPanels <- 6 ## There are alist of widths returned with every object ## but the should be the same for all 4 plots ## note: there are blank panels between every ggplot panel to ## allow spacing adjustments by the user. ie length(lwidth) = npanel*2+1 ## sets layout according to specified widths in function arguments lmLayout <- grid.layout(nrow = 2, ncol = length(lwidth)*2, widths = unit(lwidth, rep("null", length(lwidth))), heights = unit(rep(1, length(lwidth)), rep("null", length(lwidth)))) subplot <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) ## subplot is specifying which cell in the layout to populate ## i.e. subplot(row number, column number) ## again, panel number has to account for small blank ## spacing columns, hence the p*2 (where p is the panel number) grid.newpage() pushViewport(viewport(layout = lmLayout)) suppressWarnings(for(p in 1:nPanels) print(p1[[p]], vp = subplot(1, p*2))) suppressWarnings(for(p in 1:nPanels) print(p2[[p]], vp = subplot(1, nPanels*2 + 1 + p*2))) suppressWarnings(for(p in 1:nPanels) print(p3[[p]], vp = subplot(2, p*2))) suppressWarnings(for(p in 1:nPanels) print(p4[[p]], vp = subplot(2, nPanels*2 + 1 + p*2))) ################################################### ### chunk number 5: Example code 3 and fig4plot ################################################### #reproduce Figure 4 head(lake.data) names(lake.data) lake.data1<-lake.data[c(-1, -3, -9:-16)] #remove variables not needed for table2 lake.data1 library(xtable) table2<-xtable(lake.data1, label='Table2', caption='National lake assessment water chemistry estimates by disturbance categories.') print(table2, floating.environment='sidewaystable') wsa.polys <- create_map_table(WSA9) national.polys <-subset(wsa.polys, hole==0 & plug==0) #create polygon for entire US national.polys <- transform(national.polys, ID='National', region=10, poly=region*1000 + poly) wsa.polys <- rbind(wsa.polys, national.polys) mmgroupedplot(stat.data= lake.data, map.data= wsa.polys, panel.types = c('map', 'labels', 'bar_cl', 'bar_cl', 'bar_cl'), panel.data = list(NA,'Category', list('PTL.Estimate.P','PTL.LCB95Pct.P','PTL.UCB95Pct.P'), list('NTL.Estimate.P','NTL.LCB95Pct.P','NTL.UCB95Pct.P'), list('TURB.Estimate.P','TURB.LCB95Pct.P','TURB.UCB95Pct.P')), grp.by='Subpopulation', cat='Category', colors = brewer.pal(7, "RdYlBu")[c(6,4,1)], map.link = c('Subpopulation', 'ID'), map.color='lightgreen', plot.grp.spacing = 2, plot.width = 6, plot.height = 10, panel.att = list(list(1, header='Region', panel.width = 2), list(2, header='Category', text.size = .6, panel.width = 2), list(3, header='Total \n Phosphorus', graph.bgcolor='lightgray', xaxis.title='percent', xaxis.ticks = c(0,20,40,60,80), xaxis.labels = c(0,20,40,60,80)), list(4, header='Total \n Nitrogen', graph.bgcolor='lightgray', xaxis.title='percent', xaxis.ticks = c(0,20,40,60,80), xaxis.labels = c(0,20,40,60,80)), list(5, header='Turbidity', graph.bgcolor='lightgray', xaxis.title='percent', xaxis.ticks = c(0,20,40,60,80), xaxis.labels = c(0,20,40,60,80))))