#######################
### SETUP WORKSPACE ###
#######################

library(maptools)
library(raster)
library(spatstat)
library(rgdal)
library(rgeos)
library(plyr)
library(cshapes)

## Clear the workspace
rm(list=ls())

##############################################################
### NOTE: This line should be changed to point to the folder
###       where you have unzipped the contents of the
###       replication archive.
##############################################################

## Set working directory
setwd("C:/warren_jpr2015_rep")

###################################
### SET GLOBAL PARAMETER VALUES ###
###################################

#Set year for polygon-based data 
shapeyear <- 2005

#Set projection parameters
#(Lambert Azimuthal Equal Area)
proj.string <- "+proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"


#########################
### DEFINE DATA FILES ###
#########################

########################################################
### NOTE: All users of the data files referenced herein
###       should consult "citations.txt" for detailed
###       instructions.
########################################################

# #GECON file name
# gecon.fname <- "gecon40.csv"
# 
# #EPR file name
# epr.fname <- "epr-eth_v.2.0.csv"
# 
# #Geo-EPR file name
# epr.geo.fname <- "GeoEPR-ETH_2.0.shp"
# 
# #GED data file name
# ged.fname <- "ucdp-ged15.csv"
# 
# #Load external data files
# ged <- read.csv(ged.fname, stringsAsFactors=FALSE)
# gecon <- read.csv(gecon.fname, stringsAsFactors=FALSE, comment.char="#")
# epr.data <- read.csv(epr.fname, stringsAsFactors=FALSE)
# epr.geo <- readShapePoly(epr.geo.fname, delete_null_obj=TRUE)
# 
# #Save as compressed archive
# save.image("ext_dat.RData")


## Load external data sources from archive
load("ext_dat.RData")


## Load cluster-level DHS survey data file
dhs <- read.csv("dhs_cluster_summary.csv", stringsAsFactors=FALSE)


#########################################
### Generate Raster Surface Variables ###
#########################################

#Build table of ASCII raster variable names
vartab <- matrix(nrow=0, ncol=2)
vartab <- rbind(vartab, c('elev',    'af_gtopo.asc'))
vartab <- rbind(vartab, c('popdens', 'af_pdens.asc'))
vartab <- rbind(vartab, c('urban',   'af_urext.asc'))

#Load ASCII raster data
#Rescale to common template
#Convert to '.im' format
for (row in c(1:nrow(vartab))) {  
  newvarname <- paste0(vartab[row, 1], '.im')
  cat(paste0('Creating [', newvarname, '] \n'))
  dat <- readAsciiGrid(vartab[row, 2], 
                       as.image = TRUE, 
                       proj4string = CRS("+init=epsg:4326"))
  
  #Define raster template 
  #(using the first file from the raster table)
  if (row == 1) {
    dat.im <- as.im(dat)
    r <- raster(nrow=dat.im$dim[1], ncol=dat.im$dim[2])
    extent(r) <- extent(dat.im$xrange, dat.im$yrange)
    r <- projectExtent(r, proj.string)
    res(r) <- 5000
  }
  
  #Resample raster using template
  dat <- raster(dat)
  dat <- projectRaster(dat, r, method='ngb')
  
  #Convert to image via SpatialGridDataFrame
  dat <- as(dat, 'SpatialGridDataFrame')
  assign(newvarname, as.im.SpatialGridDataFrame(dat))
}
rm(vartab, newvarname, row, dat, dat.im)


#Recode urbanization values from (1,2) to (0,1)
urban.im <- eval.im(urban.im - 1)


#########################################
### Define survey years and countries ###
#########################################

#Load country shapes data
cshapes <- cshp()

#Re-project (using equal area projection)
cshapes <- spTransform(cshapes, CRS(proj.string))

#Limit country shapes to single year
cshapes$cowcode <- cshapes$COWCODE
cshapes <- cshapes[cshapes$COWSYEAR <= shapeyear & 
                   cshapes$COWEYEAR >= shapeyear, ]

#Extract unique survey country codes
dhs.cowcode <- sort(unique(dhs[, 'cowcode']))

#Extract unique survey country-years
dhs.cowcode.year <- unique(dhs[, c('dhsyear', 'cowcode')])

#Country-years w/o 2nd Uganda survey round
dhs.cy.min <- dhs.cowcode.year[dhs.cowcode.year$dhsyear <= 2010,]

#Extract shapes of surveyed countries
dhs.cshapes <- merge(cshapes, dhs.cy.min, by="cowcode", all=FALSE)

#Extract shape of Uganda
uganda.cshape <- cshapes[cshapes$cowcode == 500,]

#Convert survey locations to SpatialPointsDataFrame
coordinates(dhs) <- ~ longnum + latnum
proj4string(dhs) <- CRS("+init=epsg:4326")

#Re-project (using equal area projection)
dhs <- spTransform(dhs, CRS(proj.string))



#####################################
### Create Gridded GDP Per Capita ###
#####################################

#Limit GECON data to Africa
gecon <- gecon[gecon$LONGITUDE >= -20 & 
               gecon$LONGITUDE <=  55 & 
               gecon$LAT       >= -38 & 
               gecon$LAT       <=  38 , ]

var_stub <- "PPP"
newvarname <- "ppp_est"

#Start with 1990 values
varname <- paste0(var_stub, "1990", "_40") 
gecon[[newvarname]] <-  gecon[[varname]]

#Fill in more recent values when available
for (yr in c("1995", "2000", "2005")) {
  varname <- paste0(var_stub, yr, "_40") 
  gecon[[newvarname]] <- ifelse(
    gecon[[varname]] > 0 & is.na(gecon[[varname]])==FALSE, 
    gecon[[varname]], 
    gecon[[newvarname]])
}

#Convert to image (*.im) format
cat('Creating [gdp_ppp.im] \n')
grid <- gecon[, c('LAT', 'LONGITUDE', newvarname)]
coordinates(grid) <- ~ LONGITUDE + LAT 
gridded(grid) <- TRUE
grid <- as(grid, "SpatialGridDataFrame")
proj4string(grid) <- CRS("+init=epsg:4326") 
grid <- raster(grid)
grid <- projectRaster(grid, r, method='ngb')
grid <- as(grid, "SpatialGridDataFrame")
gdp_ppp.im <- as.im.SpatialGridDataFrame(grid)

rm(var_stub, varname, newvarname, yr, grid, gecon)


###########################################
### Convert EPR Polygons to Raster Data ###
###########################################

#Define initial projection
proj4string(epr.geo) <- CRS("+init=epsg:4326")

#Re-project (using equal area projection)
epr.geo <- spTransform(epr.geo, CRS(proj.string))

#Limit to single year
epr.data <- epr.data[epr.data$from <= shapeyear & 
                     epr.data$to   >= shapeyear, ]

epr.geo  <- epr.geo[epr.geo$STARTYEAR <= shapeyear & 
                    epr.geo$ENDYEAR   >= shapeyear, ]


#Extract geo-referenced groups
#(see GEO-EPR codebook for details)
epr.geo <- epr.geo[epr.geo$TYPE %in% c(1,3,6), ]
epr.data$COWGROUP <- epr.data$cowgroupid

#Limit to surveyed countries
epr.geo  <- epr.geo[epr.geo$COWCODE %in% dhs.cowcode, ]

#Merge EPR and GEO-EPR
epr.geo.data <- merge(epr.geo, 
                      epr.data,
                      by="COWGROUP", all=FALSE)

#Recode EPR status categories
epr.geo.data$shrpower <- 0
epr.geo.data$shrpower[epr.geo.data$status %in% 
                        c("SENIOR PARTNER",
                          "JUNIOR PARTNER")] <- 1

#Convert EPR polygons to raster image (*.im)
cat('Creating [shrpowwr.im] \n')
grid <- rasterize(epr.geo.data, r, 'shrpower', fun=max)
grid <- as(grid, "SpatialGridDataFrame")
grid <- as.im.SpatialGridDataFrame(grid) 
grid[eval.im(grid < 0)] <- 0
grid[eval.im(is.na(grid))] <- 0
shrpower.im <- grid
rm(epr.geo, epr.data, epr.geo.data, grid)


#####################################
### Generate COWCODE Raster Image ###
#####################################

cat('Creating [cowcode.im] \n')
cowcode.im <- rasterize(dhs.cshapes, r, 'cowcode')
cowcode.im <- as(cowcode.im, "SpatialGridDataFrame")       
cowcode.im <- as.im.SpatialGridDataFrame(cowcode.im)  


####################
### Save Results ###
####################

save.image(file="load_data.RData")


