The code below creates an animated GIF file showing a map of incident and prevalent ASF communes by week and an epidemic curve of the entire outbreak. The vertical bar on the epidemic curve shifts to the right according to week. Incident and prevalent communes for each week are shown on the map.
library(tidyr); library(dplyr); library(gridExtra); library(rgdal); library(ggplot2); library(animation); library(scales); library(plotly)
setwd("D:\\TEMP")
Load ASF case details, Vietnamese commune details, a map of Vietnamese provinces as an ESRI shapefile and specify the RAHO to analyse. If you want to plot data for the entire country, set sraho <- "all":
vnasf.df <- read.csv("VN_ASF_cases.csv", header = TRUE, sep = ",")
vncom.df <- read.table("VN_communes_v03.csv", header = TRUE, sep = ",")
vnadm1.ll <- readOGR(dsn = getwd(), layer = "VN_province_RAHO-LL")
sraho <- "RAHO1"
Function to format the horizontal and vertical axes labels nicely:
mformat <- function(){
function(x) format(x / 1000, digits = 2)
}
Re-project the digital map to UTM 48N:
vnadm1.utm <- spTransform(vnadm1.ll, CRS("+init=epsg:32648"))
Update the ASF case detail data frame with commune easting and northing coordinates:
vnasf.df$xcoord <- vncom.df$xcoord[match(vnasf.df$ccode, vncom.df$ccode)]
vnasf.df$ycoord <- vncom.df$ycoord[match(vnasf.df$ccode, vncom.df$ccode)]
vnasf.df$raho <- vncom.df$raho[match(vnasf.df$ccode, vncom.df$ccode)]
Format the dates. If dates misbehave check regional settings in Control Panel and change to Australian English:
vnasf.df$odate <- strptime(vnasf.df$odate, format="%d-%b-%y"); vnasf.df$odate <- as.Date(vnasf.df$odate, format = "%d-%b-%y")
vnasf.df$rdate <- strptime(vnasf.df$rdate, format="%d-%b-%y"); vnasf.df$rdate <- as.Date(vnasf.df$rdate, format = "%d-%b-%y")
vnasf.df$cdate <- strptime(vnasf.df$cdate, format="%d-%b-%y"); vnasf.df$cdate <- as.Date(vnasf.df$cdate, format = "%d-%b-%y")
Group the data by commune to return the date of onset of clinical signs of the first ASF case in each commune:
vnasft.df <- vnasf.df %>%
group_by(ccode)
vnasft.df <- data.frame(summarise(vnasft.df, raho = first(raho), edate = min(odate), xcoord = min(xcoord), ycoord = min(ycoord)))
vnasft.df <- vnasft.df[3:nrow(vnasft.df),]
Drop rows with missing values:
id <- complete.cases(vnasft.df)
vnasft.df <- vnasft.df[id,]
head(vnasft.df)
Select data for the RAHO of interest:
if(sraho != "all"){
id <- vnasft.df$raho == sraho
vnasft.df <- vnasft.df[id,]
id <- vnadm1.utm$RAHO == sraho
vnadm1.utm <- vnadm1.utm[id,]
}
vnadm1.df <- fortify(vnadm1.utm)
Set up a sequence of dates for plotting:
sdate <- seq(from = as.Date("1-Jan-2019", format = "%d-%b-%Y"), to = as.Date("30-Jun-2019", format = "%d-%b-%Y"), by = 7)
Loop through each of the dates listed in sdate and make a map of incident ASF communes (red) and prevalent ASF communes (grey):
saveGIF({for(i in 2:length(sdate)){
# Prevalent ASF communes:
id <- vnasft.df$edate < sdate[i - 1]
pvnasft.df <- vnasft.df[id,]
# Incident ASF communes (in the last week):
id <- vnasft.df$edate < sdate[i] & vnasft.df$edate >= sdate[i - 1]
ivnasft.df <- vnasft.df[id,]
# Map of incident and prevalent communes:
plot1 <- ggplot() +
geom_point(data = pvnasft.df, aes(x = xcoord, y = ycoord), col = "grey", pch = 1) +
geom_point(data = ivnasft.df, aes(x = xcoord, y = ycoord), col = "red", pch = 16) +
geom_polygon(data = vnadm1.df, aes(x = long, y = lat, group = group), col = "black", fill = "transparent") +
scale_x_continuous(name = "Easting (km)", labels = mformat()) +
scale_y_continuous(name = "Northing (km)", labels = mformat()) +
coord_fixed() +
theme_bw()
# Epidemic curve:
plot2 <- ggplot(data = vnasft.df, aes(x = edate)) +
geom_histogram(color = "grey", fill = "dark blue", binwidth = 2) +
scale_x_date(breaks = date_breaks("2 months"), labels = date_format("%d %b"), name = "Date of onset of signs") +
scale_y_continuous(breaks = seq(from = 0, to = 150, by = 10), name = "Number of communes") +
geom_vline(aes(xintercept = sdate[i]), linetype = "dashed", col = "red", lwd = 1) +
theme_bw()
# Arrange the two plots side by side:
print(grid.arrange(plot1, plot2, ncol = 2, widths = c(5,5))
)}
}, movie.name = "VN_ASF_point_map_epidemic_curve.gif", interval = 1)