# This R script was used to generate the building footprint metric layers
# in the following data release:

# Authors: Dooley, C. A. and Tatem, A.J. 
# Year: 2020. 
# Title: Gridded maps of building patterns throughout sub-Saharan Africa, version 1.0. 
# Publisher: University of Southampton: Southampton, UK. 
# Source data: Source of building Footprints "Ecopia Vector Maps Powered by Maxar Satellite Imagery"© 2020.
# doi:10.5258/SOTON/WP00666


# This script is intended as a guide to accompany the methods described in the data release
# If running this script, check file structure of building footprints data is appropriate for this code



rm(list=ls())

library(sf)         # version 0.8-1
library(rgdal)      # version 1.4-8
library(rgeos)      # version 0.5-2
library(data.table) # version 1.12.2
library(curl)       # version 3.3
library(raster)     # version 3.0-7

# set the working directory as the location of the building footprints data
setwd('ZZZ')

# output directory where the output rasters will be saved
output_filepath <- 'XXXX/XXX'

# directory where the country mastergrid can temporarily be saved
temp_filepath <- 'YYYY/YYY'

# list of countries with building footprints data
full <- list.dirs(path=paste(getwd()),full.names=FALSE,recursive=TRUE) 
gdbfolders <- full[grep(".gdb",full)]

# specify whether you want to create the 'urban' rasters
urban <- 'yes'


for(countryit in 1:length(gdbfolders)){
  
  # save country iso code
  ISO_code <- substring(gdbfolders[countryit],1,3)
  
  # list all feature classes in the file geodatabase
  fgdb <- gdbfolders[countryit]
  fc_list <- ogrListLayers(fgdb)
  
  # create centroids for each of the feature classes
  centroids_wgs <- list()
  for(i in 1:length(fc_list)){
    bf <- read_sf(dsn=fgdb,layer=paste(fc_list[i]))
    centroids_utm <- st_centroid(bf)
    rm(bf)
    centroids_wgs[[i]] <- st_transform(centroids_utm, crs ="+proj=longlat +datum=WGS84") # reproject
    rm(centroids_utm)
  }
  
  # merge centroid layers
  centroids_wgs_master <- do.call(rbind, centroids_wgs)
  rm(centroids_wgs)
  
  # read in master grid
  urlx <- paste("ftp://ftp.worldpop.org/GIS/Mastergrid/Global_2000_2020/",ISO_code,"/L0/",tolower(ISO_code),"_level0_100m_2000_2020.tif",sep='')
  utils::download.file(url = urlx,
                       destfile = paste(temp_filepath),
                       mode="wb",
                       quiet=FALSE,
                       method="libcurl")
  
  # create blank raster from master grid
  mgrid <- raster(paste(temp_filepath))
  mgrid[] <- NA
  
  # get cellID for each building
  points_IDs <- cellFromXY(mgrid, as(centroids_wgs_master, "Spatial"))
  
  # generate output rasters
  for(bf_char in c('Shape_Area','Shape_Length')){
    
    # create data table: a row for each building
    dt <- data.table(cellID = points_IDs, centroids_wgs_master[,bf_char])
    dt <- na.omit(dt)
    dt <- dt[,1:2]
    names(dt)[2] <- 'bchar'
    
    # sum building attribute for each cellID. Create and save raster
    tba <- dt[,.(tota = sum(bchar)),by=cellID]
    
    tbar <- mgrid
    tbar[tba$cellID] <- tba$tota
    if(bf_char == 'Shape_Area'){
      writeRaster(tbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_total_area.tif',sep=''))
    }
    if(bf_char == 'Shape_Length'){
      writeRaster(tbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_total_length.tif',sep=''))
    }
    rm(tbar,tba)
    
    # mean building attribute for each cellID. Create and save raster
    mba <- dt[,.(meana = mean(bchar)),by=cellID]
    
    mbar <- mgrid
    mbar[mba$cellID] <- mba$meana
    if(bf_char == 'Shape_Area'){
      writeRaster(mbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_mean_area.tif',sep=''))
    }
    if(bf_char == 'Shape_Length'){
      writeRaster(mbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_mean_length.tif',sep=''))
    }
    rm(mbar,mba)

    # cv building attribute for each cellID. Create and save raster
    cvba <- dt[,.(cva = sd(bchar)/mean(bchar)),by=cellID]
    cvba$cva[is.na(cvba$cva)] <- 0
    
    cvbar <- mgrid
    cvbar[cvba$cellID] <- cvba$cva
    if(bf_char == 'Shape_Area'){
      writeRaster(cvbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_cv_area.tif',sep=''))
    }
    if(bf_char == 'Shape_Length'){
      writeRaster(cvbar,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_cv_length.tif',sep=''))
    }
    rm(cvbar,cvba)
    
    
    if(bf_char == 'Shape_Area'){  # just want to do this section once - the fact that area is specified instead of length is irrelevant
      
      # sum buildings for each cellID. Create raster
      bc <- dt[,.N,by=cellID]
      bcr <- mgrid
      bcr[bc$cellID] <- bc$N
      
      # create px area raster
      pxarea <- area(mgrid)
      
      # create density raster: bc/px area
      bdensr <- bcr/pxarea
      
      # save building count and building density rasters
      writeRaster(bcr,filename = paste(output_filepath,'/',ISO_code,'_building_count.tif',sep=''))
      writeRaster(bdensr,filename = paste(output_filepath,'/',ISO_code,'_buildings_v1_0_density.tif',sep=''))
      
      
      if(urban == 'yes'){
        
        # identify 'clumps' (contiguous clusters of cells)
        clumps <- clump(bcr, directions=8)
        
        # data table for cells
        clumpdt <- data.table(clump_ID = clumps[],
                              bfcount = bcr[])
        clumpdt <- clumpdt[!is.na(clump_ID),]
        
        # cell count per clump
        cell_count_dt <- clumpdt[,.N,by=clump_ID]
        
        # bf count per clump
        bfcount_dt <- clumpdt[,.(sumcount=sum(bfcount)),by=clump_ID]
        
        # clumps table
        clump_sum_dt <- merge(cell_count_dt, bfcount_dt, by='clump_ID')
        
        # subset clumps that meet urban threshold
        big_clumps <- clump_sum_dt[N >= 1500 & sumcount >= 5000,]
        
        # ids for urban and rural
        sl <- clumps
        sl[!is.na(clumps)] <- 0
        sl[clumps %in% big_clumps$clump_ID] <- 1
        
        # save urban raster (urban=1, rural=0, unsettled=NA)
        :writeRaster(sl, filename=paste(output_filepath,'/',ISO_code,'_buildings_v1_0_urban.tif',sep=''), datatype='LOG1S')
        
        rm(clumps, clumpdt, clump_sum_dt, bfcount_dt, big_clumps, sl, cell_count_dt)
        
      }
      
      rm(bcr,bdensr,pxarea,bc)
      
    }
    
  }
  
  unlink(paste(temp_filepath))
  
  rm(centroids_wgs_master,dt,mgrid,points_IDs,urlx)
}