Bucket model

library("dismo")
library(EcoHydRology)
path<- "~/geoserver/data_dir/rasters/worldclim"
tmin<-getData(name = "worldclim", var = "tmin", res = 10, path=path)
tmax<-getData(name = "worldclim", var = "tmax", res = 10, path=path)
prec<-getData(name = "worldclim", var = "prec", res = 10, path=path)

clim<-stack(tmin,tmax,prec)
clim <- brick(clim)
names(clim) <- c(paste("mint", 1:12, sep = ""), paste("maxt", 1:12, sep = ""), 
    paste("prec", 1:12, sep = ""))

lat_rad <- coordinates(clim)[, 2] * pi/180

for (i in 1:12) {
    evap <- raster(clim, 1)
    Tmax <- values(subset(clim, i + 12))/10
    Tmin <- values(subset(clim, i))/10
    d <- data.frame(day = (30 * i) - 15, Tmin, Tmax, lat_rad)
    d[is.na(d)] <- 0
    Es_PET <- PET_fromTemp(Jday = d$day, Tmax_C = d$Tmax, Tmin_C = d$Tmin, lat_radians = d$lat_rad) * 
        1000
    values(evap) <- Es_PET
    if (i == 1) {
        PET <<- brick(evap)
    }
    if (i > 1) {
        PET <<- addLayer(PET, evap)
    }
}

months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", 
    "Nov", "Dec")
names(PET) <- months
plot(PET)

Bucket <- raster(PET, 1)

png("world_aet%02d.png", height = 800, width = 900)
for (n in 1:2) {
    for (i in 1:359) {
        mn <- 1 + i%/%30
        NewAET <- raster(PET, 1)
        NewBucket <- values(Bucket)
        rain <- values(subset(clim, 24 + mn))/30
        alpha <- (NewBucket - 200)/300
        evap <- values(subset(PET, mn)) * alpha * 0.8  ##A fudge factor for stomatal control.
        NewBucket <- NewBucket + (rain) - evap
        NewBucket[NewBucket > 500] <- 500
        NewBucket[NewBucket < 200] <- 200
        values(Bucket) <- NewBucket
        values(NewAET) <- evap * (NewBucket > 200)
        if (n > 1 && (i%%30) - 15 == 0) {
            print(plot(NewAET, main = months[mn]))
            if (mn == 1) {
                AET <<- brick(NewAET)
            }
            if (mn > 1) {
                AET <<- addLayer(AET, NewAET)
            }
        }
    }
}
dev.off()

plot(AET[[1]],ylim=c(-60,60))

writeRaster(AET,"world_AET.tif")
system("convert -delay 100 world*.png world_aet.gif")