NICAR 2019: Mapping with R - Sinclair Case Study

Andrew Ba Tran

3/8/2019

We’re going to recreate the graphic from The Washington Post [story about Sinclair Broadcasting]((https://www.washingtonpost.com/graphics/2018/lifestyle/sinclair-broadcasting/):

Sinclair Broadcasting

Graphics editor Chris Alcantara used a mix of R, Excel, QGIS, and Illustrator to create this map.

We’re going to try to stick to R exclusively with the help of some packages like

  • jsonlite
  • tigris
  • dplyr
  • ggplot2
  • sf
  • ggrepel
  • shadowtext

Finding the data

There’s a map on Sinclair’s website showing off all their stations.

Grab the json URL that’s running the map.

Scraping the data

We’ll use the jsonlite package

library(tidyverse)
library(stringr)
library(jsonlite)

json_url <-"http://sbgi.net/resources/assets/sbgi/MetaverseStationData.json"

# json_url <- "data/MetaverseStationData.json"

stations <- fromJSON(json_url)

Turning the JSON into a clean data frame

Turning the JSON into a clean data frame

Let’s focus on Primary stations and extract the latitude and longitude data

primary_stations <- stations %>% 
  filter(Channel=="Primary") %>% 
  mutate(
    Location=str_replace(Location, "Point \\(", ""),
    lon=str_replace(Location, " .*", ""),
    lat=str_replace(Location, ".* ", ""),
    lat=str_replace(lat, "\\)", ""))

What have we got

glimpse(primary_stations)
## Observations: 192
## Variables: 28
## $ Call_Letter            <chr> "KAAS", "KAAS-LD", "KABB", "KAEF", "KAME"…
## $ Logo_List              <chr> "", "", "kabb_fox.jpg", "/resources/asset…
## $ Logo_Map               <chr> "", "", "kabb_fox_map.jpg", "/resources/a…
## $ Web_1st_URL            <chr> "http://www.foxkansas.com", "http://www.f…
## $ Web_Address            <chr> "http://www.foxkansas.com", "http://www.f…
## $ Station                <chr> "KAAS", "KAAS-LD", "KABB", "KAEF", "KAME"…
## $ Channel                <chr> "Primary", "Primary", "Primary", "Primary…
## $ Affiliation            <chr> "FOX", "FOX", "FOX", "ABC", "Nevada Sport…
## $ DMA                    <chr> "Wichita - Hutchinson, KS", "Wichita - Hu…
## $ DMA_Code               <chr> "678", "678", "641", "802", "811", "820",…
## $ DMA_Short              <chr> "Wichita_KS", "Wichita_KS", "San Antonio_…
## $ DMA_Rank               <int> 76, 76, 31, 195, 109, 22, 57, 122, 122, 1…
## $ Station_Status         <chr> "O&O", "O&O", "O&O", "O&O", "JSA", "O&O",…
## $ Station_Address        <chr> "316 North West Street, Wichita, KS 67203…
## $ Station_City           <chr> "Wichita", "Wichita", "San Antonio", "Eur…
## $ Station_State          <chr> "KS", "KS", "TX", "CA", "NV", "OR", "AR",…
## $ Station_Zip            <int> 67203, 67203, 78229, 95501, 89502, 97232,…
## $ Station_Logo           <chr> "sbg_noimage", "sbg_noimage", "kabb_fox",…
## $ Station_URL            <chr> "http://www.foxkansas.com, http://www.fox…
## $ Station_Phone_Number   <chr> "316-942-2424", "316-942-2424", "210-366-…
## $ Station_Fax_Number     <chr> "316-942-8927", "316-942-8927", "210-377-…
## $ Actual_RF_Channel      <chr> "17", "31", "30", "22", "20", "43", "22",…
## $ News_Schedule_Weekday  <chr> "", "", "0500-0900, 1100-1200, 2100-2200"…
## $ News_Schedule_Saturday <chr> "", "", "2100-2200", "", "", "0700-0900,1…
## $ News_Schedule_Sunday   <chr> "", "", "2100-2200", "", "", "0700-0900,0…
## $ Location               <chr> "-97.388134 37.68888)", "-97.388134 37.68…
## $ lon                    <chr> "-97.388134", "-97.388134", "-98.569649",…
## $ lat                    <chr> "37.68888", "37.68888", "29.490196", "40.…

What can we do with this data?

There are two data points that we can extract from this data set:

  • latitude and longitude for each station
  • the DMA, or Designated Market Area, each station is assigned to

We can map out the station locations quickly using the sf package to visualize latitude and longitude of each one on top of a map of the country downloaded using the tigris package.

Single out the stations

station_latlon <- 
  select(primary_stations, DMA_Short, Location, lon, lat) %>%
  unique() %>% 
  mutate(lon=as.numeric(lon)) %>% 
  mutate(lat=as.numeric(lat)) %>% 
  filter(!is.na(lon))

glimpse(station_latlon)
## Observations: 121
## Variables: 4
## $ DMA_Short <chr> "Wichita_KS", "San Antonio_TX", "Eureka_CA", "Reno_NV"…
## $ Location  <chr> "-97.388134 37.68888)", "-98.569649 29.490196)", "-124…
## $ lon       <dbl> -97.38813, -98.56965, -124.16704, -119.78651, -122.643…
## $ lat       <dbl> 37.68888, 29.49020, 40.80141, 39.51247, 45.52701, 34.7…

Download map of state boundaries

library(sf)
library(tigris)

options(tigris_class = "sf")

states <- states(cb=T)

# states <- readRDS("backup_data/stats.rds")


glimpse(states)

Refining the data

Filter out some territories and states from the data and set the Projection

states <- filter(states, !STUSPS %in% c("AK", "AS", "MP", "PR", "VI", "HI", "GU"))

# Converting the projection to Albers
states <- st_transform(states, 5070)

# Changing the projection of station_latlon so it matches the states sf dataframe map

station_latlon_projected <- station_latlon %>% 
  filter(!is.na(lon)) %>% 
  st_as_sf(coords=c("lon", "lat"), crs = "+proj=longlat") %>% 
  st_transform(crs=st_crs(states)) %>% 
  st_coordinates(geometry)

station_latlon <- cbind(station_latlon, station_latlon_projected)

Map it

ggplot(states) +
  geom_sf() +
  geom_point(data=station_latlon, aes(x=X, y=Y), color="blue") +
  theme_void() +
  theme(panel.grid.major=element_line(colour="transparent")) +
  labs(title="Station locations")

We’ve added a new map layer on top of the states boundaries map with the function: geom_point()

Map it

Next, figure out the scope

# Fixing a bad data point
primary_stations$DMA_Code <- ifelse(primary_stations$DMA_Short=="Lincoln_NE", 722, primary_stations$DMA_Code)

## Summarizing by DMA
dma_totals <- primary_stations %>% 
  group_by(DMA, DMA_Code) %>% 
  count() %>% 
  arrange(desc(n)) %>% 
  ungroup() %>% 
  rename(dma_code=DMA_Code) %>% 
  mutate(dma_code=as.numeric(dma_code))

head(dma_totals)
## # A tibble: 6 x 3
##   DMA                                       dma_code     n
##   <chr>                                        <dbl> <int>
## 1 Eugene, OR                                     801     6
## 2 Lincoln - Hastings - Kearney, NE               722     6
## 3 Wichita - Hutchinson, KS                       678     6
## 4 Champaign - Springfield - Decatur, IL          648     5
## 5 Chico-Redding                                  868     5
## 6 Yakima - Pasco - Richland - Kennewick, WA      810     5

What’s the DMA footprint for each station?

Check around on the internet and you’ll find sources for the DMA shapefile.

Here’s one from Harvard’s Dataverse.

I’ve already downloaded it to our project folder.

We’ll read it in using the st_read() function from the sf package.

Mapping DMAs

geo <- st_read("data/dma_2008/DMAs.shp")

# It doesn't have a CRS so we'll assign it one
st_crs(geo) <- 4326

# Converting the projection so it's Albers
geo <- st_transform(geo, 5070)

ggplot(geo) +
  geom_sf(color="red") +
  coord_sf()

Mapping DMAs

geo <- st_read("data/dma_2008/DMAs.shp")

# It doesn't have a CRS so we'll assign it one
st_crs(geo) <- 4326

# Converting the projection so it's Albers
geo <- st_transform(geo, 5070)

ggplot(geo) +
  geom_sf(color="red") +
  coord_sf()

Mapping Sinclair DMAs

Join the DMA station count to it so we can visualize it.

# Prepping a column name to join on
geo <- geo %>% 
  mutate(dma_code=as.numeric(as.character(DMA)))

geo <- left_join(geo, dma_totals, by="dma_code") %>% 
  filter(!is.na(n)) 

ggplot() +
  geom_sf(data=states, color="red", fill=NA) +
  geom_sf(data=geo, aes(fill=n)) +
  coord_sf()

Mapping Sinclair DMAs

Data prep

# Filtering out locations based on map
cities <- c("Portland", "Seattle", "Butte", "Boise", "Reno", "Fresno", "Bakersfield", 
            "Reno", "Salt Lake City", "Las Vegas", "El Paso", "Austin", "San Antonio",
            "Corpus Christi", "Oklahoma City", "Wichita", "Lincoln", "Sioux City", 
            "Minneapolis", "Green Bay", "Milwaukee", "Des Moines", "Springfield", 
            "St. Louis", "Little Rock", "Flint", "Columbus", "Lexington", "Birmingham",
            "Mobile", "Macon", "Asheville", "Charleston", "Buffalo", "Johnstown", "Baltimore",
            "Norfolk", "Conway", "Savannah", "Gainesville", "West Palm Beach", "Albany",
            "Providence", "Portland")

station_latlon_filtered <- station_latlon %>% 
  mutate(DMA_Short= gsub("_.*", "", DMA_Short)) %>% 
  mutate(DMA_Short= gsub("Bozeman", "Butte", DMA_Short)) %>% 
  mutate(DMA_Short= gsub("Champaign", "Springfield", DMA_Short)) %>% 
  mutate(DMA_Short= gsub("Myrtle Beach", "Conway", DMA_Short)) %>% 
  mutate(DMA_Short= gsub("West Palm", "West Palm Beach", DMA_Short)) %>% 
  filter(DMA_Short %in% cities) %>% 
  group_by(DMA_Short) %>% 
  filter(row_number()==1)

Data prep

glimpse(station_latlon_filtered)
## Observations: 41
## Variables: 6
## Groups: DMA_Short [41]
## $ DMA_Short <chr> "Wichita", "San Antonio", "Reno", "Portland", "Little …
## $ Location  <chr> "-97.388134 37.68888)", "-98.569649 29.490196)", "-119…
## $ lon       <dbl> -97.38813, -98.56965, -119.78651, -122.64375, -92.2710…
## $ lat       <dbl> 37.68888, 29.49020, 39.51247, 45.52701, 34.74508, 35.3…
## $ X         <dbl> -121246.74, -249168.34, -2005457.33, -2054419.36, 3385…
## $ Y         <dbl> 1628840.1, 716202.8, 2084644.2, 2794219.9, 1304925.7, …

Style prep

geo <- geo %>% 
  mutate(bin=case_when(
    n == 1 ~ "1",
    n == 2 ~ "2",
    n == 3 ~ "3",
    n == 4 ~ "4",
    n >= 5 ~ "5+"
  ))

library(ggrepel)
library(shadowtext)

Final map

We’re layering our different shape files:

  • The state borders
  • The DMA borders we combined with the counts of them from the summarized Sinclair JSON file
  • The locations of a few dozen stations we also extracted from the Sinclair JSON file

And we’re adding a bunch of functions in there for styling, like scale_fill_brewer() and geom_shadowtext() and geom_text_repel() and various theme options.

Final map

ggplot() +
  geom_sf(data=states, color="gray", fill=NA, size=.3) +
  geom_sf(data=geo, aes(fill=bin), color="light gray", size=.2) +
  scale_fill_brewer(palette = "Oranges", name="Number of Sinclair-owned TV stations") +
  geom_point(data=station_latlon_filtered, aes(x=X, y=Y), 
             color="dark gray", fill="white", shape=21) +
  geom_shadowtext(data=station_latlon_filtered, aes(x=X, y=Y, label=DMA_Short), 
             color="black", bg.color="white", vjust=-1, size=2.5) +
  geom_text_repel() +
  coord_sf() +
  theme_void() +
  theme(panel.grid.major = element_line(colour = 'transparent'),
        legend.position="top", legend.direction="horizontal") +
  labs(title="Sinclair's nationwide reach",
       subtitle="This map only highlights stations owned outright by Sinclair. 
It does not include the many stations licensed to other operators but managed by Sinclair.",
       caption="Sources: U.S. Securities and Exchange Commission, Nielsen, Television Bureau of Advertising, GeoCommons")

ggsave("sinclair_ggplot.png", width=10, height=6, units="in")

Final map

How does it compare?