-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathget_CY_FRRecom.R
136 lines (96 loc) · 5.71 KB
/
get_CY_FRRecom.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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
####################################################################################################
####################################################################################################
## 1. get CY for the 5 FCY: Note both WLY and CY are in dry matter
#####################################################################################################
require(lpSolve)
source("QUEFTS_function.R")
## WLY
LINTUL_WLY <- readRDS("LINTUL_WLY_test.RDS")
LINTUL_WLY <- LINTUL_WLY[order(LINTUL_WLY$plantingDate, LINTUL_WLY$weekNr), ]
## INS NPK for FCY levels (is the output of random forest procedure)
SoilData_FCY1 <- readRDS("SoilData_FCY1_test.RDS")
SoilData_FCY2 <- readRDS("SoilData_FCY2_test.RDS")
SoilData_FCY3 <- readRDS("SoilData_FCY3_test.RDS")
SoilData_FCY4 <- readRDS("SoilData_FCY4_test.RDS")
SoilData_FCY5 <- readRDS("SoilData_FCY5_test.RDS")
## get CY for every FCY
CY_FCY1 <- get_CY(SoilData = SoilData_FCY1, LINTUL_WLY=LINTUL_WLY, crop = "cassava", country = "Bu")
CY_FCY2 <- get_CY(SoilData = SoilData_FCY2, LINTUL_WLY=LINTUL_WLY, crop = "cassava", country = "Bu")
CY_FCY3 <- get_CY(SoilData = SoilData_FCY3, LINTUL_WLY=LINTUL_WLY, crop = "cassava", country = "Bu")
CY_FCY4 <- get_CY(SoilData = SoilData_FCY4, LINTUL_WLY=LINTUL_WLY, crop = "cassava", country = "Bu")
CY_FCY5 <- get_CY(SoilData = SoilData_FCY5, LINTUL_WLY=LINTUL_WLY, crop = "cassava", country = "Bu")
####################################################################################################
####################################################################################################
## 2. generate fertilizer rates for the 5 FCYs: This is done on my PC because it requires multiple cores
#####################################################################################################
rm(list=ls())
require(lpSolve)
library(parallel)
require(tidyr)
require(gtools)
require(plyr)
source("QUEFTS_Functions.R")
maxInv <- 403600 ## 200 USD
rootUP <- 700000 ##347 USD/ton = 700 000 RWF /ton of cassava
areaHa <- 1
country <- "Bu"
riskAtt <- 2
crop = "cassava"
fertilizers <- data.frame(type=c("FOMI_IMBURA", "FOMI_BAGARA", "FOMI_TOTAHAZA"),
N_cont= c(0.09, 0.11, 0.21), P_cont=c(0.22*0.44, 0, 0), K_cont=c(0.04*0.83, 0.22*0.83, 0.08*0.83),
costPerBag = c(85390, 77600, 77000), bagWeight = c(50,50,5), price=c(85390/50, 77600/50, 77000/50))
### this need to be done for every FCY; the following shows how to do it for one FCY, FCY2
FCY <- "FCY2"
SoilData_FCY <- SoilData_FCY2
SoilData_FCY$location <- paste(SoilData_FCY$lat, SoilData_FCY$lon, sep="_")
### For the paper based and chabot, the recommenations are provided at LGA/region level, not for every coordinate point within the region as defined by the regula grid
### Here under, the soil and weather data is aggregated at LGA level by taking the average of vlaues. Given the recommendation is generated by lon&lat pair regardless of
### the resolution, the lon & lat are also averaged to represent an LGA. In this way , the same script can be used for any desired level of resolution.
## aggregating the data at LGA level.
CY_FCY <- unique(merge(CY_FCY2, SoilData_FCY[, c("location", "NAME_1", "NAME_2")], by="location"))
SoilData_FCY <- unique(SoilData_FCY[, c("NAME_1","NAME_2", "soilN", "soilP", "soilK")])
SoilD <- droplevels(ddply(SoilData_FCY, .( NAME_1,NAME_2), summarise,
soilN = mean(soilN, na.rm = TRUE),
soilP = mean(soilP, na.rm = TRUE), soilK = mean(soilK, na.rm = TRUE)))
length(unique(SoilData_FCY$NAME_2))
head(SoilD)
getmonths <- function(ds){
ds$plm <- as.factor(ds$weekNr)
levels(ds$plm)[levels(ds$plm) %in% 1:4] <- "January"
levels(ds$plm)[levels(ds$plm) %in% 5:8] <- "February"
levels(ds$plm)[levels(ds$plm) %in% 9:13] <- "March"
levels(ds$plm)[levels(ds$plm) %in% 14:17] <- "April"
levels(ds$plm)[levels(ds$plm) %in% 18:22] <- "May"
levels(ds$plm)[levels(ds$plm) %in% 23:26] <- "June"
levels(ds$plm)[levels(ds$plm) %in% 27:30] <- "July"
levels(ds$plm)[levels(ds$plm) %in% 31:35] <- "August"
levels(ds$plm)[levels(ds$plm) %in% 36:39] <- "September"
levels(ds$plm)[levels(ds$plm) %in% 40:43] <- "October"
levels(ds$plm)[levels(ds$plm) %in% 44:48] <- "November"
levels(ds$plm)[levels(ds$plm) %in% 49:53] <- "December"
}
ds <-getmonths(ds=CY_FCY)
dsMonthly <- ddply(ds, .(NAME_1, NAME_2, zone, daysOnField, plm), summarise,
water_limited_yield = mean(water_limited_yield), CurrentYield = mean(CurrentYield),
lat = mean(lat), long = mean(long), pl_Date = mean(pl_Date), weekNr = mean(weekNr))
dsMonthly <- dsMonthly[!dsMonthly$plm %in% c("February", "March","April","May","June", "July"), ]
head(dsMonthly)
dsMonthly$pl_Date <- round(dsMonthly$pl_Date, digits=0)
dsMonthly$weekNr <- round(dsMonthly$weekNr, digits=0)
wlyfcy <- dsMonthly
wlyfcy$location <- paste(wlyfcy$lat, wlyfcy$long, sep="_")
wlyfcy <- wlyfcy[, colnames(CY_FCY)]
SoilD <- unique(merge(SoilD, unique(wlyfcy[, c("NAME_2", "location", "lat", "long")]), by="NAME_2"))
SoilD$lon <- SoilD$long
latlon <- unique(SoilD$location)
## The parallel computation generated the FR advice for every lon-lat using the soil and weather data provided.
## output is saved in "/home/akilimo/projects/AKILIMO_Modeling/QUEFTS/Output/" by country and FCY
parallelCluster <- parallel::makeCluster(parallel::detectCores(16))
print(parallelCluster)
for(ll in latlon){
print(ll)
SoilData_ll <- SoilD[SoilD$location == ll, ]
wlyll <- wlyfcy[wlyfcy$location == ll, ]
FR_parallel_byLoc(WLY_CY_FCY_ll = wlyll, SoilData_ll = SoilData_ll, maxInv = maxInv, fertilizers=fertilizers, rootUP = rootUP,
areaHa=1, country=country, FCY = FCY, riskAtt = riskAtt, crop = crop)
}