mimsy.io

# Load up some libraries
library(dplyr, warn.conflicts = FALSE)
library(sp)
library(rgdal, warn.conflicts = FALSE)
library(maps)
library(rjson)
library(geojsonio, warn.conflicts = FALSE)
library(ggplot2)
library(broom)
library(stringr)
library(RColorBrewer)

sinclair_station_all_path <- "sinclair_station_all.csv"
sinclair_station_all_df <- read.csv(sinclair_station_all_path)

head(sinclair_station_all_df)

# Get the shapefile
dma_geo_path <- "Interactive_DMA_Map/nielsengeo.json"
dma_geo_sh <- geojson_read(dma_geo_path, what = "sp")

head(dma_geo_sh, n = 1) # Show only one Polygon in the DataFrame

# Convert the shapefile into a csv
dma_geo_sh@data$id <- rownames(dma_geo_sh@data)
dma_geo_df <- merge(tidy(dma_geo_sh), dma_geo_sh@data, by = "id")

head(dma_geo_df)

# Get some pieces of the Portland-Auburn DMA
piece1_df <- dma_geo_df[dma_geo_df$group == 500.1,]
piece2_df <- dma_geo_df[dma_geo_df$group == 500.2,]
piece34_df <- dma_geo_df[dma_geo_df$group == 500.34,]

# Get top three rows from each piece
example_portland_df <- rbind(head(piece1_df, n = 3), 
                             head(piece2_df, n = 3),
                             head(piece34_df, n = 3))
# Print
example_portland_df

# Set the size of output graphics
library(repr)
options(repr.plot.width=8, repr.plot.height=5)

# Define a function that creates a map using the Albers projection
ggalbers <- function() {
    ggplot() +
    coord_map("albers", lat0 = 37.5, lat1 = 29.5) +
    theme(plot.title = element_text(hjust = 0.1, color="#666666", face="bold", size=20), 
        axis.text = element_blank(), 
        axis.line = element_blank(), 
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.border = element_blank(),
        panel.grid = element_blank(),
        panel.background = element_blank()
    )
}

# Get only the Portland-Auburn DMA
dma_portland_df <- dma_geo_df[dma_geo_df$dma == 500,]

# Generate a graph of the Portland-Auburn DMA
ggalbers() +
    geom_polygon(data = dma_portland_df, 
                 color = "red",
                 fill = "gray30",
                 aes(x = long, y = lat, group = group)) +
    ggtitle("Portland-Auburn DMA")

# Generate a map for TV ownership by Neilson market
ggalbers() +
    geom_polygon(data = dma_geo_df, 
                 color = "gray80",
                 aes(x = long, y = lat, group = group, fill = tvperc)) +
    ggtitle("Neilson DMAs - Percent TV") +
    guides(fill = guide_legend(title = "Percentage\nTV Ownership"))

# Generate a map for cable subscriptions by Neilson market
ggalbers() +
    geom_polygon(data = dma_geo_df, 
                 color = "gray80",
                 aes(x = long, y = lat, group = group, fill = cableperc)) +
    ggtitle("Neilson DMAs - Percent Cable") +
    guides(fill = guide_legend(title = "Percentage\nCable Subs"))

# Show markets in `sinclair_station_all_df` for Illinois
IL_mask <- str_sub(sinclair_station_all_df$Market, start = -2) == "IL"
unique(sinclair_station_all_df[IL_mask, "Market"])

# Show markets in `dma_geo_df` for Illinois
IL_mask <- sapply(dma_geo_df$dma1, function(x) grepl("IL", x))
unique(dma_geo_df[IL_mask, "dma1"])

# Get 5 rows from `sinclair_station_all_df`
test_df <- na.omit(sinclair_station_all_df)
test_samples <- sample(1:nrow(test_df), 5)

test_df[test_samples,]

# Convert the data frame to a `SpatialPoints` object
test_sh <- SpatialPoints(test_df[test_samples, c("Longitude", "Latitude")])

# Make sure that the projection from `test_sh` is the same as in `dma_geo_sh`
proj4string(test_sh) <- proj4string(dma_geo_sh)

test_sh

results_df <- test_sh %over% dma_geo_sh
results_df

# A function that looks through `geo_sh` (list of polygons) for the first polygon that
# contains the coordinates in `lon` and `lat`. Once found, it extracts the wanted attribute.
attribute_of_locations <- function(geo_sh, lon, lat, attribute) {
    
    # Turn `lon`, `lat` into a `SpatialPoints` object
    points_sh <- SpatialPoints(cbind(lon, lat))
    proj4string(points_sh) <- proj4string(geo_sh)
    
    # Overlay `points_sh` on top of `geo_sh`
    results_df <- points_sh %over% geo_sh
    
    results_df[,attribute]
}

sinclair_station_dma_df <- sinclair_station_all_df %>% 
    select(c("Market", "Station", "Latitude", "Longitude", "Year.Start", "Year.End")) %>%
    filter(!is.na(Latitude) & !is.na(Longitude)) %>%
    mutate(DMA = attribute_of_locations(dma_geo_sh, Longitude, Latitude, "dma1")) %>%
    mutate(DMA.Code = attribute_of_locations(dma_geo_sh, Longitude, Latitude, "dma"))

head(sinclair_station_dma_df)

sinclair_dma_df <- sinclair_station_dma_df %>% 
    group_by(DMA, DMA.Code) %>%
    summarize(Year.Start = min(Year.Start), Year.End = max(Year.End))

head(sinclair_dma_df)

dma_attributes_path <- "Interactive_DMA_Map/tv_2017.csv"
dma_attributes_df <- read.csv(dma_attributes_path)
colnames(dma_attributes_df) <- c("DMA.Code", "Rank", "DMA", "TV.Homes", "Percent.of.US")

head(dma_attributes_df)

sinclair_dma_attr_df <- merge(dma_attributes_df[c("DMA.Code", "TV.Homes", "Percent.of.US")], 
                              sinclair_dma_df[c("DMA.Code", "Year.Start", "Year.End")], 
                              by = "DMA.Code")
sinclair_dma_attr_df$TV.Homes <- gsub(",", "", sinclair_dma_attr_df$TV.Homes)
sinclair_dma_attr_df$TV.Homes <- as.numeric(as.character(sinclair_dma_attr_df$TV.Homes))
head(sinclair_dma_attr_df)

sapply(sinclair_dma_attr_df, class)

# Get temporal limits of the data
max_year <- max(sinclair_dma_attr_df$Year.Start, na.rm = TRUE)
min_year <- min(sinclair_dma_attr_df$Year.Start, na.rm = TRUE)

# Get rid of data that has missing year values
no_na_year_df <- sinclair_dma_attr_df %>% filter(!is.na(Year.Start))

# Create aliases for the two year columns
start_year <- no_na_year_df$Year.Start
end_year <- no_na_year_df$Year.End

# Make a new data.frame and fill it using a for-loop
# This is a slow way to do it but our dataset is not huge 
sinclair_population_df <- data.frame()
for(year in min_year:max_year) {
    # Generate the station count
    # Remember to take care of possible NA values in `end_year`
    new_row <- sinclair_dma_attr_df %>%
        filter( Year.Start <= year & (Year.End >= year | is.na(Year.End)) ) %>%
        summarize(TV.Homes.Tot = sum(TV.Homes), Percent.of.US.Tot = sum(Percent.of.US))
    new_row$Year = c(year)
    sinclair_population_df <- rbind(sinclair_population_df, new_row)
}

sinclair_population_df[sample(1:nrow(sinclair_population_df), 5),]

# Create the plot
sinclair_population_plot <- ggplot(data = sinclair_population_df) +
    geom_line(aes(x = Year, y = Percent.of.US.Tot), color = brewer.pal(3, "Set2")[2]) +
    ggtitle("Sinclair Viewership") + 
    xlab("Year") + ylab("Percent of USA Population") +
    theme(plot.title = element_text(color="#666666", face="bold", size=20))

sinclair_population_plot

tail(sinclair_population_df$Percent.of.US.Tot, n = 5)

max(sinclair_population_df$Percent.of.US.Tot)

max(sinclair_population_df$TV.Homes.Tot)

head(dma_geo_df)

sinclair_final_df <- merge(sinclair_dma_df, dma_geo_df, by.x = "DMA", by.y = "dma1")
sinclair_final_df <- sinclair_final_df %>%
    select(DMA, id, long, lat, order, hole, piece, group, Year.Start, Year.End)
sinclair_final_df <- sinclair_final_df[order(sinclair_final_df$id, sinclair_final_df$order),]

head(sinclair_final_df)

plot_sinclair_dmas <- function(year) {
    # Separate the spatial data into four datasets
    plot_df <- sinclair_final_df %>%
        mutate(Acquired = case_when(
            (Year.Start < year) & (Year.End > year | is.na(Year.End)) ~ "acquired\nprevious years",
            Year.Start == year ~ "acquired\nthis year",
            Year.End == year ~ "sold\nthis year",
            TRUE ~ "not held")
        )
    
    # We don't want to show the DMAs that were not held during the given year
    held_mask <- plot_df$Acquired != "not held"
    plot_df <- plot_df[held_mask, ]
    factor_levels <- c("acquired\nprevious years", "acquired\nthis year", "sold\nthis year")
    plot_df$Acquired <- factor(plot_df$Acquired, levels = factor_levels)
    
    # Pick some nice colors
    nice_colors <- brewer.pal(3, "Set2")
    
    # Make the plot
    ggplot() +
        coord_map("albers", lat0 = 37.5, lat1 = 29.5) + 
        geom_polygon(data = dma_geo_df,
                     aes(x = long, y = lat, group = group),
                     color = "gray60") +
        geom_polygon(data = plot_df,
                     aes(x = long, y = lat, group = group, fill = Acquired),
                     color = "gray60") +
        scale_fill_manual(values = nice_colors, drop = FALSE) +
        theme(plot.title = element_text(hjust = 0.1, color="#666666", face="bold", size=20), 
            axis.text = element_blank(), 
            axis.line = element_blank(), 
            axis.ticks = element_blank(),
            axis.title = element_blank(),
            panel.border = element_blank(),
            legend.position = c(0.93, 0.5),
            legend.key.size = unit(1, "cm"),
            legend.text = element_text(size=12),
            legend.title = element_text(size=12),
            panel.grid = element_blank(),
            panel.background = element_blank()) +
        ggtitle(paste("Sinclair DMAs in", year))
}

plot_sinclair_dmas(2013)

IMAGE_OUTPUT_PATH <- "Images"

ggsave(filename = 'Sinclair_Viewership.png', 
    plot = sinclair_population_plot,
    path = IMAGE_OUTPUT_PATH,
    width = unit(8, "cm"),
    height = unit(5, "cm"))

max_year <- max(sinclair_final_df$Year.Start, na.rm = TRUE)
min_year <- min(sinclair_final_df$Year.Start, na.rm = TRUE)

for(year in min_year:max_year) {
    image <- plot_sinclair_dmas(year)
    image.name <- paste("Sinclair_DMAs", year, sep = "_")
    image.name <- paste(image.name, ".png", sep = "")
    ggsave(filename = image.name, 
        plot = image,
        path = IMAGE_OUTPUT_PATH,
        width = unit(14, "cm"),
        height = unit(10, "cm"))
}