Estimation of the number of people affected by drought in Eastern Cape, South Africa

B5a = Hectares of crops affected * average workers per hectare

B5b = Hectares of grassland affected * average workers per hectare

provided by: UNU-EHS

Outline

Part D: Estimation of SFDRR B5

Part C: Data access and preparation

First, Part C is a guide to download the required data for this next part of the Recommended Practice. Second, this part prepares the data (e.g. resizing to Eastern Cape boundaries) to be compatible for the analysis.

NOTE

Please create a folder called Data in whose subfolders you can store all the data needed for this practice.

Step 1: Download land cover map

  • Register and login in the GIS data archive of the South African Department of Environmental Affairs (DEA) via https://egis.environment.gov.za/gis_data_downloads.
  • Download the Land Cover 2013, 72 Classes - Reprojected GEO (see image below). Please store the unzipped data in the folder “./Data/LC2014”.

Step 2: Download local municipality shapefile

  • Via http://www.demarcation.org.za/site/shapefiles/ you can download shapefiles for all local municipalities of South Africa. Download the shapefile 2016 Boundaries - Local Municipalities. In part A you already downloaded the “District Municipalities”, now the lowest spatial level is needed. Please store the unzipped data in the folder “./Data/adm/lm”.

Step 3: Download agricultural household survey

  • Download Agricultural households municipal data for Eastern Cape 2011 via www.statssa.gov.za/agricultural_household_provinces/Eastern%20Cape.xlsx. This data set includes the number of households in different agricultural sectors (e.g. livestock (LS) and crop production). The data is provided by the Statistics South Africa. Please store the data in the folder “./Data/StatSA/originals”.

Step 4: Download household size data

Step 5: Load hazard map

  • Save the hazard map for the year 2015 (2015_2016_ThreeCl.tif) which was produced in Part B in the folder “./Data/H2015_250m”. The values are as follows:
Map Value Hazard class
0 Not affected (H0)
1 Damaged (H1)
2 Destroyed (H2)

Step 6: Load province boundary shapefile

  • Store the province boundary shapefile (EC.shp) which was produced in Part A in the folder “./Data/adm/pl”.

Step 7: Prepare the R code

Download the R script for this part at the beginning of this page under Part C.

NOTE

  • To write lines in R as comments use “#” before starting the sentence.
  • Do not use backslashes (’') to enter data paths but “/”.
  • Folder names should not have empty spaces because R cannot read them.
  • Load required packages using the function library().
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.2-18, (SVN revision 718)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Karen/Documents/R/win-library/3.3/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Karen/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-7



  • Set working directory to the folder Data where all data is stored.
setwd("X:/YOUR/PATH/Data")



  • Create a new folder (outputs) in the data folder where your outputs are stored.
dir.create("outputs")



  • Set memory limit to improve the computing capacity of R.
memory.limit(50000)

Step 8: Subset administrative boundaries and land cover map

Administrative boundaries and the LC map are only available on national level. Thus, the LC map and the local municipality boundaries are clipped to the size of Eastern Cape (EC).

  • Read in and subset the local municipality shapefile to the boundaries of Eastern Cape and write the output shapefile.
SA_lm<-readOGR("./adm/lm/Local_Municipalities_2016.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\UN-SPIDER\WP2\UNU-EHS\Data\adm\lm\Local_Municipalities_2016.shp", layer: "Local_Municipalities_2016"
## with 213 features
## It has 13 fields
## Integer64 fields read as strings:  OBJECTID
lm<-subset(SA_lm, PROVINCE=="EC")
writeOGR(lm,"./outputs/lm2016.shp", "lm", driver = "ESRI Shapefile")
head(SA_lm)
  • Load the EC boundary shapefile and the national land cover map, clip the land cover map to the boundaries of Eastern Cape and write the clipped map.
EC<-readOGR("./adm/pl/EC.shp")
SA_LC2014<-raster("./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 9: Resample the hazard map

As the hazard map (250m) and the landcover map (30m) have different resolutions, the hazard map needs to be resampled to the pixel size of the landcover map. Before doing that, make sure they have the same projection (in this example lat/lon).

  • Load hazard map and use a nearest neighbour interpolation to resample the hazard map to the pixel size of the landcover map.
H2015_250m<-raster("./H2015_250m/2015_2016_ThreeCl.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 10: Extract relevant land cover

The land cover map contains several LC classes including artificial, water and forest, which are not relevant for this assessment. In this example, we specifically focus on grassland for livestock and cropland for crop production. Thus, only grass- and cropland areas need to be extracted. The class value for grassland is 7, while class values for cropland are 10-31 (see also legend 72 x class land-cover legend.pdf in Report_Metadata).

  • Add a vector with all values of the LC map.
a<-c(1:72) # add a vector with all the values of the LC map



  • Add a vector in which the 7th position is 1, the rest NA for grassland. Set 1 for the positions 10-31 for cropland respectively.
  • Create matrix (recl_grassland/ recl_crop), which contains the original and the designated values for each LC class.
  • Reclassify the land cover map (grass/crop=1; other=NA).
#grassland
# grassland =7
b<-c(rep(NA,6),1,rep(NA, 65)) # add a vector in which the 7th positio is 1, th rest NA
recl_grassland<-cbind(a,b)
LC_grassland<-reclassify(LC2014,rcl = recl_grassland)
writeRaster(LC_grassland, "./outputs/LCGrass2014.tif")

#cropland
# all cultivated agricultural land=1, non-agricultural=NA (10-31=1)
c<-c(rep(NA,9), rep(1,22), rep(NA,41))
recl_crop<-cbind(a,c)
LC_crop<-reclassify(LC2014, rcl = recl_crop)
writeRaster(LC_crop, "./outputs/LCCrop2014.tif")
head(recl_grassland, n=10)
##        a  b
##  [1,]  1 NA
##  [2,]  2 NA
##  [3,]  3 NA
##  [4,]  4 NA
##  [5,]  5 NA
##  [6,]  6 NA
##  [7,]  7  1
##  [8,]  8 NA
##  [9,]  9 NA
## [10,] 10 NA

Part D: Estimation of SFDRR B5

Part D contains the final computation of outputs. The hazard information is combined with census data in order to estimate the number of affected people by drought. For Step 1 and Step 3 to Step 5 we show the process for grassland. However, in the original script cropland is also included.

Step 1: Extract spatial statistics

In order to analyse the affected grassland (or cropland) in the respective hazard class, the hazard map is overlaid with grassland (or cropland) and the number of hectares of each hazard class per local municipality are estimated.

  • Crop and mask hazard map with grassland (or cropland). Write the output.
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)
  • Extract class count (H0-H2) for each local municipality, create a data frame and write result as csv.
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
write.csv(lm@data, file = "./outputs/H_grass2015.csv")
Hgrass2015<-lm@data



  • Multiply the number of pixels (each 900m²) per class with 0.09 to compute the number of hectares within each class and the total grassland in hectares.
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
Hgrass2015

Step 2: Prepare census data

In order to estimate the number of affected people information on the livestock and crop dependent population per municipality is needed. Thus, census data needs to be reorganized to yield in the relevant information. See table below and save file in StatSA folder as AgrDepPop.csv. Alternatively, the file can also be downloaded here ##LINK###.

Column Content description Data source
MUNICNAME Name of the municipality identical to the name of the shapefile (lm.shp). Sort alphabetically.
ID ID EC1 to EC33. Please first sort the MUNICNAME aphabetically and then give IDs
agriHH2011 Number of agricultural households by local municipality StatSA (2011)
[X]Prod Number households in production activity X (Livestock, Poultry, Vegetables, Other crop production, Fodder grazing). StatSA (2011)
SumProdHH Sum of [X] Prod
RatioProd&agriHH Because there are doublecounts in SumProdHH a ratio between SumProdHH and agriHH2011 is calculated [agriHH2011/SumProdHH].
CropDepHH, LSDepHH Number of households associated with crop and livestock production, respectively summed as follows from [X]Prod: [(VegProd+OtherProd)multiplied with RatioProd&agriHH, and (LivestockProd+FodderGraz) multiplied with RatioProd&agriHH] Poultry and Other is neglected
HHSize2016 Average household size for that municipality in 2016 (HH based distribution). StatSA (2016)
CropDepPop, LSDepPop Crop- and livestock-dependent population per municipality. = CropDepHH, LSDepHH * HHSize2016

Step 3: Calculate population density

The information of the livestock (or crop) dependent population needs to be merged with the information of the affected land. Next densities of dependent people are calculated (e.g. livestock dependent population per hectare).

  • Load the prepared census data (AgrDepPop.csv) and subset it in order to keep only the total livestock and crop dependent population.
AgrDepPop<-read.csv("./StatSA/AgrDepPop.csv") # add the census data (which has been created in Excel before!)
AgrDepPop
AgrDepPop<-subset(AgrDepPop, select=c(ID, MUNICNAME, LSDepPop, CropDepPop))



  • Merge the census data with the output from Step 1 (Part D) and subset to relevant data.
GrassDepPop = merge(Hgrass2015,AgrDepPop, by ="MUNICNAME") # merge the two sets by MUNICNAME
GrassDepPop<-subset(GrassDepPop, select= c(MUNICNAME, GrassH0area, GrassH1area,GrassH2area, ID, GrassTotalArea, LSDepPop ))



  • Divide LSDepPop (or CropDepPop) by total grassland (or cropland) area to compute population density.
GrassDepPop$LSPopDens<-GrassDepPop$LSDepPop/GrassDepPop$GrassTotalArea  # LS dependent population density

Step 4: Calculate affected population per hazard class

Then the density of livestock or crop dependent population is multiplied with the affected area, which yields in the affected population per hazard class.

  • Multiply density with affected area of 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



  • Calculate the shares by dividing the population per hazard class with the total livestock dependent population.
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 5: Calculate total affected population in crop/grassland

The number of people per hazard class (H0-H2) are then added up to yield in the total affected population dependent on livestock/crops.

  • Summarize affected population (absolute and relative) over all affected classes (H1+H2) and write output.
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") 
GrassDepPop

Step 6: Calculate total affected population

In the last step the separated information on affected livestock and crop dependent population are merged into a single .csv and also linked back to the local municipality shapefile, so that the information is also spatially available.

  • Merge the results of affected livestock and crop dependent population into a single file of total affected population (relative and absolute) by the column MUNICNAME. Write output.
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")
affPop2015
  • Merge the results to the shapefile and export the shapefile.
lm<-shapefile("./outputs/lm2016.shp")
m<-merge(lm, affPop2015, by="MUNICNAME")
writeOGR(m,"./outputs/affPop2015.shp", "affPop2015", driver = "ESRI Shapefile" )