Here is a solution to drawing and coloring territories based on zip code using ggplot2
library(dplyr)
library(ggplot2)
library(stringr)
library(rgdal)
library(maptools)
shapefile <- readOGR(dsn = "cb_2018_us_zcta510_500k", layer = "cb_2018_us_zcta510_500k")
# Create a dataframe of zip codes. Be sure to have leading zeros.
zip_A <- read.csv("Territory A - Zip Codes.csv") %>%
mutate(Zipcode = str_pad(Zipcode, 5, side = "left", "0"))
head(zip_A)
## Zipcode
## 1 35004
## 2 35005
## 3 35006
## 4 35007
## 5 35010
## 6 35011
# Used for drawing state borders later
states <- map_data("state")
head(states)
## long lat group order region subregion
## 1 -87.46201 30.38968 1 1 alabama <NA>
## 2 -87.48493 30.37249 1 2 alabama <NA>
## 3 -87.52503 30.37249 1 3 alabama <NA>
## 4 -87.53076 30.33239 1 4 alabama <NA>
## 5 -87.57087 30.32665 1 5 alabama <NA>
## 6 -87.58806 30.32665 1 6 alabama <NA>
# Filter the shape file to the zip codes to be drawn.
# Create a dataframe including latitude & longitude and the order they are drawn.
uspoly <- subset(shapefile, ZCTA5CE10 %in% zip_A$Zipcode)
uspoly$group <- substr(uspoly$ZCTA5CE10, 1,3)
uspoly$ZCTA5CE10 <- droplevels(uspoly$ZCTA5CE10)
uspoly.union <- unionSpatialPolygons(uspoly, uspoly$group)
uspolyfort_A <- fortify(uspoly.union)
usmap_A <- uspolyfort_A %>%
select(long, lat, order, group) %>%
mutate(Territory_A = "In") %>%
mutate(group = as.numeric(group)) %>%
bind_rows(states %>%
select(long, lat, order, group) %>%
mutate(Territory_A = "Out") %>%
mutate(group = group + max(as.numeric(uspolyfort_A$group))))
head(usmap_A)
## long lat order group Territory_A
## 1 -72.19267 42.09154 1 1 In
## 2 -72.19304 42.08799 2 1 In
## 3 -72.19189 42.08489 3 1 In
## 4 -72.19277 42.08192 4 1 In
## 5 -72.19375 42.08154 5 1 In
## 6 -72.19307 42.07942 6 1 In
Polygons are drawn by connecting latitude and longitude in order and colored by the group
ggplot(usmap_A) +
geom_polygon(aes(x = long, y = lat, group = group, fill = Territory_A, color = Territory_A)) +
scale_colour_manual(aesthetics = "fill", values=c("cornflowerblue", NA)) +
scale_colour_manual(aesthetics = "color", values=c("cornflowerblue", "black")) + #Draws borders. (Zip border is the same as fill color for a cleaner look)
coord_fixed(1.3) +
guides(fill=FALSE) +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(),
axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_text(hjust = .5), legend.position ="none") +
labs(title = "Territory A")
# Territory B Map
zip_B <- read.csv("Territory B - Zip Codes.csv") %>%
distinct() %>%
mutate(Zipcode = str_pad(Zipcode, 5, 00000, side = "left"))
uspoly <- subset(shapefile, ZCTA5CE10 %in% zip_B$Zipcode)
uspoly$group <- substr(uspoly$ZCTA5CE10, 1,3)
uspoly$ZCTA5CE10 <- droplevels(uspoly$ZCTA5CE10)
uspoly.union <- unionSpatialPolygons(uspoly, uspoly$group)
uspolyfort_B <- fortify(uspoly.union)
usmap_B <- uspolyfort_B %>%
select(long, lat, order, group) %>%
mutate(Territory_B = "In") %>%
mutate(group = as.numeric(group)) %>%
bind_rows(states %>%
select(long, lat, order, group) %>%
mutate(Territory_B = "Out") %>%
mutate(group = group + max(as.numeric(uspolyfort_B$group)))) %>%
filter(between(lat, 20, 50)) %>% #This will exclude zip codes outside the continental USA
filter(between(long, -130, -65))
ggplot(usmap_B) +
geom_polygon(aes(x = long, y = lat, group = group, fill = Territory_B, color = Territory_B)) +
scale_colour_manual(aesthetics = "fill", values=c("firebrick2", NA)) +
scale_colour_manual(aesthetics = "color", values=c("firebrick2", "black")) +
coord_fixed(1.3) +
guides(fill=FALSE) +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(),
axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_text(hjust = .5), legend.position ="none") +
labs(title = "Territory B")
Creating a map to show overlapping territories
# This process merges the overlapping zip codes and creates 2 groups for A & B
zip_both <- inner_join(zip_A, zip_B, by = "Zipcode")
uspoly <- subset(shapefile, ZCTA5CE10 %in% zip_both$Zipcode)
uspoly$group <- substr(uspoly$ZCTA5CE10, 1,3)
uspoly$ZCTA5CE10 <- droplevels(uspoly$ZCTA5CE10)
uspoly.union <- unionSpatialPolygons(uspoly, uspoly$group)
uspolyfort_both <- fortify(uspoly.union)
zip_A_anti <- anti_join(zip_A, zip_B, "Zipcode")
uspoly <- subset(shapefile, ZCTA5CE10 %in% zip_A_anti$Zipcode)
uspoly$group <- substr(uspoly$ZCTA5CE10, 1,3)
uspoly$ZCTA5CE10 <- droplevels(uspoly$ZCTA5CE10)
uspoly.union <- unionSpatialPolygons(uspoly, uspoly$group)
uspolyfort_A_anti <- fortify(uspoly.union)
zip_B_anti <- anti_join(zip_B, zip_A, "Zipcode")
uspoly <- subset(shapefile, ZCTA5CE10 %in% zip_B_anti$Zipcode)
uspoly$group <- substr(uspoly$ZCTA5CE10, 1,3)
uspoly$ZCTA5CE10 <- droplevels(uspoly$ZCTA5CE10)
uspoly.union <- unionSpatialPolygons(uspoly, uspoly$group)
uspolyfort_B_anti <- fortify(uspoly.union)
Ordering the groups is important. The state borders must be drawn last so they lay on top of the territory outlines.
usmap_both_1 <- uspolyfort_both %>%
select(long, lat, order, group) %>%
mutate(Territory = "Territory A & B") %>%
mutate(group = as.numeric(group)) %>%
bind_rows(uspolyfort_A_anti %>%
mutate(group = as.numeric(group)) %>%
mutate(Territory = "Territory A") %>%
mutate(group = group + max(as.numeric(uspolyfort_both$group))))
usmap_both_2 <- usmap_both_1 %>%
bind_rows(uspolyfort_B_anti %>%
mutate(group = as.numeric(group)) %>%
mutate(Territory = "Territory B") %>%
mutate(group = group + max(as.numeric(usmap_both_1$group))))
usmap_both_3 <- usmap_both_2 %>%
bind_rows(states %>%
select(long, lat, order, group) %>%
mutate(Territory = "No Coverage") %>%
mutate(group = group + max(as.numeric(usmap_both_2$group)))) %>%
filter(between(lat, 20, 50)) %>%
filter(between(long, -130, -65)) %>%
mutate(Territory = factor(Territory, levels = c("Territory A", "Territory B", "Territory A & B", "No Coverage"))) %>%
select(long, lat, order, group, Territory)
head(usmap_both_3)
## long lat order group Territory
## 1 -72.19267 42.09154 1 1 Territory A & B
## 2 -72.19304 42.08799 2 1 Territory A & B
## 3 -72.19189 42.08489 3 1 Territory A & B
## 4 -72.19277 42.08192 4 1 Territory A & B
## 5 -72.19375 42.08154 5 1 Territory A & B
## 6 -72.19307 42.07942 6 1 Territory A & B
tail(usmap_both_3)
## long lat order group Territory
## 1398387 -106.3295 41.00659 15594 4653 No Coverage
## 1398388 -106.8566 41.01232 15595 4653 No Coverage
## 1398389 -107.3093 41.01805 15596 4653 No Coverage
## 1398390 -107.9223 41.01805 15597 4653 No Coverage
## 1398391 -109.0568 40.98940 15598 4653 No Coverage
## 1398392 -109.0511 40.99513 15599 4653 No Coverage
ggplot(usmap_both_3) +
geom_polygon(aes(x = long, y = lat, group = group, fill = Territory, color = Territory)) +
scale_fill_manual(values=c(`Territory A` = "cornflowerblue", `Territory B` = "firebrick2", `Territory A & B` = "purple",
`No Coverage` = NA)) +
scale_color_manual(values=c(`Territory A` = "cornflowerblue", `Territory B` = "firebrick2", `Territory A & B` = "purple",
`No Coverage` = "black", `Businesses with No Coverage` = "black")) +
coord_fixed(1.3) +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(),
axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_text(hjust = .5)) +
labs(title = "Territory Overlap")
