require(readr); require(dplyr)
require(ggplot2); require(ggpubr)
require(sf); require(leaflet)
#reads the csv file
concerts <- read_csv("TheInternetConcerts.csv")
table(concerts$Year)
## 
## 2016 2017 2018 
##   16   20   41
#Filters concerts by year and counts how many concerts occurred in each state
count2016 <- concerts %>% filter(Year == 2016) %>% count(State, name = "concerts") %>% arrange(desc(concerts))
count2017 <- concerts %>% filter(Year == 2017) %>% count(State, name = "concerts") %>% arrange(desc(concerts))
count2018 <- concerts %>% filter(Year == 2018) %>% count(State, name = "concerts") %>% arrange(desc(concerts))
countTotal <- concerts %>% count(State, name = "concerts") %>% arrange(desc(concerts))

head(count2016)
## # A tibble: 6 × 2
##   State                concerts
##   <chr>                   <int>
## 1 California                  3
## 2 New York                    2
## 3 Tennessee                   2
## 4 Washington                  2
## 5 Arizona                     1
## 6 District of Columbia        1
head(countTotal)
## # A tibble: 6 × 2
##   State      concerts
##   <chr>         <int>
## 1 California       24
## 2 New York          7
## 3 Maryland          6
## 4 Texas             6
## 5 Georgia           5
## 6 Illinois          4
#reads shapefile, excluded alaska, hawaii, and puerto rico
US <- st_read("cb_2018_us_state_20m.shp")
## Reading layer `cb_2018_us_state_20m' from data source 
##   `C:\Users\iiStr\Documents\R course\Project\cb_2018_us_state_20m.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS:  NAD83
#print(US [, c("NAME", "STATEFP")], n = Inf)
US <- US[US$STATEFP !="02" & US$STATEFP != "15" & US$STATEFP != "72", ]
US <- st_transform(US, crs = 4326)
plot(st_geometry(US))

#combines each year's concert count to the shapefile
map2016 <- merge(US, count2016, by.x = "NAME", by.y = "State", all.x = TRUE)
map2017 <- merge(US, count2017, by.x = "NAME", by.y = "State", all.x = TRUE)
map2018 <- merge(US, count2018, by.x = "NAME", by.y = "State", all.x = TRUE)
mapTotal <- merge(US, countTotal, by.x = "NAME", by.y = "State", all.x = TRUE)

#Replaces NA with 0
map2016[is.na(map2016$'concerts' ), "concerts"] <- 0
map2017[is.na(map2017$'concerts' ), "concerts"] <- 0
map2018[is.na(map2018$'concerts' ), "concerts"] <- 0
mapTotal[is.na(mapTotal$concerts ), "concerts"] <- 0

head(cbind(map2016$NAME, map2016$concerts, map2017$concerts, map2018$concerts, mapTotal$concerts))
##      [,1]          [,2] [,3] [,4] [,5]
## [1,] "Alabama"     "0"  "0"  "0"  "0" 
## [2,] "Arizona"     "1"  "0"  "0"  "1" 
## [3,] "Arkansas"    "0"  "0"  "0"  "0" 
## [4,] "California"  "3"  "7"  "14" "24"
## [5,] "Colorado"    "0"  "0"  "0"  "0" 
## [6,] "Connecticut" "0"  "0"  "0"  "0"
#head(countTotal)
#leaflet map for 2016
pal2016 <- colorNumeric(palette = "YlOrRd", domain = map2016$concerts)
leaflet(map2016) %>%
  addTiles() %>%
  addPolygons(
    fillColor = ~pal2016(concerts),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    label =  ~NAME,
    popup = ~paste("State:", NAME, "<br>Concerts:", concerts),
    highlightOptions = highlightOptions(
      weight = 3,
      color = "#fc4e2a",
      fillColor = "#fc4e2a",
      fillOpacity = 1,
      bringToFront = TRUE
    )
  ) %>%
  addLegend(
    pal = pal2016,
    values = ~concerts,
    title = paste("Concerts 2016 <br> Total:", sum(map2016$concerts)),
    position = "bottomright"
  )
#leaflet map for 2017
pal2017 <- colorNumeric(palette = "Purples", domain = map2017$concerts)
leaflet(map2017) %>%
  addTiles() %>%
  addPolygons(
    fillColor = ~pal2017(concerts),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    label =  ~NAME,
    popup = ~paste("State:", NAME, "<br>Concerts:", concerts),
    highlightOptions = highlightOptions(
      weight = 3,
      color = "#54278f",
      fillColor = "#54278f",
      fillOpacity = 1,
      bringToFront = TRUE
    )
  ) %>%
  addLegend(
    pal = pal2017,
    values = ~concerts,
    title = paste("concerts 2017 <br> Total:", sum(map2017$concerts)),
    position = "bottomright"
  )
#leaflet map for 2018
pal2018 <- colorNumeric(palette = "Blues", domain = map2018$concerts)
leaflet(map2018) %>%
  addTiles() %>%
  addPolygons(
    fillColor = ~pal2018(concerts),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    label =  ~NAME,
    popup = ~paste("State:", NAME, "<br>Concerts:", concerts),
    highlightOptions = highlightOptions(
      weight = 3,
      color = "#08519c",
      fillColor = "#08519c",
      fillOpacity = 1,
      bringToFront = TRUE
    )
  ) %>%
  addLegend(
    pal = pal2018,
    values = ~concerts,
    title = paste("concerts 2018 <br> Total:", sum(map2018$concerts)),
    position = "bottomright"
  )
#leaflet cumulative map
palTotal <- colorNumeric(palette = "RdPu", domain = mapTotal$concerts)

mapTotal$popup_html <- paste0(
  "<b>", mapTotal$NAME, "</b><br>",
  "Total: ", mapTotal$concerts, "<br>",
  "2016: ", map2016$concerts,"<div style='background:#fd8d3c;height:10px;width:", map2016$concerts * 10, "px'></div><br>", 
  "2017: ",map2017$concerts , "<div style='background:#807dba;height:10px;width:", map2017$concerts * 10, "px'></div><br>", 
  "2018: ", map2018$concerts ,"<div style='background:#2171b5;height:10px;width:", map2018$concerts * 10, "px'></div>" 
)

leaflet(mapTotal) %>%
  addTiles() %>%
  addPolygons(
    fillColor   = ~palTotal(concerts),
    fillOpacity = 0.7,
    color       = "black",
    weight      = 1,
    label =  ~NAME,
    popup       = ~popup_html,
    highlightOptions = highlightOptions(
      weight = 3,
      color = "#7a0177",
      fillColor = "#7a0177",
      fillOpacity = .9,
      bringToFront = TRUE
    )
  ) %>%
  addLegend(pal = palTotal, 
            values = ~concerts, 
            title = paste("Total concerts <br>2016-2018:", sum(mapTotal$concerts)), 
            position = "bottomright"
            )
#static combined panel of all three maps, side by side
plot2016 <- ggplot(map2016)+
  geom_sf(aes(fill = concerts), color = "black", size= 0.5)+
  scale_fill_distiller(palette = "YlOrRd", name = "Concerts")+
  labs(title = "2016")+
  theme_void()

plot2017 <- ggplot(map2017)+
  geom_sf(aes(fill = concerts), color = "black", size= 0.5)+
  scale_fill_distiller(palette = "Purples", name = "Concerts")+
  labs(title = "2017")+
  theme_void()

plot2018 <- ggplot(map2018)+
  geom_sf(aes(fill = concerts), color = "black", size= 0.5)+
  scale_fill_distiller(palette = "Blues", name = "Concerts")+
  labs(title = "2018")+
  theme_void()

panel <- ggarrange (
  plot2016, plot2017, plot2018,
  ncol = 3,
  nrow = 1,
  legend = "bottom"
)
panel

#saves
ggsave("map2016finalized.png", plot = plot2016, width = 10, height = 5, dpi = 300)
ggsave("map2017finalized.png", plot = plot2017, width = 10, height = 5, dpi = 300)
ggsave("map2018finalized.png", plot = plot2018, width = 10, height = 5, dpi = 300)
ggsave("maptotalfinalized.png", plot = , width = 10, height = 5, dpi = 300)
ggsave("panel.png", plot = panel, width = 10, height = 5, dpi = 300)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.