Estimating the number of people whose livelihoods were disrupted or destroyed, attributed to disasters Sub-indicators:

B5a = Hectares of crops affected * average workers per hectare

B5b = Livestock lost * average workers per hectare

Table of Contents

Part. C

Step 0 of 5 Technical Guidance Sub-indicator

Step 1 of 5 Installing required packages and creating directories

step 2 of 5 Data preparation

step 2.1 of 5 Subset of administrative boundaries and land cover map (LC map)

Step 2.2 of 5 Resample H2015_250m to 30m

Step 2.3 of 5 Extract relevant LC from LC2014 (grassland for livestock and cropland for crops)

Part. D

Step 3a of 5 Processing for Grassland

Step 3a.1 Overlay H2015, LC_grassland and lm to extract statistics

Step 3a.2 : Merge extracted statistics with census survey data to calculate population density

Step 4a of 5 hazard classes for grassland

Step 4a.1 of 5 hazard classes for grassland

Step 4a.2 of 5 Calculate total affected population for grassland damaged (H1) or destroyed (H2)

Step 3b.1 of 5 Overlay drought impact map with crop land cover(H2015LC_cropland) with local municipality (lm) to extract statistics

Step 3b.2: Merge extracted statistics with census survey data to calculate population density

Step 4b of 5 drought impact classes for cropland

Step 4b.1 of 5 Calculate affected population per drought impact class

Step 4b.2 of 5 Calculate total affected population for cropland damaged (H1) and destroyed (H2)

Step 5 of 5 calculate total affected population

---

Part. C

---

---

Step 0 of 4 Technical Guidance Sub-indicator

---

Sub-indicators:

B5a = Hecatares of crops affected * average workers per hectare

B5b = Livestock lost * average workers per hectare

---

Sep 1 of 5 Installing required packages and creating directories

---

install.packages("raster")
install.packages("rgdal")

library(rgdal)
library(raster)

Select the folder where all the data is stored in subfolders as your working directory

setwd("C:/Users/mduguru/Desktop/SouthAfricaWorkshop/SA/ADP/Data")

create a folder to store your outputs

dir.create("outputs")

---

step 2 of 4 Data preparation

---

---

Step 2.1 of 5 Subset of administrative boundaries and land cover map (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")

plot(SA_pl)

SA_pl

class       : SpatialPolygonsDataFrame 
features    : 9 
extent      : 16.45189, 32.89125, -34.83514, -22.12503  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 12

no non-missing arguments, returning NAno non-missing arguments, returning NA

names       : ID_0, ISO,       NAME_0, ID_1,       NAME_1, HASC_1, CCN_1, CCA_1,    TYPE_1, ENGTYPE_1, NL_NAME_1,         VARNAME_1 
min values  :  211, ZAF, South Africa,    1, Eastern Cape,  ZA.EC,     0,    EC, Provinsie,  Province,        NA, Eastern Transvaal 
max values  :  211, ZAF, South Africa,    9, Western Cape,  ZA.WC,     0,    WC, Provinsie,  Province,        NA,          Wes-Kaap 

head(SA_pl)

tail(SA_pl)

Eastern Cape

plot(EC)

head(EC)

Write the Eastern Cape subset shapefile product into your output file

# 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 land cover (LC map) to EC boundaries

# add LC map of the whole SA
SA_LC2014<-raster("./SA_LC2014/sa_lcov_2013-14_gti_geo_wgs84_vs22b.tif") 

# crop the raster with the EC shapefile
LC2014c<-crop(SA_LC2014,EC) 
LC2014<-mask(LC2014c, EC)
writeRaster(LC2014,"./outputs/LC2014_EC.tif")

plot(LC2014) 

---

Step 2.2 of 5 Resample H2015_250m to 30m

---

Add the hazard map which is to be resampled

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

Agricultural drought Impact map legend Values are discrete

Map Value	Class
0	Not affected
1	Damaged
2	Desroyed

---

Step 2.3 of 5 Extract relevant LC from LC2014 (grassland for livestock and cropland for crops)

---

The values of the classes can be found on the legend and metadata document in the LC map folder

values<-unique(values(LC2014)) # lists the values from LC

values

 [1]  0  7  6  5 41  9  3 65 67 66 24 25 40 23 72 32 33 68 12 11  2  1 36 35 52 10
[27] 15 14 13  4 71 70 69 17 16 18 34 63 62 50 49 61 59 48 58 57 51 42 37 38 43 60
[53] 64 46 44 45 47 55 53 54 56  8 22 39

Reclassification of Land cover for only grassland

recl_grassland<-matrix(c(values, c(NA,1,rep(NA,62))), ncol=2)
LC_grassland<-reclassify(LC2014,rcl = recl_grassland)

recl_grassland

      [,1] [,2]
 [1,]    0   NA
 [2,]    7    1
 [3,]    6   NA
 [4,]    5   NA
 [5,]   41   NA
 [6,]    9   NA
 [7,]    3   NA
 [8,]   65   NA
 [9,]   67   NA
[10,]   66   NA
[11,]   24   NA
[12,]   25   NA
[13,]   40   NA
[14,]   23   NA
[15,]   72   NA
[16,]   32   NA
[17,]   33   NA
[18,]   68   NA
[19,]   12   NA
[20,]   11   NA
[21,]    2   NA
[22,]    1   NA
[23,]   36   NA
[24,]   35   NA
[25,]   52   NA
[26,]   10   NA
[27,]   15   NA
[28,]   14   NA
[29,]   13   NA
[30,]    4   NA
[31,]   71   NA
[32,]   70   NA
[33,]   69   NA
[34,]   17   NA
[35,]   16   NA
[36,]   18   NA
[37,]   34   NA
[38,]   63   NA
[39,]   62   NA
[40,]   50   NA
[41,]   49   NA
[42,]   61   NA
[43,]   59   NA
[44,]   48   NA
[45,]   58   NA
[46,]   57   NA
[47,]   51   NA
[48,]   42   NA
[49,]   37   NA
[50,]   38   NA
[51,]   43   NA
[52,]   60   NA
[53,]   64   NA
[54,]   46   NA
[55,]   44   NA
[56,]   45   NA
[57,]   47   NA
[58,]   55   NA
[59,]   53   NA
[60,]   54   NA
[61,]   56   NA
[62,]    8   NA
[63,]   22   NA
[64,]   39   NA

Grassland land cover for Livestock

Write the Land cover containing just grassland to the output folder

writeRaster(LC_grassland, "./outputs/LCGrass2014.tif")

All cultivated agricultural land=1, non-agricultural=NA the values 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)

Cropland land cover

Write the Land cover containing just crops to the output folder

writeRaster(LC_crop, "./outputs/LCCrop2014.tif")

Part. D

---

Step 3 of 5 Processing for Grassland/ Cropland

---

---

Step 3a.1 Overlay H2015, LC_grassland and lm to extract statistics

---

Crop the hazard map to the grassland layer

# crop the hazard map with the grassland layer
HgrassC<- crop(H2015,LC_grassland)  
Hgrass<-mask(HgrassC,LC_grassland)

plot(Hgrass)

Write drought impact map file to output folder

writeRaster(Hgrass, "./outputs/H_grass2015.tif", overwrite=TRUE)

Extract raster values to polygons

v<-extract(Hgrass, lm)

Get class counts for each polygon (municipality)

v.counts <- lapply(v,table)

Create a data.frame where missing classes are NA

class.df <- as.data.frame(t(sapply(v.counts,'[',1:length(unique(Hgrass)))))

Add back to polygon data

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

Local municipality(lm)

plot(lm)

Hgrass2015

Calculating the area for each class in hectares

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

Sum of all grass hectare area H0+H1+H2 = Total grass area

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

---

Step 3a.2 : Merge extracted statistics with census survey data to calculate population density

---

Add agriculture dependent population

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))

AgrDepPop

** Merge the two sets by MUNICNAME**

GrassDepPop = merge(Hgrass2015,AgrDepPop, by ="MUNICNAME")  
GrassDepPop<-subset(GrassDepPop, select= c(MUNICNAME, GrassH0area, GrassH1area,GrassH2area, ID, GrassTotalArea, LSDepPop ))

LS dependent population density

GrassDepPop$LSPopDens<-GrassDepPop$LSDepPop/GrassDepPop$GrassTotalArea
GrassDepPop

---

Step 4a of 5 drought impact classes for grassland

---

---

Step 4a.1 of 5 Calculate affected population per drought impact class

---

# LS Dependent population in H0 
GrassDepPop$LSPopH0<-GrassDepPop$LSPopDens*GrassDepPop$GrassH0area 
# LS Dependent population in H1
GrassDepPop$LSPopH1<-GrassDepPop$LSPopDens*GrassDepPop$GrassH1area 
# LS Dependent population in H2
GrassDepPop$LSPopH2<-GrassDepPop$LSPopDens*GrassDepPop$GrassH2area 

Calculate the shares of pop in drought impact class/ total population for grassland

GrassDepPop$ShareLSPopH0<-GrassDepPop$LSPopH0/GrassDepPop$LSDepPop
GrassDepPop$ShareLSPopH1<-GrassDepPop$LSPopH1/GrassDepPop$LSDepPop
GrassDepPop$ShareLSPopH2<-GrassDepPop$LSPopH2/GrassDepPop$LSDepPop
GrassDepPop

---

Step 4a.2 of 5 Calculate total affected population for grassland damaged (H1) and destroyed (H2)

---

Affected population (H1+H2)

GrassDepPop$LSPopAff<-GrassDepPop$LSPopH1+GrassDepPop$LSPopH2

Calculate share of total affected pop (h1+h2)

GrassDepPop$ShareLSPopAff<-GrassDepPop$LSPopAff/GrassDepPop$LSDepPop
####### or ######
GrassDepPop$ShareLSPopAff<-GrassDepPop$ShareLSPopH1+GrassDepPop$ShareLSPopH2
GrassDepPop

Write the csv for the affected livestock dependent population (affLSPop2015.csv) to output file

write.csv(GrassDepPop,"./outputs/affLSPop2015.csv")

---

---

---

Step 3b.1. Overlay drought impact map with crop land cover(H2015LC_cropland) with local municipality (lm) to extract statistics

---

Add local municipality boundaries again, so that results from grassland analysis are removed

lm<-shapefile("./outputs/lm2016.shp")

Eastern Cape local municipalities

class       : SpatialPolygonsDataFrame 
features    : 33 
extent      : 22.73574, 30.19472, -34.21386, -30.0018  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 13
names       : OBJECTID, PROVINCE, CATEGORY,                      CAT2, CAT_B,     MUNICNAME,              NAMECODE,                        MAP_TITLE, DISTRICT, DISTRICT_N, Shape_STAr, Shape_STLe,                                       F2 
min values  :        1,       EC,        A,        Local Municipality,   BUF,     Amahlathi,     Amahlathi (EC124),     Amahlathi Local Municipality,      BUF, Alfred Nzo,  0.1226250,   2.393968,     Amahlathi Local Municipality (EC124) 
max values  :        9,       EC,        B, Metropolitan Municipality,   NMA, Walter Sisulu, Walter Sisulu (EC145), Walter Sisulu Local Municipality,      NMA,  O.R.Tambo,  2.7599461,  11.466761, Walter Sisulu Local Municipality (EC145) 

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

Extract pixels by class per municipality

v<-extract(Hcrop, lm) #Extract raster values to polygons 
v.counts <- lapply(v,table) # Get class counts for each polygon

Create a data.frame where missing classes are NA

class.df <- as.data.frame(t(sapply(v.counts,'[',1:length(unique(Hcrop))))) 

Add pixel count back to polygon data per municipality

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") 

Subset pixel counts for classes and compute the hectares in each class by multiplying by 0.09

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

multiply to get ha

Hcrop2015$CropH0area<-Hcrop2015$X0*0.09
Hcrop2015$CropH1area<-Hcrop2015$X1*0.09
Hcrop2015$CropH2area<-Hcrop2015$X2*0.09
Hcrop2015

Total crop area estimation

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

---

Step 3b.3: Merge extracted statistics with census survey data to calculate population density

---

Merge the two sets by MUNICNAME

CropDepPop = merge(Hcrop2015,AgrDepPop, by ="MUNICNAME")
CropDepPop

Subset crop required columns

CropDepPop<-subset(CropDepPop, select= c(MUNICNAME, CropH0area, CropH1area,CropH2area, ID, CropTotalArea, CropDepPop ))
CropDepPop

Crop dependent population density

CropDepPop$CropPopDens<-CropDepPop$CropDepPop/CropDepPop$CropTotalArea
CropDepPop

---

Step 4b of 5 drought impact classes for cropland

---

---

Step 4b.1 of 5 Calculate affected population per drought impact class

---

# crop  Dependent population in H0
CropDepPop$CropPopH0<-CropDepPop$CropPopDens*CropDepPop$CropH0area
# crop  Dependent population in H1
CropDepPop$CropPopH1<-CropDepPop$CropPopDens*CropDepPop$CropH1area
# crop  Dependent population in H2
CropDepPop$CropPopH2<-CropDepPop$CropPopDens*CropDepPop$CropH2area 

Calculate the shares of pop in drought impact class/ total population for cropland

CropDepPop$ShareCropPopH0<-CropDepPop$CropPopH0/CropDepPop$CropDepPop 
CropDepPop$ShareCropPopH1<-CropDepPop$CropPopH1/CropDepPop$CropDepPop
CropDepPop$ShareCropPopH2<-CropDepPop$CropPopH2/CropDepPop$CropDepPop
CropDepPop

---

Step 4b.2 of 5 Calculate total affected population for cropland damaged (H1) and destroyed (H2)

---

Affected population (H1+H2)

CropDepPop$CropPopAff<-CropDepPop$CropPopH1+CropDepPop$CropPopH2

Calculate share of total affected pop (h1+h2)

CropDepPop$ShareCropPopAff<-CropDepPop$CropPopAff/CropDepPop$CropDepPop
####### or ######
CropDepPop$ShareCropPopAff<-CropDepPop$ShareCropPopH1+CropDepPop$ShareCropPopH2
CropDepPop

Write the csv for the affected crop dependent population (affCropPop2015.csv) to output file

write.csv(CropDepPop,"./outputs/affCropPop2015.csv")

---

Step 5 of 5 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)
affPop2015

write.csv(affPop2015,"./outputs/affPop2015.csv")

*** Additionaly: If you want to have a shapefile with these information, the output can be merged with the local municipality (lm) shapefile ***

lm<-shapefile("./outputs/lm2016.shp")
m<-merge(lm, affPop2015, by="MUNICNAME")

writeOGR(m,"./outputs/affPop2015.shp", "affPop2015", driver = "ESRI Shapefile" )