## note ########################################################################################### ### Data source and run LINTUL ########################################################################################### ## define the country and the year LINTUL is going to run Country_year <- rbind(data.frame(country="Rwanda", Q3Year = 2010, Q2Year = 2019, Q1Year = 2016), data.frame(country="Burundi", Q3Year = 2015, Q2Year = 2011, Q1Year = 2016), data.frame(country="DRC", Q3Year = 2012, Q2Year = 2019, Q1Year = 2016) ) ## define planting and harvest dates as the count of number of days in the year. ## for cassava we assume it can have weekly planting throughout the year and ## can stay in he field up to 15 months = 455 growing days Planting_Harvest_days <- data.frame(startDate=seq(1, 365, 7), endDate = (seq(1, 365, 7) + 454), weekNr = seq(1:53), months = c(rep("Jan", 5), rep("Feb", 4), rep("Mar", 5), rep("Apr", 4), rep("May", 4), rep("Jun", 5), rep("Jul", 4), rep("Aug", 5), rep("Sep", 4), rep("Oct",4), rep("Nov",5), rep("Dec", 4))) ## given in Rwanda planting happens in Sep, Oct, Feb and Mar we need to filter Planting_Harvest_RW <- Planting_Harvest_days[Planting_Harvest_days$months %in% c("Sep", "Oct", "Feb", "Mar"),] Planting_Harvest_BU <- Planting_Harvest_days[!Planting_Harvest_days$months %in% c("May", "Jun", "Jul", "Aug"),] Planting_Harvest_DRC <- Planting_Harvest_BU ### fertilizers used NPK 171717 (31USD/ 50 kg) ### Cassava gate price/ ton: Farmer-farmer (69USD/ton) and Farmer- factory (89USD/ton) ### given the sever does not have multiple cores to parallel the work, ### the following is done on my pc setwd("/home/akilimo/projects/SoilNPK/lintul") source('lINTUL2_CASSAVA.R') runLINTUL_server(year=2016, country = "Rwanda", Planting_Harvest_days = Planting_Harvest_RW) runLINTUL_server(year=2019, country = "Rwanda", Planting_Harvest_days = Planting_Harvest_RW) runLINTUL_server(year=2010, country = "Rwanda", Planting_Harvest_days = Planting_Harvest_RW) runLINTUL_server(year=2015, country = "Burundi", Planting_Harvest_days = Planting_Harvest_RW) runLINTUL_server(year=2011, country = "Rwanda", Planting_Harvest_days = Planting_Harvest_RW) runLINTUL_server(year=2016, country = "Rwanda", Planting_Harvest_days = Planting_Harvest_RW) ## checking output require(ggplot2) require(tidyr) setwd("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2019") wlyFiles <- list.files(path=getwd()) wly1 <- readRDS(wlyFiles[100]) summary(wly1) wly1 <- wly1[, c("plantingDate", "WLY_242","WLY_270", "WLY_305", "WLY_333", "WLY_361", "WLY_396", "WLY_424", "WLY_452")] names(wly1) <- c("plantingDate","M8", "M9","M10", "M11", "M12", "M13", "M14", "M15" ) wly1 <- gather(wly1, Hdate, WLY, M8:M15) wly1$Hdate <- factor(wly1$Hdate, levels = c("M8", "M9","M10", "M11", "M12", "M13", "M14", "M15")) head(wly1) ggplot(wly1, aes(plantingDate, WLY, col=Hdate))+ geom_point()+ facet_grid(.~Hdate) ############################################################################### ##################### select max of the 3 years################################ require(gtools) require(plyr) setwd("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2010") latlongQ3 <- mixedsort(list.files(path = getwd(), pattern = ".RDS")) latlongQ3 <- substr(latlongQ3, 5, 30) latlongQ3 <- unlist(strsplit(latlongQ3, ".RDS")) str(latlongQ3) setwd("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2016") latlongQ1 <- mixedsort(list.files(path = getwd(), pattern = ".RDS")) latlongQ1 <- substr(latlongQ1, 5, 30) latlongQ1 <- unlist(strsplit(latlongQ1, ".RDS")) str(latlongQ1) setwd("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2019") latlongQ2 <- mixedsort(list.files(path = getwd(), pattern = ".RDS")) latlongQ2 <- substr(latlongQ2, 5, 30) latlongQ2 <- unlist(strsplit(latlongQ2, ".RDS")) str(latlongQ2) all_locs <- unique(c(latlongQ3, latlongQ1, latlongQ2)) str(all_locs) #' For every latlong, this fucntion takes the max of the 3 years WLY, for every planting day and corresponding harvesing days between 8 and 15 month #' @param ll #' @return #' #' @author Meklit #' @export WLY_3years <- function(ll, Planting_Harvest, latlongQ3, latlongQ1, latlongQ2){ # Planting_Harvest <- data.frame(st=seq(1, 365, 7), en = (seq(1, 365, 7) + 454), PlweekNr = seq(1:53)) # WLY_Year1 <- NULL WLY_Year2 <- NULL WLY_Year3 <- NULL WLY_Year1_hd <- NULL WLY_Year2_hd <- NULL WLY_Year3_hd <- NULL harvest8_15 <- paste("WLY_", seq(235, 455, 7), sep="") WLY_8to15 <- NULL for(hdays in harvest8_15){ if(ll %in% latlongQ3){ WLY_Year1 <- readRDS(file=paste("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2010/WLY_", ll, ".RDS", sep="")) WLY_Year1_hd <- WLY_Year1[, c("lat", "long", "plantingDate", hdays)] colnames(WLY_Year1_hd)[4] <- "WLY_yearQ3" WLY_Year1_hd$WLY_yearQ3 <- WLY_Year1_hd$WLY_yearQ3/0.3 ## dry matter to fresh matter } if(ll %in% latlongQ1){ WLY_Year2 <- readRDS(file=paste("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2019/WLY_", ll, ".RDS", sep="")) WLY_Year2_hd <- WLY_Year2[, c("lat", "long", "plantingDate", hdays)] colnames(WLY_Year2_hd)[4] <- "WLY_yearQ2" WLY_Year2_hd$WLY_yearQ2 <- WLY_Year2_hd$WLY_yearQ2/0.3 } if(ll %in% latlongQ2){ WLY_Year3 <- readRDS(file=paste("/home/akilimo/projects/SoilNPK/lintul/Rwanda/LINTUL_WLY_2016/WLY_", ll, ".RDS", sep="")) WLY_Year3_hd <- WLY_Year3[, c("lat", "long", "plantingDate", hdays)] colnames(WLY_Year3_hd)[4] <- "WLY_yearQ1" WLY_Year3_hd$WLY_yearQ1 <- WLY_Year3_hd$WLY_yearQ1 / 0.3 } if(!is.null(WLY_Year2_hd) & !is.null(WLY_Year1_hd)){ WLY_2year <- merge(WLY_Year2_hd, WLY_Year1_hd, by=c("lat", "long", "plantingDate")) }else if(is.null(WLY_Year2_hd) & !is.null(WLY_Year1_hd)){ WLY_2year <- WLY_Year1_hd }else if(!is.null(WLY_Year2_hd) & is.null(WLY_Year1_hd)){ WLY_2year <- WLY_Year2_hd } if(!is.null(WLY_Year3_hd)){ WLY_3year <- merge(WLY_2year, WLY_Year3_hd, by=c("lat", "long", "plantingDate")) }else{ WLY_3year <- WLY_2year } data_long <- gather(WLY_3year, wly_years, measurement, names(WLY_3year)[4]:names(WLY_3year)[ncol(WLY_3year)], factor_key=TRUE) data_long <- ddply(data_long, .(lat, long, plantingDate), summarize, WLY=max(measurement)) colnames(data_long) <- c("lat","long","plantingDate", hdays) if(hdays==harvest8_15[1]){ WLY_8to15 <- data_long }else{ WLY_8to15 <- merge(WLY_8to15, data_long, by=c("lat","long","plantingDate") ) } } WLY_8to15 <- merge(WLY_8to15, Planting_Harvest[, c("startDate", "weekNr", "months")], by.x="plantingDate", by.y="startDate") return(WLY_8to15) } Max3Years_WLY <- do.call("rbind",lapply(all_locs[1:length(all_locs)], function(ll) WLY_3years(ll, Planting_Harvest = Planting_Harvest_RW, latlongQ3=latlongQ3, latlongQ1=latlongQ1, latlongQ2=latlongQ2))) str(Max3Years_WLY) Max3Years_WLY$location <- paste(Max3Years_WLY$lat, Max3Years_WLY$long, sep="_") length(unique(Max3Years_WLY$location)) head(Max3Years_WLY) summary(Max3Years_WLY) summary(Max3Years_WLY$WLY_361) unique(Max3Years_WLY[Max3Years_WLY$WLY_361 < 1000, ]$location) summary(Max3Years_WLY[Max3Years_WLY$location == "-2.625_28.925", ]) setwd("/home/akilimo/projects/SoilNPK/lintul/Rwanda") saveRDS(Max3Years_WLY, "Max3Years_WLY_Rwanda.RDS") Max3Years_WLY <- readRDS("Max3Years_WLY_Rwanda.RDS") summary(Max3Years_WLY) ### change the dry matter from LINTUL to fresh matter (/0.3) ## try to use Patricia's model on later stage