Data Preparation

We obtained the geographical coordinates of Nepal’s administrative units (District Level: 77 Districts) in GeoJSON format from the Open Knowledge Nepal data portal (http://localboundries.oknp.org/).

Define Libraries

library("stringr")
library("dplyr")
library("reshape2")
library("geojsonio")
library("broom")
library("ggplot2")
library("ggthemes")
library("RColorBrewer")

Define Path

dir.wrk <- getwd()
dir.data <- file.path(dir.wrk, "data/data_household")
dir.annot <- file.path(dir.wrk, "data/data_annotations")
dir.output <- file.path(dir.wrk, "data/data_processed")
dir.maps <- file.path(dir.wrk, "data/data_maps")

Define Files

file.geo <- file.path(dir.maps, "nepal_district.geojson")
file.dat1 <- file.path(dir.output, "maps_tbl_district_total_household.tsv")
file.dat2 <- file.path(dir.output, "maps_tbl_district_fueltype_ratio.tsv")
file.dat3 <- file.path(dir.output, "maps_tbl_district_ethnicity_ratio.tsv")

Load Household Population by District

dat1 <- read.delim(file.dat1, header = TRUE, stringsAsFactors = FALSE)

head(dat1)
##               id  Freq
## 1        Dhading 86345
## 2        Dolakha 70495
## 3         Gorkha 75883
## 4 Kavrepalanchok 91895
## 5      Makwanpur 88365
## 6        Nuwakot 75429

Load Maps

geo <- geojsonio::geojson_read(file.geo, parse=TRUE, what="sp")
geo77 <- broom::tidy(geo, region="DISTRICT")
geo11 <- subset(geo77, geo77$id %in% dat1$id)

head(geo11)
## # A tibble: 6 x 7
##    long   lat order hole  piece group     id     
##   <dbl> <dbl> <int> <lgl> <fct> <chr>     <chr>  
## 1  84.7  27.9 68752 FALSE 1     Dhading.1 Dhading
## 2  84.7  27.9 68753 FALSE 1     Dhading.1 Dhading
## 3  84.7  27.9 68754 FALSE 1     Dhading.1 Dhading
## 4  84.7  27.9 68755 FALSE 1     Dhading.1 Dhading
## 5  84.7  27.9 68756 FALSE 1     Dhading.1 Dhading
## 6  84.7  27.9 68757 FALSE 1     Dhading.1 Dhading

Plot Choropleth Maps: Nepal with all Districts

# PREPARE DATA ---
d <- geo77
d$Status <- 0
d$Status[which(d$id %in% dat1$id)] <- 1
d$Status <- as.factor(d$Status)

# DEFINE COLORS ---
cpalette.grp <- c("#FFFFFF","#BDBDBD")

# PLOT
map77 <- ggplot(data=d, aes(x=long, y=lat, group=group)) +
          geom_path() +
          geom_polygon(aes(fill=Status), color="#000000") +
          scale_fill_manual(values=cpalette.grp) +
          coord_equal() +
          theme_map() +
          theme(
            plot.title = element_text(size = 10, color="#000000", hjust=0.5),
            #aspect.ratio = 1,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            axis.ticks = element_blank(),   
            strip.text = element_text(size = 10, color="#000000", hjust=0.5),
            strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
            panel.background = element_rect(fill="#FFFFFF", color=NA),
            legend.text = element_text(size = 10, color="#000000"),
            legend.title = element_blank(),
            legend.key.size = unit(0.5, "cm"),
            legend.position = "none") +
          guides(fill=guide_legend(title="No. of Households"))

map77

Load Household Population by Earthquake-affected Districts

dat1 <- read.delim(file.dat1, header = TRUE, stringsAsFactors = FALSE)

# CONVERT TO DISCRETE DATA ---
dm <- mutate(dat1, Group = cut(Freq, seq(10000, 1e+05, 10000), include.lowest = TRUE)) %>% 
    group_by(Group)
group_ratio <- levels(dm$Group)

# MERGE DATA ---
df <- merge(geo11, dm, by = "id")
head(df)
##        id     long      lat order  hole piece     group  Freq         Group
## 1 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1 86345 (8e+04,9e+04]
## 2 Dhading 84.73970 27.94059 68753 FALSE     1 Dhading.1 86345 (8e+04,9e+04]
## 3 Dhading 84.73969 27.94087 68754 FALSE     1 Dhading.1 86345 (8e+04,9e+04]
## 4 Dhading 84.73980 27.94221 68755 FALSE     1 Dhading.1 86345 (8e+04,9e+04]
## 5 Dhading 84.73977 27.94336 68756 FALSE     1 Dhading.1 86345 (8e+04,9e+04]
## 6 Dhading 84.73981 27.94370 68757 FALSE     1 Dhading.1 86345 (8e+04,9e+04]

Plot Choropleth Maps: Total No. of Households

map1 <- ggplot(data=df, aes(x=long, y=lat, group=group, fill=Group)) +
          geom_path() +
          geom_polygon(aes(fill=Group), color="#000000") +
          scale_fill_brewer(palette = "Blues") +
          coord_equal() +
          theme_map() +
          theme(
            plot.title = element_text(size = 10, color="#000000", hjust=0.5),
            aspect.ratio = 1,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            axis.ticks = element_blank(),   
            strip.text = element_text(size = 10, color="#000000", hjust=0.5),
            strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
            panel.background = element_rect(fill="#FFFFFF", color=NA),
            legend.text = element_text(size = 10, color="#000000"),
            legend.title = element_text(size = 10, color="#000000"),
            legend.key.size = unit(0.5, "cm"),
            legend.position = "right") +
          guides(fill=guide_legend(title="No. of Households", ncol=1))

map1

Load Household FuelType Usage Proportions by Earthquake affected Districts

dat2 <- read.delim(file.dat2, header = TRUE, stringsAsFactors = FALSE)

head(dat2)
##               id Electricity   Gobar.Gas   Kerosene     LP.Gas     Others
## 1        Dhading  0.08465608 0.205064154 0.05714286 0.10945537 0.06768559
## 2        Dolakha  0.14285714 0.005109572 0.11904762 0.08824108 0.13537118
## 3         Gorkha  0.08994709 0.174520268 0.07619048 0.09196958 0.10698690
## 4 Kavrepalanchok  0.07936508 0.199500397 0.15238095 0.19404179 0.04148472
## 5      Makwanpur  0.08994709 0.227546270 0.11428571 0.24897251 0.16157205
## 6        Nuwakot  0.11111111 0.045304871 0.09523810 0.07829843 0.05676856
##         Wood
## 1 0.11518469
## 2 0.09625774
## 3 0.10178191
## 4 0.11323104
## 5 0.10059783
## 6 0.10454018
# CONVERT TO DISCRETE DATA ---
dm <- reshape2::melt(dat2, id.vars = "id", variable.name = "FuelType", value.name = "Ratio")
dm <- mutate(dm, RatioGroup = cut(Ratio, seq(0, 0.25, 0.05), include.lowest = TRUE)) %>% 
    group_by(RatioGroup)

group_ratio <- levels(dm$RatioGroup)
type_fuel <- c("Wood", "LP.Gas", "Gobar.Gas", "Kerosene", "Electricity", "Others")

# MERGE DATA ---
df <- merge(geo11, dm, by = "id")

# FACTORIZE DATA ---
df$FuelType <- factor(df$FuelType, levels = type_fuel)
df$RatioGroup <- factor(df$RatioGroup, levels = group_ratio)

head(df)
##        id     long      lat order  hole piece     group    FuelType      Ratio
## 1 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1 Electricity 0.08465608
## 2 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1      Others 0.06768559
## 3 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1        Wood 0.11518469
## 4 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1    Kerosene 0.05714286
## 5 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1      LP.Gas 0.10945537
## 6 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1   Gobar.Gas 0.20506415
##   RatioGroup
## 1 (0.05,0.1]
## 2 (0.05,0.1]
## 3 (0.1,0.15]
## 4 (0.05,0.1]
## 5 (0.1,0.15]
## 6 (0.2,0.25]

Plot Choropleth Maps: Proportion of Households using six FuelTypes

map2 <- ggplot(data=df, aes(x=long, y=lat, group=group, fill=RatioGroup)) +
          geom_path() +
          geom_polygon(aes(fill=RatioGroup), color="#000000") +
          scale_fill_brewer(palette = "Blues") +
          facet_wrap(~FuelType, ncol=3, drop=FALSE,  strip.position="top") +
          coord_equal() +
          theme_map() +
          theme(
               plot.title = element_text(size = 10, color="#000000", hjust=0.5),
               aspect.ratio = 1,
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),    
               strip.text = element_text(size = 10, color="#000000", hjust=0.5),
               strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
               panel.background = element_rect(fill="#FFFFFF", color=NA),
               legend.text = element_text(size = 10, color="#000000"),
               legend.title = element_text(size = 10, color="#000000"),
               legend.key.size = unit(0.5, "cm"),
               legend.position = "right") +
          guides(fill=guide_legend(title="Proportion of Households", ncol=1))

map2

Load Household Population by Ethnicity and District

dat3 <- read.delim(file.dat3, header = TRUE, stringsAsFactors = FALSE)

head(dat3)
##               id    Brahman Chepang.Thami    Chettri      Dalit Gurung.Magar
## 1        Dhading 0.13183869  1.992073e-01 0.09895938 0.13744536   0.18856305
## 2        Dolakha 0.06860788  3.570013e-01 0.18636705 0.06896213   0.02672287
## 3         Gorkha 0.11886990  5.221076e-02 0.06745067 0.17466477   0.30926197
## 4 Kavrepalanchok 0.19493279  8.884029e-04 0.10279011 0.10281694   0.04920577
## 5      Makwanpur 0.12246981  2.473860e-01 0.07565513 0.05270396   0.06055718
## 6        Nuwakot 0.13789929  6.833869e-05 0.07101362 0.06762366   0.04849707
##      Madhesi Muslim.Others      Newar   Rai.Limbu     Tamang
## 1 0.06368421    0.09718067 0.11869743 0.030502227 0.08889155
## 2 0.02500000    0.02683347 0.10158448 0.170972201 0.05719957
## 3 0.12394737    0.31012601 0.08608953 0.039318077 0.01146000
## 4 0.03947368    0.08258544 0.16717928 0.007617877 0.15064451
## 5 0.27210526    0.05439217 0.08283894 0.072277684 0.20192754
## 6 0.08394737    0.06055661 0.08180400 0.105882353 0.14936267
# CONVERT TO DISCRETE DATA ---
dm <- reshape2::melt(dat3, id.vars = "id", variable.name = "Ethnicity", value.name = "Ratio")
dm <- mutate(dm, RatioGroup = cut(Ratio, seq(0, 0.4, 0.05), include.lowest = TRUE)) %>% 
    group_by(RatioGroup)

group_ratio <- levels(dm$RatioGroup)
type_ethnicity <- c("Newar", "Brahman", "Madhesi", "Chettri", "Muslim.Others", "Gurung.Magar", 
    "Dalit", "Rai.Limbu", "Tamang", "Chepang.Thami")

# MERGE DATA ---
df <- merge(geo11, dm, by = "id")

# FACTORIZE DATA ---
df$Ethnicity <- factor(df$Ethnicity, levels = type_ethnicity)
df$RatioGroup <- factor(df$RatioGroup, levels = group_ratio)

head(df)
##        id     long      lat order  hole piece     group     Ethnicity
## 1 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1       Brahman
## 2 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1 Muslim.Others
## 3 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1     Rai.Limbu
## 4 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1       Madhesi
## 5 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1        Tamang
## 6 Dhading 84.73974 27.94036 68752 FALSE     1 Dhading.1  Gurung.Magar
##        Ratio RatioGroup
## 1 0.13183869 (0.1,0.15]
## 2 0.09718067 (0.05,0.1]
## 3 0.03050223   [0,0.05]
## 4 0.06368421 (0.05,0.1]
## 5 0.08889155 (0.05,0.1]
## 6 0.18856305 (0.15,0.2]

Plot Choropleth Maps: Proportion of Households by Ethnicity

map3 <- ggplot(data=df, aes(x=long, y=lat, group=group, fill=RatioGroup)) +
          geom_path() +
          geom_polygon(aes(fill=RatioGroup), color="#000000") +
          scale_fill_brewer(palette = "Blues") +
          facet_wrap(~Ethnicity, ncol=4, drop=FALSE,  strip.position="top") +
          coord_equal() +
          theme_map() +
          theme(
               plot.title = element_text(size = 10, color="#000000", hjust=0.5),
               aspect.ratio = 1,
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),    
               strip.text = element_text(size = 10, color="#000000", hjust=0.5),
               strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
               panel.background = element_rect(fill="#FFFFFF", color=NA),
               legend.text = element_text(size = 10, color="#000000"),
               legend.title = element_text(size = 10, color="#000000"),
               legend.key.size = unit(0.5, "cm"),
               legend.position = "right") +
          guides(fill=guide_legend(title="Proportion of Households", ncol=1))

map3