gistfile1.r

#source MASTER
setwd("~/your/working/directory") #insert your working directory

#load libraries
library(rgdal)
library(maptools)
library(maps)
library(ggplot2)
library(plyr)
library(grid)
library(sp)
library(deldir)
library(Hmisc) #for capitalize function

#choose type of healthcare site to make voronoi polygons
typeOfSite <- "hospital"

##Load population data by ward and turn into dotmap
SAWards <- readShapePoly("./Ward_Pop2/Ward_Pop2.shp")   #Polygon shapefile with field for population called Cen_2011
dots.pop <- dotsInPolys(SAWards, SAWards$Cen_2011/100)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
pop.df <- data.frame(coordinates(dots.pop)[,1:2]) # extract the dataframe for ggplot

#FUNCTION to create voronoi polygons
#The following function takes a SpatialPointsDataFrame as input,
#and returns a SpatialPolygonsDataFrame that represents the Voronoi tessellation
#of the input point layer.
voronoipolygons = function(layer) {
  crds = layer@coords
  z = deldir(crds[,1], crds[,2], plotit=TRUE)
  w = tile.list(z)
  polys = vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP = SpatialPolygons(polys)
  voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                         y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
                                                                                      function(x) slot(x, 'ID'))))
}

#load clinic data
sites <- read.csv("./mesedi.csv")  #this data will be used to make the voronoi polygons.
#This script uses the fields latitude and longitude for the analysis later

#Filter and clean up clinic data
sites <- subset(sites, sites$typeOfSite ==typeOfSite) #filter by user defined variable typeOfSite (see line 5)
sites <- subset(sites, Country=="South Africa") #filter for only South Africa
sites <- sites[!duplicated(sites$latitude),] #remove sites with duplicate lat/long, which would break our voronoi algorithm

#add dummy sites
#it can be useful to add clincis outside the bounds of your shape, like for example South Africa
  #to prevent extra large polygons on the periphery of your voronoi diagram.
  #I've commented them out here, but this is one simple way of doing it.
#sites[dim(sites)[1]+1,c("latitude","longitude")]=c(-24.15,28.8) #add dummy point further north
#sites[dim(sites)[1]+1,c("latitude","longitude")]=c(-25.5,30) #add dummy point further east

#create data.frame 'coords' and run voronoi script
coords <- as.data.frame(cbind(sites$longitude, sites$latitude)) #coordinates of each clinic
colnames(coords) = c("lat", "long") #update column names for coords
points <- SpatialPoints(coords) #turn this into type SpatialPoints
spdf <- SpatialPointsDataFrame(points,sites) #Make a data.frame that has the points (as SpatialPoints) and other fields about the clinic
voronoiPolys <- voronoipolygons(spdf) # run our voronoipolygons function on this data.frame

#remove dummy points
#if you added dummy point, now is a good time to remove them so they don't get added to the map.
#voronoiPolys <- voronoiPolys[1:110,] #removes the two dummy points added in the section above.

#map sites
clinicMap <- ggplot() +
  geom_point(data=coords, size=2.5, colour="red", aes(x = lat, y = long)) + #define style for the sites
  labs(title = "Locations of Sites") + #optional, set the title text
  theme(plot.title = element_text(size = rel(2))) + #optional, set the size of the title
  coord_cartesian(xlim = c(14.85297, 34.38551), ylim = c(-35.89747, -21.10941)) #optional, setting the bounding box to be mapped

  #save sites map as png
  ggsave(plot=clinicMap,filename="./clinicMap.png", width=12,height=9)

#map sites and Voronoi
sitesAndVoronoi <- ggplot(voronoiPolys, aes(x = long, y = lat)) +
  geom_polygon(aes(group = group), colour = I("grey65"),
               size=0.4, fill = "white", alpha = 0.3) + coord_equal() + #add voronoi polygons
  labs(title = paste("Voronoi Diagram of "),capitalize(typeOfSite), "s", sep="") + #pull in the type of site into the title
  theme(plot.title = element_text(size = rel(2))) + # set title font size
  geom_point(data=coords, size=1.5, colour="red", aes(x=lat,y=long)) + #add sites
  expand_limits(x = NULL, y = NULL) + #remove buffer on the edges of the map
  coord_cartesian(xlim = c(14.85297, 34.38551), ylim = c(-35.89747, -21.10941)) #optional, set bounds for map

  #save sitesAndVoronoi map as png
  ggsave(plot=sitesAndVoronoi,filename="./sitesAndVoronoi.png", width=12,height=9)

#map Voronoi diagram and population dots
voronoiAndDots <- sitesAndVoronoi + geom_point(data=pop.df, aes(x=x,y=y), size=0.5, colour="blue", alpha=0.15) +
  geom_point(data=coords, size=1.5, colour="red") +
  labs(title = "Voronoi Diagram and Population Density (in dots)") +
  theme(plot.title = element_text(size = rel(2)))

  #Save this map as png
  ggsave(plot=voronoiAndDots,filename="./voronoiAndDots.png", width=12,height=9)

##The following steps are for counting the number of dots in each thessian polgon in the voronoi diagram

##Perform Count using 'over' and then bind this data back to the polygons
popCount <- sapply(over(voronoiPolys, SpatialPoints(pop.df), returnList = TRUE), length)
voronoiPolys <- spCbind(voronoiPolys, as.data.frame(popCount))

#you have to tweak the SpatialPolygonsDataFrame into a data.frame to map it using ggplot
# more on this here: https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles
voronoiPolys@data$id <- rownames(voronoiPolys@data)
voronoiPolys.points <- fortify(voronoiPolys, region="id")
voronoiPolys.df <- join(voronoiPolys.points, voronoiPolys@data, by="id")

#do the same thing with a shapefile to display the country boundary of SA in the final plot.
SABorder <- readShapePoly("./ZAF_adm/ZAF_adm0.shp")
SABorder@data$id <- rownames(SABorder@data)
SABorder.points <- fortify(SABorder, region="id")
SABorder.df <- join(SABorder.points, SABorder@data, by="id")

#remove an island off the southern coast of SA, to improve the final map.
SABorder.df <- subset(SABorder.df, lat>(-40))

#rename popCount column, so the legend is nicer
names(voronoiPolys.df)[names(voronoiPolys.df)=="popCount"]  <- "Population"

#create voronoiWithPopMap
voronoiWithPopMap <- ggplot(voronoiPolys.df, aes(x=long, y=lat)) + coord_equal() +
  geom_polygon(colour="black", size=0.1, aes(group=group, fill=Population)) + #colored voronoi
  geom_point(data=coords, size=0.7, colour="red", alpha=.7, aes(x=lat,y=long)) + #sites
  geom_path(data=SABorder.df, colour="white", size=0.3, alpha=0.5, aes(group=group)) + #SA Border
  labs(title = paste("Estimated Population Served per",capitalize(typeOfSite))) +
  theme(plot.title = element_text(size = rel(2)),
        legend.justification=c(1,0), legend.position=c(0.98,0.02), legend.key.size=unit(.05,"npc")) +
  expand_limits(x = NULL, y = NULL) +
  coord_cartesian(xlim = c(14.85297, 34.38551), ylim = c(-35.89747, -21.10941))

  #show voronoiWithPopMap map in plot panel
  voronoiWithPopMap

  #save voronoiWithPopMap map
  ggsave(plot=voronoiWithPopMap,filename=paste("./voronoiWithPopMap-",typeOfSite,".png", sep=""), width=12,height=9)