###########################################################################################
##                                                                                       ##
## R code for pre-processing database downloads for Malawi population density 		 ##
## clickable example application at http://sysbio.mrc-bsu.cam.ac.uk/click_dynmap_example ##
##                                                                                       ##
###########################################################################################

library(raster)

##################################### GRUMP ##############################################

## Downloaded data details:-
## ncols         600
## nrows         1200
## xllcorner     32
## yllcorner     -18
## cellsize      0.0083333333333
## NODATA_value  -9999

## read in matrix

grump.mwi <- as.matrix(read.table(file="mwiuds00ag.asc", skip=6))  

## set nodata to NA 

grump.mwi[which(grump.mwi==-9999)] <- NA 

## create raster object with the correct WGS84 projection and extent

grump.mwi.utm <-raster(as.matrix(grump.mwi), crs="+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0", xmn=32, ymn=-18, xmx=32+0.0083333333333*600, ymx=-18+0.0083333333333*1200)

## save WGS84 raster

save(grump.mwi.utm, file="grump_mwi_utm.Rdata")

## project WGS84 raster to Google projection

grump.mwi.goog <- projectRaster(grump.mwi.utm,  crs="+init=epsg:3857")

## save Google raster

save(grump.mwi.goog, file="grump_mwi_goog.Rdata")


##################################### GPW ################################################

## Downloaded data details:-
## ncols         120
## nrows         240
## xllcorner     32
## yllcorner     -18
## cellsize      0.0416666666667
## NODATA_value  -9999

## read in matrix

gpw.mwi <- as.matrix(read.table(file="mwids00ag.asc", skip=6))

## set nodata to NA

gpw.mwi[which(gpw.mwi==-9999)] <- NA

## create raster object with the correct WGS84 projection and extent

r <- raster(as.matrix(gpw.mwi), crs="+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0", xmn=32, ymn=-18, xmx=32+0.0416666666667*120, ymx=-18+0.0416666666667*240)

## each GPW raster cell is 5x5 GRUMP cells so split each GPW cell into 25 identical cells

r.dis<-disaggregate(r, fact=5)
gpw.mwi.utm <- r.dis

## save WGS84 raster

save(gpw.mwi.utm, file="gpw_mwi_utm.Rdata")

## project WGS84 raster to Google projection

gpw.mwi.goog <- projectRaster(r.dis,  crs="+init=epsg:3857")

## save Google raster

save(gpw.mwi.goog, file="gpw_mwi_goog.Rdata")


##################################### Afripop ############################################

### File "apmwi10v4.flt" was opened in SAGA-GIS using the gdal module and then saved as an xyz file "apmwi.xyz"

## read xyz file

apmwi<-read.table(file="apmwi.xyz")

## remove nodata entries

apmwi<-apmwi[which(apmwi[,3]!=-9999),]

## create raster from xyz data

r<-rasterFromXYZ(apmwi)

## Afripop has a 10x10 higher resolution than GRUMP so aggregate cells to reduce to same resolution as GRUMP

afpop.mwi.utm.ag <- aggregate(r, fact=10, fun=sum)

## set the projection as WGS84

projection(afpop.mwi.utm.ag)<-"+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

## project raster to Google projection

afpop.mwi.goog.ag <- projectRaster(afpop.mwi.utm.ag,  crs="+init=epsg:3857")

## xy limits of Afripop raster are not the same as for GRUMP and GPW so resample onto the same area as these two datasets

afpop.mwi.goog.res <- resample(afpop.mwi.goog.ag, grump.mwi.goog)
afpop.mwi.utm.res <- resample(afpop.mwi.utm.ag, grump.mwi.utm)

## save Google raster

save(afpop.mwi.goog.res, file="afpop_mwi_goog_res.Rdata")

## save WGS84 raster

save(afpop.mwi.utm.res, file="afpop_mwi_utm_res.Rdata")