## ----plantalgo1---------------------------------------------------------------
plant_harvest <- function(x, min_prec=30, max_prec=90, max_len=5) {

## check input
# must have 12 values (months)
	stopifnot(length(x) == 12)
# cannot have any NAs	
	if (any(is.na(x))) return (matrix(NA, nrow=12, ncol=2))

# max_len is the number of months
	stopifnot(max_len > 0 & max_len <= 12)
# min precipitation threshold below max threshold
	stopifnot(min_prec < max_prec)

## compute threshold 
# median of the rainfall
	med <- stats::median(x)
# median clamped between min_prec and max_prec	
	prec_threshold <- min(max_prec, max(min_prec, med))

# make a time series of 24 months, from July to July, to make it easier 
# to do 'circular' computation across the year boundary. That is, we 
# cannot ignore that (January comes after December)	
	x24 <- c(x[7:12], x, x[1:6])
	
# which months are above the threshold?	
	above_th <- x24 >= prec_threshold

# cumulate successive months above the threshold
	y <- cumsum(above_th)

# remove months below the threshold and reset 
# the beginning of a sequence to 1
# (perhaps the most obscure step)
	wet <- (y - cummax(y * !above_th))

# go back to a single 12 months (Jan to Dec) year
	wet <- wet[7:18]

# set up output	
	planting <- harvest <- rep(0, 12)	

# find the length of the growing season
	m <- min(max_len, max(wet))
	
# growing season must be at least 3 months	
	if (m > 2) {
	# harvest months	
		harvest[wet >= m] <- 1
	# planting months
		p <- which(wet >= m) - m + 1
	# 1 month before January -> December
		p[p < 1] <- 12 + p[p < 1]
		planting[p] <- 1
	}
	return( cbind(planting=planting, harvest=harvest) )
}


## ----plantalgo3---------------------------------------------------------------
rain <- c(0,0,10,40,50,60,60,50,40,10,0,0)
agro::plant_harvest(rain)


## ----printfun-----------------------------------------------------------------
printSeason <- function(rain) {
	x <- agro::plant_harvest(rain)
	rownames(x) <- month.abb
	t(x)
}


## ----printfun1----------------------------------------------------------------
printSeason(rain)


## ----plantalgofun-------------------------------------------------------------
plotSeason <- function(rain) {
	x <- agro::plant_harvest(rain)
	barplot(rain, names.arg=month.abb, las=2)
	barplot(x[,1]*rain, add=T, col="red", axes=FALSE, density=15)
	barplot(x[,2]*rain, add=T, col="blue", axes=FALSE, density=15, angle=315)
	legend("topright", c("rain", "plant", "harvest"), fill=c("gray", "red", "blue"), density=c(-1,25,25))
}


## ----plantalgo10--------------------------------------------------------------
plotSeason(rain)


## ----plantalgo20--------------------------------------------------------------
rain <- c(0,rep(50,4), 0, 0, rep(50,4), 0)
rain
printSeason(rain)
plotSeason(rain)


## ----plantalgo25--------------------------------------------------------------
rain <- c(4,4,5,5,6,6,6,5,5,4,4,4) * 10
rain
printSeason(rain)
plotSeason(rain)


## ----plantalgo30--------------------------------------------------------------
rain <- c(1:6, 6:1)
printSeason(rain)
plotSeason(rain)


## ----plantalgo100, message=FALSE----------------------------------------------
library(geodata)
rain <- geodata::worldclim_country("Senegal", var="prec", path=".")
adm <- geodata::gadm("Senegal", level=1, path=".")
rain <- mask(rain, adm)
names(rain) <- month.abb


## ----plantalgo105-------------------------------------------------------------
plot(sum(rain))
lines(adm)


## ----plantalgo107, fig.width=10, fig.height=7---------------------------------
plot(rain)


## ----plantalgo110-------------------------------------------------------------
arain <- aggregate(rain, 10, mean, na.rm=TRUE)


## ----plantalgo120-------------------------------------------------------------
x <- app(arain, agro::plant_harvest)
x


## ----plantalgo122-------------------------------------------------------------
plant <- x[[1:12]]
names(plant) <- month.abb
harv <- x[[13:24]]
names(harv) <- month.abb


## ----plantalgo130, fig.width=10, fig.height=7---------------------------------
plot(plant)


## ----plantalgo160-------------------------------------------------------------
plantna <- classify(plant, cbind(0, NA))
p <- which.max(plantna)
plot(p)


## ----plantalgo180-------------------------------------------------------------
harvna <- classify(harv, cbind(0, NA))
h <- which.max(harvna)
plot(h)