Commit 4ae6642b authored by Rich FitzJohn's avatar Rich FitzJohn
Browse files

First version of package

parent 5ec5b812
Package: traitcompetition
Title: Trait Competition Workshop
Version: 0.1
Date: 2014-05-21
Author: Georges Kunstler
Maintainer: Georges Kunstler <georgesk@gmail.com>
Description: Code for the trait competition workshop
License: GPL (>=2)
Imports:
sp,
boot
PACKAGE := $(shell grep '^Package:' DESCRIPTION | sed -E 's/^Package:[[:space:]]+//')
all:
make -C src
document: roxygen
roxygen:
@mkdir -p man
Rscript -e "library(methods); devtools::document()"
install: roxygen
R CMD INSTALL .
clean:
make -C src clean
build:
R CMD build .
check: build
R CMD check --no-manual `ls -1tr ${PACKAGE}*gz | tail -n1`
@rm -f `ls -1tr ${PACKAGE}*gz | tail -n1`
@rm -rf ${PACKAGE}.Rcheck
# test:
# make -C tests/testthat
.PHONY: attributes document roxygen install clean build check
# Generated by roxygen2 (4.0.1): do not edit by hand
S3method(plot,biome.map)
S3method(plot,ellips.on.biome.map)
S3method(plot,points.on.biome.map)
S3method(plot,points.on.biome.map.diff.tropical)
S3method(plot,points.on.biome.map.ecocode)
S3method(plot,points.on.biome.map.panel)
export(GetClimate)
export(NearestPlot)
export(axis.log10)
export(f.quantile)
export(f.quantile.boot)
export(fun.boxplot.breaks)
export(fun.circles.plot)
export(fun.convert.type.B)
export(fun.convert.type.I)
export(fun.ellipsoid.hull)
export(fun.overly.plot.on.biomes)
export(fun.plot.colors.density)
export(fun.points.colors.density)
export(fun.quadrat)
export(fun.quadrat.cluster)
export(fun.sp.in.ecoregion)
export(fun_div)
export(function.perc.dead)
export(function.perc.dead2)
export(is.wholenumber)
export(label)
export(make.transparent)
export(na.clean)
export(niceColors)
export(prepare.plotting.info)
export(to.dev)
export(to.pdf)
export(trim.negative.growth)
export(trim.positive.growth)
export(trim.trailing)
import(sp)
###################################################
###################################################
###################################################
##### FUNCTION TO FORMAT DATA FOR THE WORKSHOP
#### G. Kunstler 11/09/2013
############################
## FUNCTION remove trailing white space
##' @export
trim.trailing <- function (x) sub("\\s+$", "", x)
#############################################
## FUN to clean species name for FRENCH NFI
fun.clean.species.tab <- function(species.tab){
species.tab2 <- species.tab[!is.na(species.tab$Latin_name), ]
### species IFN reformat names
## clean species names and synonyme names
species.tab2$Latin_name <- (gsub("_", " ", species.tab2$Latin_name))
species.tab2$Latin_name_syn<- (gsub("_", " ", species.tab2$Latin_name_syn))
## remove trailing white space
species.tab2$Latin_name_syn<- trim.trailing(species.tab2$Latin_name_syn)
species.clean <- species.tab2[!duplicated(species.tab2$Latin_name),
c("code", "Latin_name",
"Exotic_Native_cultivated")]
names(species.clean) <- c("sp", "Latin_name", "Exotic_Native_cultivated")
return(species.clean)
}
#####################################################
### compute quantile 99% and sd with a bootstrap
##' .. compute quantile 99% and sd with a bootstra..
##'
##' @title f.quantile.boot
##' @param i subset (species) for which to compute quantile
##' @param x vector of height
##' @param fac vector of subset
##' @param R number of resampling
##' @param probs probability of quantile to compute
##' @return a matrix with # col and 1 row with mean sd adn nobs
##' @author Kunstler
##' @export
f.quantile.boot <- function(i, x, fac, R, probs = 0.99){
if(length( na.exclude(x[fac == i])) > 0 ){
quant.boot <- boot::boot(x[fac == i], f.quantile, R = R, probs = probs)
return(as.matrix(c(mean = mean(quant.boot$t), sd = sd(quant.boot$t),
nobs = length(na.exclude(x[fac == i]))),
ncol = 3, nrow = 1))
}else{
return(as.matrix(c(mean = NA, sd = NA, nobs = NA), ncol = 3, nrow = 1))
}
}
##' @export
f.quantile <- function (x, ind, probs){
quantile(x[ind], probs = probs, na.rm = TRUE)
}
#######################
### function to compute number of dead per plot
##' @export
function.perc.dead <- function(dead){
sum(dead)/length(dead)} # == mean(dead)
## function to deal with missing value
##' @export
function.perc.dead2 <- function(dead) {
out <- sum(dead, na.rm = T)/length(dead[!is.na(dead)])
if(!is.finite(out)) out <- NA
return(out)
}
## function to convert var factor in character or numeric
##' @export
fun.convert.type.I <- function(data.tree){
character.var <- c("sp", "sp.name", "ecocode")
numeric.var <- c("D", "G", "BA.G", "dead",
"Lon", "Lat", "weights", "MAT", "MAP")
factor.to.convert.var <- c("obs.id", "tree.id", "cluster", "plot", "census")
for (i in factor.to.convert.var)
data.tree[[i]] <- as.integer(factor(data.tree[[i]]))
for (i in character.var)
data.tree[[i]] <- as.character(data.tree[[i]])
for (i in numeric.var)
data.tree[[i]] <- as.numeric(data.tree[[i]])
return(data.tree)
}
##' @export
fun.convert.type.B <- function(data.tree){
character.var <- c("sp", "sp.name", "ecocode")
numeric.var <- c("D", "G", "BA.G", "dead",
"x", "y")
factor.to.convert.var <- c("obs.id", "tree.id", "cluster", "plot", "census")
for (i in factor.to.convert.var) data.tree[[i]] <- as.integer(factor(data.tree[[i]]))
for (i in character.var) data.tree[[i]] <- as.character(data.tree[[i]])
for (i in numeric.var) data.tree[[i]] <- as.numeric(data.tree[[i]])
return(data.tree)
}
######
######
## FUNCTION TO PLOT MAP OF TREE with function of DBH
##' .. Function to plot map of tree with circle function of their dbh..
##'
##' @title Plot Tree
##' @param plot.select the plot for which draw the map
##' @param x
##' @param y
##' @param plot vectore of plot id for each tree
##' @param D diameter in cm
##' @param inches controling the circle size
##' @param fg.l color circle
##' @param ...
##' @author Kunstler
##' @export
fun.circles.plot <- function(plot.select, x, y, plot, D, inches, fg.l, ...){
x.t <- x[plot == plot.select]
y.t <- y[plot == plot.select]
D.t <- D[plot == plot.select]
fg <- fg.l[plot == plot.select]
D.t[is.na(D.t)] <- 0
symbols(x.t, y.t, circles = D.t , main = plot.select, inches = inches,
fg = fg, ...)
}
##' .. remove too negative growth based on Condit R package
##' with param fitted to BCI ..
##'
##' .. taken from trim.growth function in CTFS.R ..
##' @title trim.negative.growth
##' @param dbh1 in mm
##' @param dbh2 in mm
##' @param slope not to be changed
##' @param intercept
##' @param err.limit
##' @return a vector TRUE FALSE with FALSE outlier to be removed
##' @author Kunstler
##' @export
trim.negative.growth <- function(dbh1, dbh2, slope = 0.006214,
intercept = .9036, err.limit = 5){
stdev.dbh1 <- slope*dbh1 + intercept
bad.grow <- which(dbh2<=(dbh1-err.limit*stdev.dbh1))
accept <- rep(TRUE, length(dbh1))
accept[bad.grow] <- FALSE
accept[is.na(dbh1) | is.na(dbh2) | dbh2<=0 | dbh1<=0] <- FALSE
return(accept)
}
##' .. remove too high growth ..
##'
##' .. taken from trim.growth in Condit CTFS R package ..
##' @title trim.positive.growth
##' @param growth in mm
##' @param maxgrow in mm
##' @return TRUE FALSE vector with FALSE outlier
##' @author Kunstler
##' @export
trim.positive.growth <- function(growth, maxgrow = 75){
bad.grow <- which(growth>maxgrow)
accept <- rep(TRUE, length(growth))
accept[bad.grow] <- FALSE
accept[is.na(growth)] <- FALSE
return(accept)
}
##################### CREATE PLOT BASED ON square.sxsquare.s m cell from X Y
##' @export
fun.quadrat <- function(x, y, square.s) {
if(sum(!is.na(x))>10){
vec.x <- seq(0, max(x, na.rm = T) + 20, by = square.s)
vec.y <- seq(0, max(y, na.rm = T) + 20, by = square.s)
x.cut <- cut(x, breaks = vec.x, include.lowest = TRUE)
y.cut <- cut(y, breaks = vec.y, include.lowest = TRUE)
out <- apply(cbind(x.cut, y.cut), 1, paste, collapse = ".")
out[is.na(x.cut)] <- NA
out[is.na(y.cut)] <- NA
return(unclass(out))
}else{
return(rep(NA, length(x)))
}
}
## function to apply per cluster
##' @export
fun.quadrat.cluster <- function(cluster.id, cluster, tree.id, x, y, square.s){
temp <- fun.quadrat(x[cluster == cluster.id], y[cluster == cluster.id],
square.s = square.s)
temp2 <- paste((cluster[cluster == cluster.id]), temp, sep = ".")
temp2[is.na(temp)] <- NA
return(data.frame(make.quad = temp2, tree.id = tree.id[cluster == cluster.id]))
}
##' @export
fun_div <- function(sp){
p_i <- table(factor(sp))/sum(table(factor(sp)))
shannon <- -sum(p_i*log(p_i))
simpson <- 1-sum(p_i^2)
return(data.frame(shannon=shannon,simpson=simpson))
}
##' @export
fun.sp.in.ecoregion <- function(data, thres.prop=0.05){
names(table(data[["sp"]]))[table(data[["sp"]])/length(data[["sp"]]) > thres.prop]
}
#Returns the index of the nearest plot(x2,y2) to each of the target plots (x1,y1)
##' @export
NearestPlot <- function(x1,y1,x2,y2) {
res = numeric(length(x1))
for (i in 1:length(x1)) {
dists = sqrt((x1[i] - x2)^2 + (y1[i] - y2)^2)
res[i] = which.min(dists)
}
return(res)
} #NearestPlot
#looks up MAT and MAP for given lat/lon values
##' @export
GetClimate <-function(lats,lons) {
#create raster grid for tiles
xNum<-0:11
yNum<-0:4
xLon<-seq(from=-165,to=165,by=30)
yLat<-seq(from=75,to=-45,by=-30)
NumMat<-outer(yNum,xNum,paste,sep="")
NumLookup<-as.vector(NumMat)
NumRast <- raster(matrix(1:60,nrow=5),xmn=-180,xmx=180,ymn=-60,ymx=90)
LonLatLookup<-expand.grid(Lat=yLat,Lon=xLon)
row.names(LonLatLookup) <- NumMat
#get the set of tiles corresponding to the lons/lats arguments
tiles <- extract(NumRast,cbind(lons,lats))
unique.tiles <- unique(tiles)
plot.mat <- numeric(length(tiles))
plot.map <- numeric(length(tiles))
#download each of the tiles
for (i.tile in unique.tiles) {
temp.tile <- getData('worldclim',var='bio',res=0.5,lon=LonLatLookup[i.tile,"Lon"],lat=LonLatLookup[i.tile,"Lat"])
plot.mat[tiles==i.tile] <- extract(temp.tile,cbind(lons[tiles==i.tile],lats[tiles==i.tile]),layer=1,nl=1)/10 #retrieve MAT, divide by 10 for deg C
plot.map[tiles==i.tile] <- extract(temp.tile,cbind(lons[tiles==i.tile],lats[tiles==i.tile]),layer=12,nl=1) #retrieve MAP
nas <- which(tiles==i.tile & is.na(plot.mat))
#assign climate of nearest plot to all nas
if (length(nas)>0) {
good <- which(tiles==i.tile & !is.na(plot.mat))
near <- NearestPlot(lons[nas],lats[nas],lons[good],lats[good])
plot.mat[nas] <- plot.mat[good[near]]
plot.map[nas] <- plot.map[good[near]]
}
}
return(list(MAT=plot.mat,MAP=plot.map))
} #GetClimate
##' @export
to.pdf <- function(expr, filename, ..., verbose=TRUE) {
if(!file.exists(dirname(filename)))
dir.create(dirname(filename), recursive=TRUE)
if ( verbose )
cat(sprintf("Creating %s\n", filename))
pdf(filename, ...)
on.exit(dev.off())
eval.parent(substitute(expr))
}
##' @export
to.dev <- function(expr, dev, filename, ..., verbose=TRUE){
if(!file.exists(dirname(filename)))
dir.create(dirname(filename), recursive=TRUE)
if ( verbose )
cat(sprintf("Creating %s\n", filename))
dev(filename, ...)
on.exit(dev.off())
eval.parent(substitute(expr))
}
##' @export
na.clean <-function(x){x[!is.na(x)]}
#returns up to 80 unique, nice colors,
#generated using http://tools.medialab.sciences-po.fr/iwanthue/
# Starts repeating after 80
##' @export
niceColors<-function(n=80){
cols<-rep(c("#75954F","#D455E9","#E34423","#4CAAE1","#451431","#5DE737",
"#DC9B94","#DC3788","#E0A732","#67D4C1","#5F75E2","#1A3125",
"#65E689","#A8313C","#8D6F96","#5F3819","#D8CFE4","#BDE640",
"#DAD799","#D981DD","#61AD34","#B8784B","#892870","#445662",
"#493670","#3CA374","#E56C7F","#5F978F","#BAE684","#DB732A",
"#7148A8","#867927","#918C68","#98A730","#DDA5D2","#456C9C",
"#2B5024","#E4D742","#D3CAB6","#946661","#9B66E3","#AA3BA2",
"#A98FE1","#9AD3E8","#5F8FE0","#DF3565","#D5AC81","#6AE4AE",
"#652326","#575640","#2D6659","#26294A","#DA66AB","#E24849",
"#4A58A3","#9F3A59","#71E764","#CF7A99","#3B7A24","#AA9FA9",
"#DD39C0","#604458","#C7C568","#98A6DA","#DDAB5F","#96341B",
"#AED9A8","#55DBE7","#57B15C","#B9E0D5","#638294","#D16F5E",
"#504E1A","#342724","#64916A","#975EA8","#9D641E","#59A2BB",
"#7A3660","#64C32A"),
ceiling(n/80))
cols[1:n]
}
##' @export
make.transparent <- function(col, opacity=0.5) {
tmp <- col2rgb(col)/255
rgb(tmp[1,], tmp[2,], tmp[3,], alpha=opacity)
}
## Position label at a fractional x/y position on a plot
##' @export
label <- function(px, py, lab, ..., adj=c(0, 1)) {
usr <- par("usr")
text(usr[1] + px*(usr[2] - usr[1]),
usr[3] + py*(usr[4] - usr[3]),
lab, adj=adj, ...)
}
##' @export
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5)
abs(x - round(x)) < tol
##' @export
axis.log10 <- function(side=1, horiz=FALSE, labels=TRUE, baseAxis = TRUE,
wholenumbers=T, labelEnds=T,las=1, at=NULL) {
fg <- par("fg")
if(is.null(at)){
#get range on axis
if(side ==1 | side ==3) {
r <- par("usr")[1:2] #upper and lower limits of x-axis
} else {
r <- par("usr")[3:4] #upper and lower limits of y-axis
}
#make pertty intervals
at <- pretty(r)
#drop ends if desirbale
if(!labelEnds)
at <- at[at > r[1] & at < r[2]]
}
#restrict to whole numbers if desriable
if(wholenumbers)
at<-at[is.wholenumber(at)]
lab <- do.call(expression, lapply(at, function(i) bquote(10^.(i))))
#convert at if
if(baseAxis)
at<-10^at
#make labels
if ( labels )
axis(side, at=at, lab, col=if(horiz) fg else NA,
col.ticks=fg, las=las)
else
axis(side, at=at, FALSE, col=if(horiz) fg else NA,
col.ticks=fg, las=las)
}
##
##' @export
fun.plot.colors.density <- function(x,y,...){
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
plot(x,y, col=colors.dens, pch=20,...)
}
##' @export
fun.points.colors.density <- function(x,y,...){
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
points(x,y, col=colors.dens, pch=20,...)
}
##' @export
fun.boxplot.breaks <- function(x,y,Nclass=15,add=TRUE,...){
x1 <- x[!is.na(x) & !is.na(y)]
y1 <- y[!is.na(x) & !is.na(y)]
x <- x1
y <- y1
breakss <- seq(min(x),max(x),length=Nclass+1)
x.cut <- cut(x,breakss)
mid.points <- breakss[-(Nclass+1)]+abs(breakss[2]-breakss[1])/2
if(add){boxplot(y~x.cut,at=mid.points,outline=FALSE,add=add,names=NULL,
col=grey(1-table(x.cut)/length(x.cut)),
boxwex=(breakss[2]-breakss[1])*0.9,
pars=list(xaxt='n',axes=FALSE),...)
}else{
boxplot(y~x.cut,at=mid.points,outline=FALSE,names=round(mid.points,0),
boxwex=(breakss[2]-breakss[1])*0.9,
col=grey(1-table(x.cut)/length(x.cut)),...)
}
}
## function to overlay points on biomes
##' @export
fun.overly.plot.on.biomes <- function(MAP, MAT, names.vec, poly = poly.DF){
xy = cbind(MAP, MAT)
dimnames(xy)[[1]] = names.vec
pts = SpatialPoints(xy)
return(over(pts, poly))
}
## We need to load sp here for the plot.poly method.
##' @export
##' @import sp
plot.biome.map <- function(poly = poly.DF){
# MAKE EMPTY PLOT
plot(0,0, type = 'n', xlim = c(0, 450), ylim = c(30, -15),
xlab = "Mean annual precipitation (cm)",
ylab = "Mean annual temperature (deg C)")
# MAKE POLYGONS
plot(poly,col = make.transparent(c("navajowhite3", "darkgoldenrod1",
"sienna","darkolivegreen4",
"darkseagreen3", "forestgreen",
"darkgreen","olivedrab","gray"),
0.9),
border = NA,
add = TRUE)
}
##' @export
plot.points.on.biome.map.diff.tropical <- function(MAP, MAT, set,ecocode,
cols = niceColors(),
add.legend=TRUE){
plot.biome.map()
col <- c()
pch <- c()
cex <- c()
sets <- unique(set)
for(i in 1:length(sets)){
j <- (set == sets[i])
if(unique(ecocode[j])[1]=='tropical') {
col <- c(col,make.transparent(cols[i],0.75))
pch <- c(pch,16)
cex <- c(cex,1.25)
points(MAP[j], MAT[j], pch = pch[i], col = col[i], cex=cex[i])
}else{
col <- c(col,make.transparent(cols[i],0.25))
pch <- c(pch,16)
cex <- c(cex,.25)
points(MAP[j], MAT[j], pch = pch[i], col = col[i], cex=cex[i])
}
}
if(add.legend)
legend("topright", pch = pch, legend = sets, bty = "n", col = col, cex = 1)
}
##' @export
plot.points.on.biome.map.ecocode <- function(MAP, MAT, MAP.sd, MAT.sd, set,
ecocode, cols = niceColors(),
add.legend = TRUE){
plot.biome.map()
col <- c()
pch <- c()
cex <- c()
sets <- unique(set)
for(i in 1:length(sets)){
j <- (set == sets[i])
col <- c(col,make.transparent(cols[i],0.99))
pch <- c(pch,16)
cex <- c(cex,1.25)
points(MAP[j], MAT[j], pch = pch[i], col = col[i], cex=cex[i])