Cartogram for disability and internet access in Indonesia

#choro-carto-prov-kab.R
library(tmap)
library(maptools)
library(ggplot2)
library(rgeos)
library(foreign)
library(Rcartogram)
library(getcartr)

setwd('~/_codeIDNgeo')
idnprov <- readRDS('IDN_adm1.rds')  #idnkab: IDN_adm2.rds  #plot(idnprov)

disab2015 <- read.dta('~/_negeriku/tmpdisabbyidgeo.dta')
  #prepared by ana-dddisab.do ...

net2015   <- read.dta('~/_negeriku/tmpnetbyidgeo.dta')
  # ...which generated idgeo

idnprov@data = data.frame(idnprov@data,
  disab2015[match(idnprov@data[, "ID_1"], disab2015[, "idgeo"]), ])

idnprov@data = data.frame(idnprov@data, 
  net2015[match(idnprov@data[, "ID_1"], net2015[, "idgeo"]), ])

summary(idnprov$disab)
summary(idnprov$net)

disab.carto <-quick.carto(idnprov,idnprov$disab, blur=1)
spplot(disab.carto, "disab", col.regions=brewer.pal(3, "Greens"),
  main = "Cartogram of disability")

net.carto   <-quick.carto(idnprov,idnprov$net, blur=1)
spplot(idnprov,   "netcat", col.regions=brewer.pal(3, "Greens"),
  main = "Original projection")
spplot(net.carto, "netcat", col.regions=brewer.pal(3, "Greens"),
  main = "Cartogram of internet access")

pdf('~/_negeriku/net2015.pdf',width=6,height=4,paper='special') 
  par(mfrow=c(2,1),mar=c(0.5,0.5,3,0.5))
  plot(idnprov, main = "Original Projection")
  plot(net.carto, main = "Cartogram of internet access ")
dev.off()

####################################PROVINSI LAGI
idn <- rgdal::readOGR('IDN_adm1.shp')
#plot(idn)  #periksa: bagus.

idnY <- read.dta('~/_negeriku/kerjafml2014.dta')
  #theme: ID_1, formal job, prov
idnY$ID_1 <- as.factor(idnY$ID_1)   #prov 34 values

idnm <- append_data(idn, idnY, key.shp = 'ID_1', key.data = 'ID_1',
                    ignore.na = TRUE, ignore.duplicates = TRUE)
#map data, ID_1 34 values!

choromap <- qtm(idnm, fill='Formal', 
  legend.position = c("left", "bottom"))
choromap

#sekarang buat kartogram
idn.carto <- quick.carto(idnm, idnm@data$Formal, blur=0)
  #bbrp menit kemudian...
#quick.carto doesn't like exactly zero values; naturally!
#  But it can cope with missing!

idn.buffer <- gBuffer(idn.carto, width=0, byid=TRUE)
idn.f <- fortify(idn.buffer, region="ID_1") 
  #ID_1 -> id; convert to a regular data frame: KABU or ID2014_2
idn.f$id <- as.factor(idn.f$id)  
  #my peace of mind that 'id' and 'ID_1" same-typed
idn.f <- merge(idn.f, idnm@data, by.x="id", by.y="ID_1") 
  #merge cartogram with map

carto.fjob <- ggplot(idn.f, aes(long, lat, group=group,
    fill=idn.f$Formal)) + 
  labs(fill="Formal jobs") + geom_polygon() + 
  ggtitle("Cartogram formal jobs 2014")
carto.fjob

#KABUPATEN-KOTA
idn <- rgdal::readOGR('petabps2014.shp')   
#plot(idn)  #periksa: bagus.

idnY  <- read.dta('~/_negeriku/ssn14rawat.dta') #theme var
idnY$KABU <- as.factor(idnY$KABU)   #KABU 497 values
summary(idnY)
  #catch zeroes because quick.carto abhors them.
  # Better set it to missing.
idnY2 <- read.table('~/_negeriku/cataskabu50.tsv', 
   sep='\t', header=TRUE)  #another theme var
idnY2$cdkabkota <- as.factor(idnY2$cdkabkota)   # 497 values
summary(idnY2)  
  #catch zeroes because quick.carto abhors 'em.
  # Drop 'em or give 'em an average value
idnY2 <- idnY2[idnY2$catas>0,]

idnY <- merge(idnY, idnY2, by.x='KABU', by.y='cdkabkota')

idnm <- tmaptools::append_data(idn, idnY, 
           key.shp = 'ID2014_2', key.data = 'KABU',
           ignore.na = TRUE, ignore.duplicates = TRUE) 
#map data, ID2014_2 514 values!
#Please check whether the variables are renamed with .x or .y

choromap <- qtm(idnm, fill='rj')  #peta tema
choromap #takes time!

chorocatas <- qtm(idnm, fill='catas')  #gambar peta tema (with renamed var? CHECK THIS NAME!)
chorocatas #takes time!

#sekarang buat kartogram
idn.carto <- quick.carto(idnm, idnm@data$rj, blur=0)  
  #bbrp menit kemudian...
  #quick.carto doesn't like exactly zero values; naturally!
  # But it can cope with missing!

idn.buffer <- gBuffer(idn.carto, width=0, byid=TRUE)  
  ##47 warnings for missing values!
idn.f <- fortify(idn.buffer, region="ID2014_2")  
  #convert to a regular data frame: KABU or ID2014_2
idn.f <- merge(idn.f, idnm@data, by.x="id", by.y="ID2014_2") 
  #merge cartogram with map

cartorj <- ggplot(idn.f, aes(long, lat, group=group, 
     fill=idn.f$rj)) + 
  labs(fill="Rawat jalan") + geom_polygon() + 
  ggtitle("Cartogram rawat jalan 2014")
cartorj #neat!

### buat kartogram catastrophic health expenditure
idn.carto2 <- quick.carto(idnm, idnm@data$catas, blur=0)  
  #bbrp menit kemudian...
idn.buffer2 <- gBuffer(idn.carto2, width=0, byid=TRUE)  
  ##47 warnings for missing values
idn.f2 <- fortify(idn.buffer2, region="ID2014_2")  
  #convert to a regular data frame: KABU or ID2014_2
idn.f2 <- merge(idn.f2, idnm@data, by.x="id", by.y="ID2014_2") #merge cartogram with map
cartocatas <- ggplot(idn.f2, aes(long, lat, group=group,
     fill=idn.f2$catas)) + 
  labs(fill="Catastrophic health exp") + geom_polygon() +
  ggtitle("Catastrophic health exp")
cartocatas #neat!