forked from SFrav/Ruminant-feed-balance-2024
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathextractResults.R
96 lines (83 loc) · 4.02 KB
/
extractResults.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
# Extract results at state
library(terra)
#library(stars)
library(sf)
library(exactextractr)
library(rnaturalearth)
library(dplyr)
library(tidyr)
library(stars)
# avoid scientific notation
options(scipen = 999)
# root folder
root <- "."
country <- "Nigeria"
# paths
spatialDir <- paste0(root, "/src/3Balance-estimates/", country, "/SpatialData")
Results_dir <- paste0(root, "/src/3Balance-estimates/", country, "/Results")
aoi2 <-read_sf(paste0(spatialDir, "/inputs/aoi2.shp"))
# extract DM
feed_list <- c("crop", "grass", "browse")
tFeed <- rast(paste0(spatialDir, "/outputs/Feed_", feed_list, "DM2023.tif"))
feedDM <- exact_extract(tFeed, aoi2, "sum")
feedDM_stats <- cbind(LGA = aoi2$NAME_2, feedDM)
colnames(feedDM_stats)[2:4] <- c("cropDM_kg", "grassDM_kg", "browseDM_kg")
write.csv(feedDM_stats, paste0(Results_dir, "/FeedDM_LGA.csv"), row.names = FALSE)
# extract diet proportions
feed_list <- c("crop", "grass", "browse", "after")
tFeed <- rast(paste0(spatialDir, "/outputs/Feed_", feed_list, "DM2023.tif"))
feedDM <- exact_extract(tFeed, aoi2, "sum")
feedDM_stats <- cbind(LGA = aoi2$NAME_2, feedDM)
colnames(feedDM_stats)[2:4] <- c("cropDM_kg", "grassDM_kg", "browseDM_kg")
write.csv(feedDM_stats, paste0(Results_dir, "/FeedDM_LGA.csv"), row.names = FALSE)
# extract cropping and non-cropping days
season_list <- c("cropping", "dry")
tSeasons <- rast(paste0(spatialDir, "/outputs/", season_list, "Days_2023.tif"))
seasonLen <- exact_extract(tSeasons, aoi2, "mean")
seasonLen_stats <- cbind(LGA = aoi2$NAME_2, seasonLen)
colnames(seasonLen_stats)[2:3] <- c("croppingDays", "dryDays")
seasonLen_stats <- seasonLen_stats %>% mutate_if(is.numeric, round, digits = 0)
write.csv(seasonLen_stats, paste0(Results_dir, "/SeasonLength.csv"), row.names = FALSE)
# extract livestock population
lvCategories <- c("Bull", "Cow", "Steer", "Heifer", "Calf", "Sheep", "Lamb", "Goat", "Kid")
tLv <- rast(paste0(spatialDir, "/inputs/GLW4/Dissag_GLW4_2020_",lvCategories,".tif"))
names(tLv) <- c("Bull", "Cow", "Steer", "Heifer", "Calf", "Sheep", "Lamb", "Goat", "Kid")
tSumLv <- exact_extract(tLv, aoi2, "sum")
tSumLv <- cbind(LGA = aoi2$NAME_2, tSumLv)
colnames(tSumLv)[2:10] <- c("Bulls", "Cows", "Steers", "Heifers", "Calves", "Sheep", "Lambs", "Goats", "Kids")
tSumLv <- tSumLv %>% mutate_if(is.numeric, round, digits = 0)
# tSumLv$Cattle <- tSumLv$Bulls+tSumLv$Cows+tSumLv$Steers+tSumLv$Heifers+tSumLv$Calves
# sum(tSumLv$Cattle)
write.csv(tSumLv, paste0(Results_dir, "/LivestockPopulation.csv"), row.names = FALSE)
# extract land use types and proportions
luClasses <- c("grass", "shrub", "tree", "crops")
tLU <- rast(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/Nigeria/Feed_DrySeason/LandUse/LU",luClasses,"300.tif"))
#tLUcrop <- rast(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/Nigeria/Feed_DrySeason/LandUse/LUcrops300DEA.tif"))
#tLU <- c(tLU, tLUcrop)
tMeanLU <- exact_extract(tLU, aoi2, "mean")
tMeanLU <- cbind(LGA = aoi2$NAME_2, tMeanLU)
colnames(tMeanLU)[2:5] <- c("grass", "shrub", "tree", "crop")
tMeanLU <- tMeanLU %>% mutate(across(where(is.numeric), ~ . * 100))
write.csv(tMeanLU, paste0(Results_dir, "/LandUseProportion.csv"), row.names = FALSE)
# extract burning incidences
luClasses <- c("Crop", "Grass")
tBurn <- rast(paste0(root, "/src/3Balance-estimates/Nigeria/SpatialData/inputs/Burned/2023/burn",luClasses,"Months.tif"))
tBurn <- exact_extract(tBurn, aoi2, "frac")
tBurnEvents <- cbind(LGA = aoi2$NAME_2, tBurn)
tBurnEvents <- tBurnEvents %>%
mutate(burnEventsCrop = case_when(
frac_0.burnCropMonths >= 1 ~ 0,
frac_1.burnCropMonths > 0 |
frac_2.burnCropMonths > 0 |
frac_3.burnCropMonths > 0 |
frac_4.burnCropMonths > 0 ~ 1,
TRUE ~ NA_real_),
burnEventsGrass = case_when(
frac_0.burnGrassMonths >= 1 ~ 0,
frac_1.burnGrassMonths > 0 |
frac_2.burnGrassMonths > 0 |
frac_3.burnGrassMonths > 0 |
frac_4.burnGrassMonths > 0 ~ 1,
TRUE ~ NA_real_)) %>%
dplyr::select(LGA, burnEventsCrop, burnEventsGrass)
write.csv(tBurnEvents, paste0(Results_dir, "/BurnEventsBoolean.csv"), row.names = FALSE)