We also end up using `magrittr` and `tidyr` but only for one function, so you’ll see those with `::` in the code.<\/p>\n
Grabbing the USDA page is pretty straightforward:<\/p>\n
url <- \"http:\/\/www.aphis.usda.gov\/wps\/portal\/aphis\/ourfocus\/animalhealth\/sa_animal_disease_information\/sa_avian_health\/ct_avian_influenza_disease\/!ut\/p\/a1\/lVJbb4IwFP41e1qwFZDLI-oUnGgyswm8kAMUaAaFQNG4X7-ibnEPYtakDz3nO_kupyhAHgoYHGgGnFYMiv4daOFqa8vjKZad5c58wc7mY-Eaa13Z2qoA-AKA7xwL_53fvjpaP_-Gp_Z8jHcK2qMABTHjNc-RD3VO2zCuGCeMhwWNGmhOT7iFsOqaMK3irj2_gNESijAnUPD8tpLQlkBLQsrSqinPJi7tAwX2i4_5tSBgRUfYF_wM9mLqmCbIj2QzxZpMJMUYg6TGkSLBBCaSPEnSJIljXVH0q_kBdw_CO5sXkNnSslV9LQJTDRk7czGumy7GjnYFDOTrCw36XRJTRbt_mlo9VD1FgfuctqrVByA37szNBAPwXOpzR97gPi7tm30gb2AfQkxWVJH4ifvZLavFIsUQrA1JSUOaUV61HHnH43HUtQmMsuqA6vK9NJST9JluNlKwyPr7DT6YvRs!\/?1dmy&urile=wcm%3apath%3a%2Faphis_content_library%2Fsa_our_focus%2Fsa_animal_health%2Fsa_animal_disease_information%2Fsa_avian_health%2Fsa_detections_by_states%2Fct_ai_pacific_flyway\"\r\n\r\n#' read in the data, extract the table and clean up the fields\r\n#' also clean up the column names to since they are fairly nasty\r\n\r\npg <- html(url)<\/pre>\nIf you poke at the source for the page you'll see there are two tables in the code and we only need the first one. Also, if you scan the rendered table on the USDA page by eye you'll see that the column names are horrible for data analysis work and they are also inconsistent in the values used for various columns. Furthermore, there are commas in the flock counts and it would be handy to have the date as an actual date type. We can extract the table we need and clean all that up in a reasonably-sized `dplyr` pipe:<\/p>\n
pg %>%\r\n html_nodes(\"table\") %>%\r\n magrittr::extract2(1) %>%\r\n html_table(header=TRUE) %>%\r\n filter(`Flock size`!=\"pending\") %>%\r\n mutate(Species=str_replace(tolower(Species), \"s$\", \"\"),\r\n `Avian influenza subtype*`=str_replace_all(`Avian influenza subtype*`, \" \", \"\"),\r\n `Flock size`=as.numeric(str_replace_all(`Flock size`, \",\", \"\")),\r\n `Confirmation date`=as.Date(mdy(`Confirmation date`))) %>%\r\n rename(state=State, county=County, flyway=Flyway, flock_type=`Flock type`,\r\n species=Species, subtype=`Avian influenza subtype*`, date=`Confirmation date`,\r\n flock_size=`Flock size`) -> birds<\/pre>\nLet's take a look at what we have:<\/p>\n
glimpse(birds)\r\n\r\n## Observations: 202\r\n## Variables:\r\n## $ state (chr) \"Iowa\", \"Minnesota\", \"Minnesota\", \"Iowa\", \"Minnesota\", \"Iowa\",...\r\n## $ county (chr) \"Sac\", \"Renville\", \"Renville\", \"Hamilton\", \"Kandiyohi\", \"Hamil...\r\n## $ flyway (chr) \"Mississippi\", \"Mississippi\", \"Mississippi\", \"Mississippi\", \"M...\r\n## $ flock_type (chr) \"Commercial\", \"Commercial\", \"Commercial\", \"Commercial\", \"Comme...\r\n## $ species (chr) \"turkey\", \"chicken\", \"turkey\", \"turkey\", \"turkey\", \"turkey\", \"...\r\n## $ subtype (chr) \"EA\/AM-H5N2\", \"EA\/AM-H5N2\", \"EA\/AM-H5N2\", \"EA\/AM-H5N2\", \"EA\/AM...\r\n## $ date (date) 2015-06-04, 2015-06-04, 2015-06-04, 2015-06-04, 2015-06-03, 2...\r\n## $ flock_size (dbl) 42200, 415000, 24800, 19600, 37000, 26200, 17200, 1115700, 159...<\/pre>\nTo make an animated map of cumulative flock totals by week, we'll need to<\/p>\n
- group the `birds` data frame by week and state
\n- calculate the cumulative sums
\n- fill in the gaps where there are missing state\/week combinations
\n- carry the last observations by state\/week forward in this expanded data frame
\n- make breaks for data ranges so we can more intelligently map them to colors<\/p>\n
This ends up being a longer `dplyr` pipe than I usuall like to code (I think very long ones are hard to follow) but it gets the job done and is still pretty readable:<\/p>\n
\r\nbirds %>%\r\n mutate(week=as.numeric(format(birds$date, \"%Y%U\"))) %>%\r\n arrange(week) %>%\r\n group_by(week, state) %>%\r\n tally(flock_size) %>%\r\n group_by(state) %>%\r\n mutate(cum=cumsum(n)) %>%\r\n ungroup %>%\r\n select(week, state, cum) %>%\r\n mutate(week=as.Date(paste(week, 1), \"%Y%U %u\")) %>%\r\n left_join(tidyr::expand(., week, state), .) %>%\r\n group_by(state) %>%\r\n do(na.locf(.)) %>%\r\n mutate(state_abb=state.abb[match(state, state.name)],\r\n cum=as.numeric(ifelse(is.na(cum), 0, cum)),\r\n brks=cut(cum,\r\n breaks=c(0, 200, 50000, 1000000, 10000000, 50000000),\r\n labels=c(\"1-200\", \"201-50K\", \"50k-1m\",\r\n \"1m-10m\", \"10m-50m\"))) -> by_state_and_week<\/pre>\nNow, we perform the standard animation steps:<\/p>\n
- determine where we're going to break the data up
\n- feed that into a loop
\n- partition the data in the loop
\n- render the plot to a file
\n- combine all the individual images into an animation<\/p>\n
For this graphic, I'm doing something a bit extra. The color ranges for the hexbin choropleth go from very light to very dark, so it would be helpful if the titles for the states went from very dark to very light, matching the state colors. The lines that do this check for state breaks that fall in the last two values and appropriately assign `\"black\"` or `\"white\"` as the color.<\/p>\n
\r\ni <- 0\r\n\r\nfor (wk in unique(by_state_and_week$week)) {\r\n\r\n # filter by week\r\n\r\n by_state_and_week %>% filter(week==wk) -> this_wk\r\n\r\n # hack to let us color the state labels in white or black depending on\r\n # the value of the fill\r\n\r\n this_wk %>%\r\n filter(brks %in% c(\"1m-10m\", \"10m-50m\")) %>%\r\n .$state_abb %>%\r\n unique -> white_states\r\n\r\n centers %>%\r\n mutate(txt_col=\"black\") %>%\r\n mutate(txt_col=ifelse(id %in% white_states, \"white\", \"black\")) -> centers\r\n\r\n # setup the plot\r\n\r\n gg <- ggplot()\r\n gg <- gg + geom_map(data=us_map, map=us_map,\r\n aes(x=long, y=lat, map_id=id),\r\n color=\"white\", fill=\"#dddddd\", size=2)\r\n gg <- gg + geom_map(data=this_wk, map=us_map,\r\n aes(fill=brks, map_id=state_abb),\r\n color=\"white\", size=2)\r\n gg <- gg + geom_text(data=centers,\r\n aes(label=id, x=x, y=y, color=txt_col), size=4)\r\n gg <- gg + scale_color_identity()\r\n gg <- gg + scale_fill_brewer(name=\"Combined flock size\\n(all types)\",\r\n palette=\"RdPu\", na.value=\"#dddddd\", drop=FALSE)\r\n gg <- gg + guides(fill=guide_legend(override.aes=list(colour=NA)))\r\n gg <- gg + coord_map()\r\n gg <- gg + labs(x=NULL, y=NULL,\r\n title=sprintf(\"U.S. Avian Flu Total Impact as of %s\\n\", wk))\r\n gg <- gg + theme_bw()\r\n gg <- gg + theme(plot.title=element_text(face=\"bold\", hjust=0, size=24))\r\n gg <- gg + theme(panel.border=element_blank())\r\n gg <- gg + theme(panel.grid=element_blank())\r\n gg <- gg + theme(axis.ticks=element_blank())\r\n gg <- gg + theme(axis.text=element_blank())\r\n gg <- gg + theme(legend.position=\"bottom\")\r\n gg <- gg + theme(legend.direction=\"horizontal\")\r\n gg <- gg + theme(legend.title.align=1)\r\n\r\n # save the image\r\n\r\n # i'm using \"quartz\" here since I'm on a Mac. Use what works for you system to ensure you\r\n # get the best looking output png\r\n\r\n png(sprintf(\"output\/%03d.png\", i), width=800, height=500, type=\"quartz\")\r\n print(gg)\r\n dev.off()\r\n\r\n i <- i + 1\r\n\r\n}<\/pre>\nWe could use one of the R animation packages to actually make the animation, but I know ImageMagick pretty well so I just call it as a `system` command:<\/p>\n
system(\"convert -delay 60 -loop 1 output\/*png output\/avian.gif\")<\/pre>\nAll that results in:<\/p>\n
![\"avian\"](\"https:\/\/i0.wp.com\/rud.is\/b\/wp-content\/uploads\/2015\/06\/avian.gif?ssl=1\")
<\/a><\/p>\n