Economic Value of direct maize production losses
Sub-indicator:– Direct Crop Loss (C2c)

Table of Contents

Part. C

Step 0 of 4 Technical Guidance Sub-indicator

Step 1 of 4 Installing required packages and creating directories

Step 2 of 4 Data preparation

Step 2.1 of 4 Resample H2015_250m to 30m

Step 2.2 of 4 Maize2015: Extract relevant LC from LC2015 (maize in LC map=5)

Step 2.3 of 4 subsetting Kiev region from the Ukraine administrative area

Part. D

Step 3 of 4 Processing

Step 3.1 of 4 Overlay the agricultural drought impact map (H2015), maize land cover for 2015 (LC2015) and Kiev map (adm2) to extract statistics

Step 3.2 of 4 Calculate expected and actual harvested ha to get deltaHa

Step 3.3 of 4 Calculate change in yield (deltaY)

Step 4 of 4 Outputs and final results

Step 4 of 4.1 Loss on Destroyed Area (C2C2) in U.S. Dollars

Step 4 of 4.2 Loss on Damaged Area (C2C1) in U.S. Dollars

Step 4 of 4.3 Maize production losses for 2015

---

Part. C

---

---

Step 0 of 4 Technical Guidance Sub-indicator

---

Economic Value of direct maize production losses

Sub-indicator: Direct Crop Loss (C2c)

Economic loss due to crops damaged and destroyed by a disaster (drought)

• C2c1 represents damged maize fields

• C2c2 represents destroyed maize fields

C2c = Loss in annual crop stocks + Loss in perennial crop stocks + Annual crop production loss + Perennial crop production loss + Crop assets loss (complete and partial)

---

Step 1 of 4 Installing required packages and creating directories

---

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

library(rgdal)
library(raster)
library(foreign)

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

setwd("C:/Users/mduguru/Desktop/Ukraine_workshop/UKR/EconLoss/Data")

create a folder to store your outputs

dir.create("outputs")

---

Step 2 of 4 Data preparation

---

---

Step 2.1 of 4 Resample H2015_250m to 30m

---

Read the land use raster into R 30m

LC2015<-raster("./LC2015/Kyiv_2015.tif")
LC2015

class       : RasterLayer 
dimensions  : 8719, 6690, 58330110  (nrow, ncol, ncell)
resolution  : 30.00059, 29.99943  (x, y)
extent      : 239770.9, 440474.8, 5452339, 5713904  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\mduguru\Desktop\Ukraine_workshop\UKR\EconLoss\Data\LC2015\Kyiv_2015.tif 
names       : Kyiv_2015 
values      : 0, 13  (min, max)

Land cover with available classes

plot(LC2015, col=colors[0:length(classes)+1], colNA=colors[1], axes=FALSE, box=FALSE)

CODE	Land Cover
0	No value
1	Artificial
2	Winter wheat
3	Winter rapeseed
4	Spring crops (wheat, barley)
5	Maize
6	Sugar beet
7	Sunflower
8	Soybeans
9	Other cereals
10	Forest
11	Grassland
12	Bare land
13	Water

Read the agricultural drought impact to be resampled from 250m to 30 m resolution of the Land cover map

H2015_250m<-raster("./H2015_250m/2015_droughtThreeCl.tif")
H2015_250m

class       : RasterLayer 
dimensions  : 1175, 1441, 1693175  (nrow, ncol, ncell)
resolution  : 0.002083333, 0.002083333  (x, y)
extent      : 29.22708, 32.22917, 49.13125, 51.57917  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\mduguru\Desktop\Ukraine_workshop\UKR\EconLoss\Data\H2015_250m\2015_droughtThreeCl.tif 
names       : X2015_droughtThreeCl 
values      : 0, 2  (min, max)

Check the projection of the agricultural impact map above in case the raster is not in UTM, the coordinate reference system needs to be changed to UTM just like the Land cover map before resampling

H2015_250m <- projectRaster (H2015_250m, crs="+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0", method = "ngb")
 H2015_250m

class       : RasterLayer 
dimensions  : 1251, 1579, 1975329  (nrow, ncol, ncell)
resolution  : 148, 232  (x, y)
extent      : 218873.7, 452565.7, 5436549, 5726781  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X2015_droughtThreeCl 
values      : 0, 2  (min, max)

Resample now

H2015<-resample(H2015_250m, LC2015,method="ngb")

H2015

class       : RasterLayer 
dimensions  : 8719, 6690, 58330110  (nrow, ncol, ncell)
resolution  : 30.00059, 29.99943  (x, y)
extent      : 239770.9, 440474.8, 5452339, 5713904  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X2015_droughtThreeCl 
values      : 0, 2  (min, max)

Write the resampled product into your output file

writeRaster(writeRaster(H2015,"./outputs/2015HUTM36_30m.tif")

---

Step 2.2 of 4 Maize2015: Extract relevant LC from LC2015 (maize in LC map=5)

---

lists the values from LC

values<-unique(values(LC2015))
values

 [1]  0 11 10 13 12  3  8  2  5  7  1  4  6

recl_m<-matrix(c(values, c(rep(NA,8),1, rep(NA,4))), ncol=2) #reclassify rasterlayer: maize=1, other=NA (see legend in data folder)
recl_m

      [,1] [,2]
 [1,]    0   NA
 [2,]   11   NA
 [3,]   10   NA
 [4,]   13   NA
 [5,]   12   NA
 [6,]    3   NA
 [7,]    8   NA
 [8,]    2   NA
 [9,]    5    1
[10,]    7   NA
[11,]    1   NA
[12,]    4   NA
[13,]    6   NA

LC_m<-reclassify(LC2015, rcl = recl_m)

A Land Cover map of Maize

plot(LC_m)

writeRaster(LC_m, "./outputs/Maize2015.tif")

---

Step 2.3 of 4 subsetting Kiev region from the Ukraine administrative area

---

adm2UKR<-shapefile("./adm2UKR/ukr_cod_admbnda_adm2_q2_2017.shp")
plot(adm2UKR)

head(adm2UKR)

adm2UKR

class       : SpatialPolygonsDataFrame 
features    : 680 
extent      : 22.13691, 40.22079, 44.38641, 52.37929  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 18
names       :  a2Pcode,             a2NameUa,            a2NameRus,      a2NameLat,            codeNameA2, codeTopoA2,              POL_ADM2,  PcodeOld2, a1Pcode,               a1NameUa,              a1NameRus,                 a1NameLat,            codeNameA1, codeTopoA1, a0Pcode, ... 
min values  : UKR01101, Жданівська, Ждановский,      Alchevska, міська рада,   81300000, міська рада, 0110100000,   UKR01, Житомирська, Житомирская, Avtonomna Respublika Krym, міська рада,   81200000,     UKR, ... 
max values  : UKR85369,   Острозька, Острожский, Zvenyhorodskyi,            район,   81400000, міська рада, 8536900000,   UKR85,         Одеська,       Одесская,               Zhytomyrska,        область,   81230100,     UKR, ... 

adm2<-subset(adm2UKR, a1NameLat=="Kyivska")
adm2<-subset(adm2, a2NameLat!="Kyivska")

Kiev Region

plot(adm2)

head(adm2)

Check the projection of the map make sure that adm2 and H2015 have the same projected coordinate system if not reproject the Kiev administrative area adm2

adm2<-spTransform(adm2, CRS("+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 "))

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

To Read in the attribute table of the shapefile for Kiev

adm2_attr<-read.dbf("./outputs/adm2.dbf") # read in the attribute table of the shapefile

adm2_attr

at<-adm2_attr[order(adm2_attr$a2NameLat),] # order the names alphabetically of the rayons

at$ID<-c(1:37) # add an ID

at<-at[,c(19,4)] # only select the names and the ID from the attribute table

write.dbf(at,"./outputs/adm2.dbf") # add attribute table back to adm2

Load the updated shapfile again

adm2<-shapefile("./outputs/adm2.shp") # load the updated shapfile again

---

Part. D

---

---

Step 3 of 4 Processing

---

---

Step 3.1 of 4 Overlay the agricultural drought impact map (H2015), maize land cover for 2015 (LC2015) and Kiev map (adm2) to extract statistics

---

Hmc<- crop(H2015,LC_m)
Hm<-mask(Hmc,LC_m)

writeRaster(Hm, "./outputs/2015Hmaize.tif", overwrite=TRUE)

Extract raster values to polygons like below

v<-extract(Hm, adm2)

Get class counts for each polygon

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

class.df

Add back to polygon data

adm2@data <- data.frame(adm2@data, class.df)

write the drought impact for maize to outputs folder

write.csv(adm2@data, file = "./outputs/H_maize.csv")

---

Step 3.2 of 4 Calculate expected and actual harvested ha to get deltaHa

---

H_maize<-read.csv("./outputs/H_maize.csv")

H_maize$X0[is.na(H_maize$X0)]<-0 # replaces all NA with 0 --> 0 ha affected
H_maize$X1[is.na(H_maize$X1)]<-0 # replaces all NA with 0 --> 0 ha affected
H_maize$X2[is.na(H_maize$X2)]<-0 # replaces all NA with 0 --> 0 ha affected

To clean data remove all data except for the ones in bracket as below

H_maize<-subset(H_maize, select = c(ID,a2NameLat, X0,X1, X2))

Calculate expected and actual harvested ha to get deltaHa

Adds column with the sum of all maize pixels X0, X1 and X2 which are in column 3-5

H_maize$pixelSum =rowSums(H_maize[, 3:5], na.rm = T)

Area in ha, 1 pixel is 900m², thus the number of pixel needs to be multiplied with 0.09 to get ha

H_maize$H0Area=H_maize$X0*0.09 
H_maize$H1Area=H_maize$X1*0.09 
H_maize$H2Area=H_maize$X2*0.09

The expected area to be harvested is the total maize crop area(H0, H1, and H2)

Sign	Drought Impact class
H0	Not Affected Maize Land
H1	Damaged Maize Land
H2	Destryed Maize Land

H_maize$expected=H_maize$H0Area+ H_maize$H1Area+ H_maize$H2Area 

The actual area is the harvested area (H0 and H1)

H_maize$actual=H_maize$H0Area+ H_maize$H1Area

Compute The difference between the expected havest hectare and the actual harvested hectare

H_maize$deltaHa= H_maize$expected - H_maize$actual

---

Step 3.3 of 4 Calculate change in yield (deltaY)

---

Note: Add the yield data (including the average yield 2010-2014)

maizeYield2015 <-read.csv("./Y2004-2015/YieldStatMaize.csv")

compute for average if the average yield was not yet calculated

maizeYield2015$avg<-((maizeYield2015$X2010 + maizeYield2015$X2011+ maizeYield2015$X2012i+ maizeYeld2015$X2013+ maizeYield2015$X2014)/5)

Remove all yield data from other years than 2015, avg and the ID which is required

maizeYield2015<-subset(maizeYield2015, select = c(ID,X2015,avg))

multiply with 0.1 for conversion from centner to tons

maizeYield2015$deltaY<-(maizeYield2015$avg - maizeYield2015$X2015)*0.1

Set all values <0 ==0 to avoid positive losses

maizeYield2015$deltaY<-replace(maizeYield2015$deltaY, maizeYield2015$deltaY<0, 0)

Add the yield data to maize drought impact

H_maizeYield = merge(maizeYield2015,H_maize, by ="ID")

look up the price from the page *** Agriculturaloutput prices*** from the price excel document UA_WP1_consolidated-database_USD the table for the year 2015 and add it as a column

H_maizeYield$P=127.7

---

Step 4 of 4 Outputs and final results

---

---

Step 4 of 4.1 Loss on Destroyed Area (C2C2) in U.S. Dollars

---

*** Please inspect the tables for clearity on the developed methodology***

H_maizeYield$C2C2= H_maizeYield$P* H_maizeYield$deltaHa* maizeYield2015$avg

---

Step 4 of 4.2 Loss on Damaged Area (C2C1) in U.S. Dollars

---

H_maizeYield$C2C1= H_maizeYield$P* H_maizeYield$deltaY* H_maizeYield$H1Area

H_maizeYield

---

Step 4 of 4.3 Maize production losses for 2015

---

Sum of damaged and destroyed maize per Rayon

H_maizeYield$maizeLoss2015= H_maizeYield$C2C1+H_maizeYield$C2C2

Total loss in maize production

overallLossMaize <-sum (H_maizeYield $ maizeLoss2015, na.rm= T)

Overall Economic Value of maize Loss in U.S.D $291493992

overallLossMaize

[1] 291493992

Store the output CSV H_maizeYield as H_LossMaize.csv in the output file

write.csv(H_maizeYield, "./outputs/H_LossMaize.csv")

Optionaly : if you want to have a shapefile with these information, the output can be merged with the adm shapefile

m<-merge(adm2, H_maizeYield, by="ID")

class       : SpatialPolygonsDataFrame 
features    : 37 
extent      : 239632.3, 440394.5, 5452195, 5713929  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 23
names       : ID,  a2NameLat.x,  X0.x,  X1.x,  X2.x, X2015,  avg, deltaY,  a2NameLat.y,  X0.y,  X1.y,  X2.y, pixelSum,  H0Area,  H1Area, ... 
min values  :  1, Baryshivskyi,     2,     1,     1,    27, 29.6,   0.00, Baryshivskyi,     0,     0,     0,        0,    0.00,    0.00, ... 
max values  : 37,  Zghurivskyi, 88683, 38765, 55152,   139, 87.6,   2.48,  Zghurivskyi, 88683, 38765, 55152,   176752, 7981.47, 3488.85, ... 

write Shapefile with all information to outputs folder

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

---

---

END