# R function by Lindsay Clark <lvclark@illinois.edu>, February 17, 2015.
# Requires download of NAPPFAST Global Plant Hardiness Maps:
# http://www.nappfast.org/Plant_hardiness/2012/PLANT_HARDINESS_GeoTIFF.zip
# See http://www.nappfast.org/Plant_hardiness/2012/2012%20ph_index.htm
# for more information about the data set.
# This function creates a colored map of USDA plant hardiness zones for
# any geographic range worldwide.

# xlim and ylim should each be numeric vectors of length 2, giving the
# extent of longitude and latitude, respectively, to plot.
# GeoTIFF is the path to the GeoTiff file to use.
# res is how many pixels to plot in each dimension.
# col is a character vector of colors to use, in order of zone number
# starting with 1.
# database is passed to the maps function.
# cex is character expansion factor for the pixel that is drawn.
# legend is the position for a legend (see ?legend) or NA if no legend is drawn.

plotHardiness <- function(xlim, ylim,
                          GeoTIFF="~/PLANT_HARDINESS_GeoTIFF/PLANT_HARDINESS10YR.geo.tif",
                          res = round(abs(11 * (ylim[2] - ylim[1]))),
                          col = rgb(c(214,171,228,171,115, 69,172,
                                      236,224,235,239,212,133),
                                    c(211,172,144,105,162,184,215,
                                      216,183,155, 82, 90, 51),
                                    c(254,220,240,255,252, 67,108,
                                      127, 69, 58, 31, 87, 40), maxColorValue=255),
                          database = "worldHires", cex=0.5,
                          legend="bottomright"){
    require(raster)
    require(rgdal)
    require(maps)
    if(database %in% c("china","japan","nzHires","worldHires","world2Hires")){
        require(mapdata)
    }

    # import the geotiff file
    myraster <- raster(GeoTIFF)
    # get the points to extract from the dataset
    Xseq <- seq(xlim[1], xlim[2], length.out=res)
    Yseq <- seq(ylim[1], ylim[2], length.out=res)

    # draw initial map
    map(database, xlim=xlim, ylim=ylim)
    # get cells from raster object
    thesecells <- cellFromXY(myraster, cbind(rep(Xseq, each=res), rep(Yseq, times=res)))
    # get the hardiness zone for each point
    thesezones <- extract(myraster, thesecells)
    # plot the points
    points(rep(Xseq, each=res), rep(Yseq, times=res), pch=15, cex=cex,
           col=col[thesezones])
    # put map on top
    map(database, xlim=xlim, ylim=ylim, add=TRUE)
    # add legend
    if(!is.na(legend)){
        minzone <- min(thesezones, na.rm=TRUE)
        maxzone <- max(thesezones, na.rm=TRUE)
        legend(x=legend, legend=minzone:maxzone, fill=col[minzone:maxzone],
               title="Hardiness zones", bg="white")
    }
}