[R-SIG-Mac] Contour lines not plotting in correct position

[Duncan Murdoch]

On 12-10-24 10:42 PM, Thomas Veltman wrote:
> Hey gang!
>
> This is my first post to the list, but I will try to be as precise as possible.  I have googled around and can't find an easy explanation for why this is happening, so I figured I had better get some help from the experts, so thanks in advance!  My code is supplied below.  I am running R 2.15.1 on Mountain Lion, 64-bit.  OSX is 10.8.2
>
> I want to plot iso-efficiency lines for a hypothetical device.  I have a fully-parameterized expression, and have checked its sanity by hand, and it appears to operate as it ought to.  However, when I use the contour() function, the contour lines it produces don't lie in the correct place.  I generally believe their shape, but they appear to be translated by some amount.  I know that the analytical expression produces a given value, but when I look for that point on the plot, the iso-efficiency line doesn't cut through the point.  This is illustrated in the following code by the blue circle.  That point should fall on the red line, but clearly it does not.  Nevertheless, the governing equation (as defined in the R code) produces the correct result (if you check the array "Electrolyzer", the numbers match what you would get from simply applying the appropriate values to the expression, so it isn't a problem in the iterative population of the array).  I tried varying the mesh !
>   size (i.e. the number of points in the vectors etaELEC and FE), but of course this didn't have any effect besides smoothing the curves out a bit more.  I also thought that it might be a problem with the png output, but I tried generating the figure in the quartz window, and the error persists.  I also tried closing R, not saving the workspace, reopening, making a completely new working directory, and then re-entering the functions.  In all cases, the plot still ends up funky.
>
> Is there something fundamental that I am missing about the contour plot function?  Also, for those playing along at home, the code includes an "easy" way to add subscripts after special characters in a plot title (which is typically a useful feature).  Relevant code is pasted below:
>
>
> Ecap <- function(etaELEC,FE,etaFC,etaRP,ENH3,etaCOMP,etaST) {(etaELEC)/((1+((1-FE)/(FE*1.43))*(1-etaST*etaELEC+((1.2*1.43*etaELEC)/(etaRP*ENH3)))+((7.3*etaELEC)/(2*etaCOMP*ENH3))+((25*etaELEC)/(etaRP*ENH3))))}
>
>
> etaELEC <- seq(from=0.05,to=1, by=0.01)
> FE <- seq(from=0.05,to=1, by=0.01)
>
>
> Electrolyzer <- array(0,dim=c(length(etaELEC),length(FE)))
> for (i in 1:length(etaELEC)){
> for (j in 1:length(FE)){
> Electrolyzer[i,j]<-Ecap(FE=FE[j],etaST=0,etaFC=0,etaELEC=etaELEC[i],etaRP=0.6,etaCOMP=0.95,ENH3=340)}}
> png(file="Standalone.png",bg="transparent",height=5400,width=5400,res=864)
> par(family="serif",tck=0.01)
> plot(0, type="n", xlim=c(0,1),ylim=c(0,1), cex.lab=1.2, main="Energy Capture Coefficient for Standalone Device", xlab="Electrolyzer Voltage Efficiency", ylab="Faradaic Efficiency", cex.sub=0.75, sub=expression(paste(eta,""[FC],"=0, ",eta,""[Comp],"=0.95, ",eta,""[Sterling],"=0, ",eta,""[RP],"=0.60, "))) #make blank plot box with correct axes, etc...
> contour(x=seq(0,1,length.out=(nrow(Electrolyzer))),y=seq(0,1,length.out=(ncol(Electrolyzer))),Electrolyzer,axes=FALSE,levels=c(-0.1,0,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1))
> contour(x=seq(0,1,length.out=(nrow(Electrolyzer))),y=seq(0,1,length.out=(ncol(Electrolyzer))),Electrolyzer,axes=FALSE,levels=c(0.102),col="red")
> points(0.5,0.1547,col="blue")
> dev.off()
>
>

You use a different definition for the x in the call to contour, i.e. 
seq(0,1,length.out=(nrow(Electrolyzer))), then you used when you 
computed it, i.e. seq(from=0.05,to=1, by=0.01).

A good practice is to compute the x and y vectors just once, and use the 
outer() function to compute the array.  You may need to use Vectorize() 
to make sure your function works properly in outer().  For example,

Electrolyzer <- outer(etaELEC, FE, function(x, y) 
Ecap(FE=y,etaST=0,etaFC=0,etaELEC=x,etaRP=0.6,etaCOMP=0.95,ENH3=340))

then use etaELEC and FE in the calls to contour.  I would also use 
"add=TRUE" in your calls to contour() if you want those plots overlaid: 
  that will work on any device, not just png():

contour(etaELEC, 
FE,Electrolyzer,axes=FALSE,levels=c(-0.1,0,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), 
add=TRUE)
contour(etaELEC, 
FE,Electrolyzer,axes=FALSE,levels=c(0.102),col="red",add=TRUE)

Another good practice is to simplify your question as much as possible 
before posting it.  Often you'll discover your own error in the simple 
version, and won't need to post at all.

And finally, as David said, this sort of question should go to R-help, 
not R-sig-mac, unless you think there's something Mac-specific about the 
problem.  But posting a working example is much better than most "first 
posts"; thanks for doing that.

Duncan Murdoch