Step-by-step plotting choropleth map of Nepal

map.choropleth.nepal-updated.R

library(rgdal)
library(ggplot2)
library(dplyr)

# clone NepalMaps from https://github.com/anjesh/NepalMaps 

# read shapefile
nepal_shp <- readOGR(dsn="NepalMaps/baselayers/NPL_adm", layer="NPL_adm3", stringsAsFactors = FALSE)
# fortify shapefile data to data frame
shp_df <- fortify(nepal_shp, region = "NAME_3")

# update the HPI data in districts.csv as per http://data.opennepal.net/content/human-poverty-index-value-districts-2011
hpi.csvdata <- read.csv("https://github.com/opennepal/odp-poverty/raw/master/Human%20Poverty%20Index%20Value%20by%20Districts%20(2011)/data.csv", stringsAsFactors = FALSE)
hpi.data <- hpi.csvdata %>% 
  filter(Sub.Group == "HPI") %>% 
  select(id = District, HPI = Value)

# identify the mismatched districts by uncommenting 2 lines below
# unique(hpi.data$id[!hpi.data$id %in% unique(nepal.adm3.shp.df$id)])
# unique(unique(nepal.adm3.shp.df$id)[!unique(nepal.adm3.shp.df$id) %in% hpi.data$id])

# fix the mismatched districts
hpi.data$id[hpi.data$id == "Darchaula"] <- "Darchula"
hpi.data$id[hpi.data$id == "Kavre"] <- "Kavrepalanchok"

# merge hpi data with shapefile dataframe
shp_df <- shp_df %>% 
  left_join(hpi.data, by="id")

# finding the centers of all the district polygons. 
centroids <- rgeos::gCentroid(nepal_shp, byid = TRUE, id = unique(nepal_shp$NAME_3))
centroids_df <- centroids %>% 
  as.data.frame() %>% 
  mutate(label = row.names(.))

# showing district names for which HPI>40
centroids.selected <- centroids_df[centroids_df$label %in% (hpi.data[hpi.data$HPI>40,]$id),]

# setting up bare theme
theme_bare <- theme(
  axis.line = element_blank(), 
  axis.text.x = element_blank(), 
  axis.text.y = element_blank(),
  axis.ticks = element_blank(), 
  axis.title.x = element_blank(), 
  axis.title.y = element_blank(),
  legend.text=element_text(size=7),
  legend.title=element_text(size=8),
  panel.background = element_blank(),
  panel.border = element_rect(colour = "gray", fill=NA, size=0.5)
)

# plot map
ggplot(data = shp_df,aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = HPI), color = 'gray', size = 0.1) +
  ggtitle("Human Poverty Index Map") +
  guides(fill=guide_colorbar(title="HP Index")) + 
  scale_fill_gradient(high = "#e34a33", low = "#fee8c8", guide = "colorbar") +
  coord_fixed(1.3) +
  theme(legend.justification=c(0,-0.1), legend.position=c(0.05,0)) +
  with(centroids.selected, annotate(geom="text", x = x, y = y, label=label, size=2)) +
  theme_bare

map.choropleth.nepal.R

library(rgdal)
library(ggplot2)
library(dplyr)

# clone NepalMaps from https://github.com/anjesh/NepalMaps 

# read shapefile
nepal.adm3.shp <- readOGR(dsn="./NepalMaps/baselayers/NPL_adm", layer="NPL_adm3", stringsAsFactors = FALSE)

# fortify shapefile data to data frame
nepal.adm3.shp.df <- fortify(nepal.adm3.shp, region = "NAME_3")

# write districts to csv file
nepal.adm3.shp.df %>%
  distinct(id) %>%
  write.csv("districts.csv", row.names = FALSE)

# update the HPI data in districts.csv as per http://data.opennepal.net/content/human-poverty-index-value-districts-2011
hpi.data <- read.csv("districts.csv")
colnames(hpi.data) <- c("id","HPI")

# merge hpi data with shapefile dataframe
nepal.adm3.shp.df <- merge(nepal.adm3.shp.df, hpi.data, by ="id")

# finding the centers of all the district polygons. Thanks http://stackoverflow.com/questions/28962453/how-can-i-add-labels-to-a-choropleth-map-created-using-ggplot2
centroids <- setNames(do.call("rbind.data.frame", by(nepal.adm3.shp.df, nepal.adm3.shp.df$group, function(x) {Polygon(x[c('long', 'lat')])@labpt})), c('long', 'lat')) 
centroids$label <- nepal.adm3.shp.df$id[match(rownames(centroids), nepal.adm3.shp.df$group)]

# showing district names for which HPI>40
centroids.selected <- centroids[centroids$label %in% (hpi.data[hpi.data$HPI>40,]$id),]

# setting up bare theme
theme_bare <- theme(
  axis.line = element_blank(), 
  axis.text.x = element_blank(), 
  axis.text.y = element_blank(),
  axis.ticks = element_blank(), 
  axis.title.x = element_blank(), 
  axis.title.y = element_blank(),
  legend.text=element_text(size=7),
  legend.title=element_text(size=8),
  panel.background = element_blank(),
  panel.border = element_rect(colour = "gray", fill=NA, size=0.5)
)

# plot map
ggplot(data = nepal.adm3.shp.df,aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = HPI), color = 'gray', size = 0.1) +
  ggtitle("Human Poverty Index Map") +
  guides(fill=guide_colorbar(title="HP Index")) + 
  scale_fill_gradient(high = "#e34a33", low = "#fee8c8", guide = "colorbar") +
  coord_fixed(1.3) +
  theme(legend.justification=c(0,-0.1), legend.position=c(0.05,0)) +
  with(centroids.selected, annotate(geom="text", x = long, y = lat, label=label, size=2)) +
  theme_bare

You might need to install extra libraries as well.

install.packages('maptools')
install.packages('rgeos')
install.packages('rgdal')