Install Packages

install.packages(“sf”) install.packages(“raster”) install.packages(“spData”) install.packages(“spDataLarge”) install.packages(“tmap”) install.packages(“mapview”) install.packages(“shiny”) install.packages(“leaflet”)

Load Packages

library(sp)
library(sf)
library(tmap)
library(raster)
library(spData)
library(spDataLarge)
library(leaflet)
library(tidyverse)
library(haven)
library(rsconnect)

Load Data

# Replication Esther Duflo
duflo2001 <- read_dta("G:/My Drive/Dataset/Susenas/dta/additional/duflo2001.dta")

# SD Inpres Data, unique by ROB
inpres <- read_dta("G:/My Drive/Dataset/Susenas/dta/additional/esther_data.dta")

# Indonesia Map 1995 (Source IPUMS International)
kab90 <- st_read("G:/My Drive/Dataset/Susenas/dta/map/1990/geo2_id1990.shp")
## Reading layer `geo2_id1990' from data source `G:\My Drive\Dataset\Susenas\dta\map\1990\geo2_id1990.shp' using driver `ESRI Shapefile'
## Simple feature collection with 298 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 95.00708 ymin: -10.99251 xmax: 141.0222 ymax: 6.08125
## geographic CRS: WGS 84
as_tibble(kab90)
## # A tibble: 298 x 7
##    CNTRY_NAME ADMIN_NAME CNTRY_CODE IPUM1990 REGY1990 PARENT
##    <chr>      <chr>      <chr>      <chr>    <chr>    <chr> 
##  1 Indonesia  Lake Toba  360        088088   8888     088   
##  2 Indonesia  Aceh Sela~ 360        011001   1101     011   
##  3 Indonesia  Aceh Teng~ 360        011002   1102     011   
##  4 Indonesia  Aceh Timur 360        011003   1103     011   
##  5 Indonesia  Aceh Teng~ 360        011004   1104     011   
##  6 Indonesia  Aceh Barat 360        011005   1105     011   
##  7 Indonesia  Aceh Besar 360        011006   1106     011   
##  8 Indonesia  Pidie      360        011007   1107     011   
##  9 Indonesia  Aceh Utara 360        011008   1108     011   
## 10 Indonesia  Kodya. Sa~ 360        011072   1172     011   
## # ... with 288 more rows, and 1 more variable: geometry <MULTIPOLYGON [°]>
# Delete Timor Leste
kab90 <- subset(kab90, kab90$PARENT!="054")
  
# Create ROB
kab90 <- kab90 %>%
  mutate(ROB=kab90$REGY1990)

# Outer join (Merge, keep all)
inpreskab <- st_as_sf(merge(x=inpres, y=kab90, by="ROB", all=TRUE))

Choropleth of SD Inpres Intensity by District

# Indonesia Map
leaflet(data = kab90) %>%
  addPolygons(label = ~ADMIN_NAME,
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE))
# Interactive Choropleth with Program intensity popups
# Add Color Palette based on predefined point of a variable's vaue
bins <- c(0,2,4,6,8,Inf)
pal <- colorBin("YlOrRd", domain = inpres$nin, bins = bins)

# Map with program intensity as its popup
inpreskab %>%
  mutate(popup = str_c("<strong>", ADMIN_NAME, "</strong>",
                       "<br/>",
                       "Program Intensity: ", nin) %>%
           map(htmltools::HTML)) %>%
  leaflet() %>%
  addProviderTiles("CartoDB") %>%
  addPolygons(label = ~popup,
              fillColor = ~pal(nin),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "15px",
                direction = "auto")) %>%
  addLegend(pal = pal,
            values = ~n,
            opacity = 0.7,
            title = NULL,
            position = "bottomright")
## Warning in is.na(values): is.na() applied to non-(list or vector) of type
## 'closure'