######Recommended Practice: Number of people whose livelihoods were disrupted or destroyed, attributed to drought in EC
# Agricultural Dependent population EC

# indicator: B5 --> Number of people whose livelihoods were disrupted or destroyed, attributed to disaster
# required data: C-2Ca Number of hectares of crops damaged or destroyed by disasters.(to be used to establish the statistic of Number of Workers affected)
# 

# sub-indicators: 
# B5a: hecatares of crops affected * average workers per hectare
# B5b: Livestock lost * average workers per hectare

# UNISDR technical guidance p. 20 - 22


### Step1: Getting Started an data acquisition
library(raster)
library(rgdal)

# select the folder where all the data is stored in subfolders as your working directory
setwd("X:/Data") 

# create a folder to store your outputs
dir.create("outputs")

### Step 2: Data preparation
#   Step 2.1: Subset of administrative boundaries and LC map

# subset province boundaries to EC
SA_pl<-readOGR("./adm/pl/ZAF_adm1.shp") # select the level of provinces which is adm1
EC<-subset(SA_pl, NAME_1=="Eastern Cape")

# subset local municipalities to EC
SA_lm<-readOGR("./adm/lm/Local_Municipalities_2016.shp")
lm<-subset(SA_lm, PROVINCE=="EC")

# regarding coordinate reference system: the data here is all in "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,"0" 
# thus there is no need to reproject the data 
writeOGR(lm,"./outputs/lm2016.shp", "lm", driver = "ESRI Shapefile")

# clip LC map to EC boundaries
SA_LC2014<-raster("./SA_LC2014/sa_lcov_2013-14_gti_geo_wgs84_vs22b.tif") # add LC map of the whole SA
LC2014c<-crop(SA_LC2014,EC) # crop the raster with the EC shapefile
LC2014<-mask(LC2014c, EC)
writeRaster(LC2014,"./outputs/LC2014_EC.tif")


#   Step 2.2: Resample H2015_250m to 30m
H2015_250m<-raster("./H2015_250m/2015_2016_droughtThreecl.tif")  # add the hazard map which is to be resampled
H2015<-resample(H2015_250m, LC2014,method="ngb") # resample with nearest neighbour method to not loose category values
writeRaster(H2015,"./outputs/2015H_30m.tif")

#   Step 2.3: Extract relevant LC from LC2014 (grassland for livestock and cropland for crops)
# the values of the classes can be found in the metadata document in the LC map folder 
values<-unique(values(LC2014)) # lists the values from LC
# livestock (grasslands (7))
values # shows the values of th raster
recl_grassland<-matrix(c(values, c(NA,1,rep(NA,62))), ncol=2)
LC_grassland<-reclassify(LC2014,rcl = recl_grassland)
writeRaster(LC_grassland, "./outputs/LCGrass2014.tif")

# all cultivated agricultural land=1, non-agricultural=NA (10-31=1)
recl_crop<-matrix(c(values, c(rep(NA,10),1,1,NA,1,rep(NA,6),1,1,rep(NA,3),rep(1,4),rep(NA,4),rep(1,3), rep(NA,26),1, NA)), ncol=2)
LC_crop<-reclassify(LC2014, rcl = recl_crop)
writeRaster(LC_crop, "./outputs/LCCrop2014.tif")


# Step 3: Processing
  #   Step 3.1. Overlay H2015, LC_grassland/LC_cropland and lm and extract statistics
    ####### grass

HgrassC<- crop(H2015,LC_grassland) # crop the hazard map with the grassland layer 
Hgrass<-mask(HgrassC,LC_grassland)
writeRaster(Hgrass, "./outputs/H_grass2015.tif", overwrite=TRUE)

v<-extract(Hgrass, lm) #Extract raster values to polygons 
v.counts <- lapply(v,table) # Get class counts for each polygon
class.df <- as.data.frame(t(sapply(v.counts,'[',1:length(unique(Hgrass)))))  # Create a data.frame where missing classes are NA

lm@data <- data.frame(lm@data, class.df)# Add back to polygon data
head(lm@data)
write.csv(lm@data, file = "./outputs/H_grass2015.csv")
Hgrass2015<-lm@data

Hgrass2015<-subset(Hgrass2015, select = c(OBJECTID,MUNICNAME,X0, X1, X2))

Hgrass2015$GrassH0area<-Hgrass2015$X0*0.09 # multiply to get ha
Hgrass2015$GrassH1area<-Hgrass2015$X1*0.09
Hgrass2015$GrassH2area<-Hgrass2015$X2*0.09

Hgrass2015$GrassTotalArea<-Hgrass2015$GrassH0area+ Hgrass2015$GrassH1area+Hgrass2015$GrassH2area

  # Step 3.3: Merge extracted statistics with census survey data to calculate population density
AgrDepPop<-read.csv("./StatSA/AgrDepPop.csv") # add agr dependent population (which has been created in Excel before!)
AgrDepPop<-subset(AgrDepPop, select=c(ID, MUNICNAME, LSDepPop, CropDepPop))
GrassDepPop = merge(Hgrass2015,AgrDepPop, by ="MUNICNAME") # merge the two sets by MUNICNAME
GrassDepPop<-subset(GrassDepPop, select= c(MUNICNAME, GrassH0area, GrassH1area,GrassH2area, ID, GrassTotalArea, LSDepPop ))
GrassDepPop$LSPopDens<-GrassDepPop$LSDepPop/GrassDepPop$GrassTotalArea  # LS dependent population density  

# Step 4: Outputs
  #Step4 4.1.: Calculate affected population per hazard class 
GrassDepPop$LSPopH0<-GrassDepPop$LSPopDens*GrassDepPop$GrassH0area # LS Dependent population in H0 
GrassDepPop$LSPopH1<-GrassDepPop$LSPopDens*GrassDepPop$GrassH1area # LS Dependent population in H1
GrassDepPop$LSPopH2<-GrassDepPop$LSPopDens*GrassDepPop$GrassH2area # LS Dependent population in H2
GrassDepPop$ShareLSPopH0<-GrassDepPop$LSPopH0/GrassDepPop$LSDepPop # calculate the shares of pop in hazard class/ total population
GrassDepPop$ShareLSPopH1<-GrassDepPop$LSPopH1/GrassDepPop$LSDepPop
GrassDepPop$ShareLSPopH2<-GrassDepPop$LSPopH2/GrassDepPop$LSDepPop


  #Step 4.2.: Calculate total affected population of affected crop/grassland
GrassDepPop$LSPopAff<-GrassDepPop$LSPopH1+GrassDepPop$LSPopH2 # affected population (H1+H2)
GrassDepPop$ShareLSPopAff<-GrassDepPop$LSPopAff/GrassDepPop$LSDepPop # calculate share of total affected pop (h1+h2)
write.csv(GrassDepPop,"./outputs/affLSPop2015.csv")


# Step 3: Processing
#   Step 3.1. Overlay H2015, LC_grassland/LC_cropland and lm and extract statistics
     #######crop

# add local municipality boundaries again, so that results from grassland analysis are removed
lm<-shapefile("./outputs/lm2016.shp")

HcropC<- crop(H2015,LC_crop)
Hcrop<-mask(HcropC,LC_crop)
writeRaster(Hcrop, "./outputs/H_crop2015.tif", overwrite=TRUE)

v<-extract(Hcrop, lm) #Extract raster values to polygons 
v.counts <- lapply(v,table) # Get class counts for each polygon
class.df <- as.data.frame(t(sapply(v.counts,'[',1:length(unique(Hcrop)))))  # Create a data.frame where missing classes are NA

lm@data <- data.frame(lm@data, class.df)# Add back to polygon data
head(lm@data)
write.csv(lm@data, file = "./outputs/H_crop2015.csv")
Hcrop2015<-lm@data
Hcrop2015<-read.csv("./outputs/H_crop2015.csv") 

Hcrop2015<-subset(Hcrop2015, select = c(OBJECTID,MUNICNAME,X0, X1, X2))

Hcrop2015$CropH0area<-Hcrop2015$X0*0.09 # multiply to get ha
Hcrop2015$CropH1area<-Hcrop2015$X1*0.09
Hcrop2015$CropH2area<-Hcrop2015$X2*0.09

Hcrop2015$CropTotalArea<-Hcrop2015$CropH0area+Hcrop2015$CropH1area+Hcrop2015$CropH2area

  # Step 3.3: Merge extracted statistics with census survey data to calculate population density
CropDepPop = merge(Hcrop2015,AgrDepPop, by ="MUNICNAME") # merge the two sets by MUNICNAME
CropDepPop<-subset(CropDepPop, select= c(MUNICNAME, CropH0area, CropH1area,CropH2area, ID, CropTotalArea, CropDepPop ))
CropDepPop$CropPopDens<-CropDepPop$CropDepPop/CropDepPop$CropTotalArea  # crop dependent population density  

###Step 4: Outputs
  #Step4.1.: Calculate affected population per hazard class 
CropDepPop$CropPopH0<-CropDepPop$CropPopDens*CropDepPop$CropH0area # crop  Dependent population in H0
CropDepPop$CropPopH1<-CropDepPop$CropPopDens*CropDepPop$CropH1area # crop  Dependent population in H1
CropDepPop$CropPopH2<-CropDepPop$CropPopDens*CropDepPop$CropH2area # crop  Dependent population in H2
CropDepPop$ShareCropPopH0<-CropDepPop$CropPopH0/CropDepPop$CropDepPop 
CropDepPop$ShareCropPopH1<-CropDepPop$CropPopH1/CropDepPop$CropDepPop
CropDepPop$ShareCropPopH2<-CropDepPop$CropPopH2/CropDepPop$CropDepPop

  #Step4.2.: Calculate total affected population of affected crop/grassland

CropDepPop$CropPopAff<-CropDepPop$CropPopH1+CropDepPop$CropPopH2 # affected population (H1+H2)
CropDepPop$ShareCropPopAff<-CropDepPop$CropPopAff/CropDepPop$CropDepPop # calculate share of total affected pop (h1+h2)
write.csv(CropDepPop,"./outputs/affCropPop2015.csv")

  # Step4.3.: calculate total affected population

# merge the two datasets
affPop2015<-merge(CropDepPop,GrassDepPop, by="MUNICNAME")
affPop2015$affPop<-affPop2015$CropPopAff+affPop2015$LSPopAff
affPop2015$ShareaffPop<-affPop2015$affPop/(affPop2015$CropDepPop+affPop2015$LSDepPop)
write.csv(affPop2015,"./outputs/affPop2015.csv")

## additional: if you want to have a shapefile with these information, the output can be merged with the lm shapefile
lm<-shapefile("./outputs/lm2016.shp")
m<-merge(lm, affPop2015, by="MUNICNAME")
writeOGR(m,"./outputs/affPop2015.shp", "affPop2015", driver = "ESRI Shapefile" )