[R-sig-eco] colorkey for image plots

Tomas tdomingu at staffmail.ed.ac.uk
Wed Feb 4 16:20:56 CET 2009


Hi
I'm developing a model of leaf photosynthesis based on leaf nitrogen (N) 
and phosphorus (P). I've created a figure with a large range of possible 
combinations of N and P (x and y axes) with the outputs of my model as 
colors from the image function.
I couldn't include on the figure a color key that would give the range 
of the model outputs. I tried other functions (filled.contour, 
levelplot, wireframe) but couldn't get around the need for ascending 
order of values.

My code follows below. I'd appreciate any comments. Also, ideas on how 
to include a third leaf trait on this.

Cheers,
Tomas

###################################
# figure with the variation of Vcmax with N and P

rm(list=ls(all=TRUE))
detach()

##the libraries
library(akima)
library(stats)

## the data
data1 <- read.csv(file.name <- choose.files(),na.strings = "NA" , header 
= TRUE)
attach(data1)
# the following prints the names of the variables from the dataset, 
check if file was correctly loaded
names(data1)


n<-seq(0,4,length=25)
p<-seq(0,0.4,length=26)
N<-rep(n,length(p))
P<-rep(p,length(n))


vn <- lm(Vcmax_Ci_area_25C~Narea_gm.2)
vp <- lm(Vcmax_Ci_area_25C~Parea_gm.2)
np<-lm(Narea_gm.2~Parea_gm.2)
jn <- lm(Jmax_Ci_area_25C~Narea_gm.2)
jp <- lm(Jmax_Ci_area_25C~Parea_gm.2)
vnp <- lm(Vcmax_Ci_area_25C~Narea_gm.2*Parea_gm.2)



#############################################################
# model of Vcmax_Ci_area_25C based on both N and P
fitVJ <- function(x){
	a <- x[1]		# variable 1 (Intercept of Vcmax_Ci_area_25C~N.area.g.m-2)
	b <- x[2]		# variable 2 (Slope of Vcmax_Ci_area_25C~N.area.g.m-2)
	c <- x[3]		# variable 3 (Intersept of Vcmax_Ci_area_25C~P.area.g.m-2)
	d <- x[4]		# variable 4 (Slope of Vcmax_Ci_area_25C~P.area.g.m-2)
	e <- x[5]		# variable 1 (Intercept of Vcmax_Ci_area_25C~N.area.g.m-2)
	f <- x[6]		# variable 2 (Slope of Vcmax_Ci_area_25C~N.area.g.m-2)
	g <- x[7]		# variable 3 (Intersept of Vcmax_Ci_area_25CC~P.area.g.m-2)
	h <- x[8]		# variable 4 (Slope of Vcmax_Ci_area_25C~P.area.g.m-2)
	
	Vmodeled <- pmin((a+b*(Narea_gm.2)),(c+d*(Parea_gm.2)))
	Jmodeled <- pmin((e+f*(Narea_gm.2)),(g+h*(Parea_gm.2)))
	Vdiff <- sum(((Vcmax_Ci_area_25C) - Vmodeled)^2 + ((Jmax_Ci_area_25C) 
- Jmodeled)^2)			
	}
guessvj <- c(coef(vn)[[1]],coef(vn)[[2]], 
coef(vp)[[1]],coef(vp)[[2]],coef(jn)[[1]],coef(jn)[[2]], 
coef(jp)[[1]],coef(jp)[[2]])
g<- c(1,1,1,1,1,1,1,1)
Fitv<- optim(guessvj,fitVJ,control = list (maxit =99000), method = "BFGS")

V <- pmin((Fitv$par[1]+Fitv$par[2]*N),(Fitv$par[3]+Fitv$par[4]*P))
V2 <- pmin((Fitv$par[5]+Fitv$par[6]*N),(Fitv$par[7]+Fitv$par[8]*P))


par(mfrow=c(1,2))
D1<-interp(P,N,V)
image(D1,col = terrain.colors(50),xlim=c(0,0.35),ylim=c(0,3.9))
title("Vcmax-modelled");title(xlab="Parea, g m-2"); title(ylab="Narea, g 
m-2")
points(Narea_gm.2~Parea_gm.2)

D2<-interp(P,N,V2)
image(D2,col = terrain.colors(50),xlim=c(0,0.35),ylim=c(0,3.9))
title("Jmax-modelled");title(xlab="Parea, g m-2"); title(ylab="Narea, g 
m-2")
points(Narea_gm.2~Parea_gm.2)

### end ###

-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



More information about the R-sig-ecology mailing list