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

Hi,
Glad to have found out the long needed tutorial.

Few things I am not able to understand:

• I am confused for line 14 to 23.
• In line 19, is it the same districts.csv data you are reading that you wrote in line 16 ? If its the case, I don't think we can give 2 column names for data with 1 column.
• I am not able to understand the merge step.

Given the increasing importance of .geojson file now, another tutorial with this format would be helpful.

while using line 26 and 27 produces errors message and deleting 75 districts records but did indirectly
dd=read.csv(("C:/Users/yagyarimal/Desktop/publication/NepalMaps-master/districts.csv"))
centroids=data.frame(centroids,dd)
colnames(centroids)=c("long","lat","label")
head(centroids)

I have simplified the script https://gist.github.com/anjesh/64cd61e3f94352ce5249ef368082df0c#file-map-choropleth-nepal-updated-r You should be able to make changes to the code to make the centroids work as per your need

Thankx a lot Anjesh ! I am a beginner and am self-learning and this helps me a lot.

Thanks so much Anjesh. I worked through the steps and I got it right to a point where I got stocked.. I am working on Africa and I am only selecting 20 countries out of Africa.. I used a percapita GDP of these countries for a particular year and I set the other countries which is not among my study to zero.. How ever when I tried to plot the cloroplethr with label name.. It couldn't imprint names on the countries .. Showing geom_text deleted 768 missing and it returned to the normal map without labels..

You might need to install extra libraries as well.

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