R code for creating maps of M.ovi disease distribution from WHDB

Movi_map.r

library(ggplot2)
library(rgdal)
library(maptools)
library(dplyr)
library(scales)
library(ggmap)

nvranges <- readShapeSpatial("data/NV_Ranges/NV_MountainRanges.shp")
n.map <- fortify(nvranges, region = "RANGE")
state <- readShapeSpatial("data/NV_State/NV_Admin_State.shp")
mstate <- fortify(state)

# MAP OF NV MOUNTAIN RANGES. SIMPLE, BUT SCALE ISN'T CORRECT
ggplot() +
  geom_path(data = n.map, aes(x = long, y = lat, group = group), color = "black", size = .25)

# MAP OF NV MOUNTAIN RANGES. CORRECT SCALE. HOWEVER COORD_MAP() DOESN'T WORK FOR ME???
ggplot() +
  geom_path(data = n.map, aes(x = long, y = lat, group = group), color = "black", size = .25) +
  coord_equal()

# MAP OF NV RANGES, WITH STATE BORDER. 
ggplot() +
  geom_path(data = n.map, aes(x = long, y = lat, group = group), color = "darkred", size = .25) +
  geom_path(data = mstate, aes(x = long, y = lat, group = group), color = "black", size = .25) +
  coord_equal()

# READING IN MOVI ELISA DATA
elisa <- read.csv("data/bhs_elisa.csv", strip.white = T)

gbrange <- elisa %>%
  group_by(RANGE) %>%
  summarise(N = n(),
            Elisa_Pos = sum(result == "Detected"),
            Elisa_Prv = Elisa_Pos/N)

range <- elisa %>%
  group_by(RANGE, SampleYear) %>%
  summarise(N = n(),
            Elisa_Pos = sum(result == "Detected"),
            Elisa_Prv = Elisa_Pos/N)

range2014 <- subset(range, SampleYear == 2014)
range2013 <- subset(range, SampleYear == 2013)
range2012 <- subset(range, SampleYear == 2012)
range2011 <- subset(range, SampleYear == 2011)
range2010 <- subset(range, SampleYear == 2010)

# TESTING THE CONTINUOUS SCALE, IT ISN'T WORKING ON THE MAP
range2014['f_test'] <- c(20:42)

ggplot(range2014, aes(x = f_test, y = N, fill = f_test)) +
  geom_point() +
  scale_fill_continuous(low = "green", high = "blue")
  

# MAPING PREVALANCE DATA WITH GGPLOT

ggplot() +
  geom_map(data = range2014, aes_string(map_id = "RANGE", fill = range2014$Elisa_Prv), map = n.map) + 
  geom_path(data = n.map, aes(x = long, y = lat, group = group), color = "darkred", size = .25) +
  geom_path(data = mstate, aes(x = long, y = lat, group = group), color = "black", size = .25) +
  scale_fill_gradient(low = "white", high = "blue", limits = c(0.0, 1.0)) +
  coord_equal()

# USE GEOM POLYGON INSTEAD
# RIGHT JOIN INSTEAD OF LEFT BECAUSE I ONLY WANT THE MOUNTAIN RANGE ROWS THAT MATCH THE DATAFRAME
plotdata <- right_join(n.map, range2014, by = c("id" = "RANGE"))

ElisaMap <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plotdata, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2014 Sample Year", fill = "")
  
ElisaMap

ggsave(ElisaMap, file = "Elisa.png", width = 6, height = 6, type = "cairo-png")

# MAPS FOR EACH YEAR

range2014 <- subset(range, SampleYear == 2014)
range2013 <- subset(range, SampleYear == 2013)
range2012 <- subset(range, SampleYear == 2012)
range2011 <- subset(range, SampleYear == 2011)
range2010 <- subset(range, SampleYear == 2010)

plot2014 <- right_join(n.map, range2014, by = c("id" = "RANGE"))
plot2013 <- right_join(n.map, range2013, by = c("id" = "RANGE"))
plot2012 <- right_join(n.map, range2012, by = c("id" = "RANGE"))
plot2011 <- right_join(n.map, range2011, by = c("id" = "RANGE"))
plot2010 <- right_join(n.map, range2010, by = c("id" = "RANGE"))
all <- right_join(n.map, gbrange, by = c("id" = "RANGE"))

mapall <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = all, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence All Years", fill = "")
map2014 <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plot2014, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2014 Sample Year", fill = "")
map2013 <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plot2013, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2013 Sample Year", fill = "")
map2012 <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plot2012, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2012 Sample Year", fill = "")
map2011 <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plot2011, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2011 Sample Year", fill = "")
map2010 <- ggplot() +
  geom_polygon(data = mstate, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = .1) +
  #geom_polygon(data = n.map, aes(x = long, y = lat, group = group), fill = NA, color = "darkred", size = .1) +
  geom_polygon(data = plot2010, aes(x = long, y = lat, group = group, fill = Elisa_Prv), color = "darkred", size = .1) +
  coord_equal() +
  scale_fill_gradient2(low = muted("slategray2"), high = "slategray4") +
  theme_nothing(legend = T) +
  labs(title = "M.ovi Prevalence 2010 Sample Year", fill = "")


map2010
map2011
map2012
map2013
map2014
mapall

ggsave(map2010, file = "2010.png", width = 5, height = 4.5, type = "cairo-png")
ggsave(map2011, file = "2011.png", width = 5, height = 4.5, type = "cairo-png")
ggsave(map2012, file = "2012.png", width = 5, height = 4.5, type = "cairo-png")
ggsave(map2013, file = "2013.png", width = 5, height = 4.5, type = "cairo-png")
ggsave(map2014, file = "2014.png", width = 5, height = 4.5, type = "cairo-png")
ggsave(mapall, file = "all.png", width = 5, height = 4.5, type = "cairo-png")