Graph to compare two matrices in R

I have two matrices (of approximately 300 x 100) and I would like to plot a graph to see the parts of the first one that are higher than those of the second.

I can do, for instance:

# Calculate the matrices and put them into m1 and m2
# Note that the values are between -1 and 1
par(mfrow=c(1,3))
image(m1, zlim=c(-1,1))
image(m2, zlim=c(-1,1))
image(m1-m2, zlim=c(0,1))

This will plot only the desired regions in the 3rd plot but I would like to do something a bit different, like putting a line ar开发者_StackOverflow社区ound those areas over the first plot in order to highlight them directly there.

Any idea how I can do that?

Thank you nico

How about:

par(mfrow = c(1, 3))
image(m1, zlim = c(-1, 1))
contour(m1 - m2, add = TRUE)
image(m2, zlim = c(-1, 1))
contour(m1 - m2, add = TRUE)
image(m1 - m2, zlim = c(0, 1))
contour(m1 - m2, add = TRUE)

This adds a contour map around the regions. Sort of puts rings around the areas of the 3rd plot (might want to fiddle with the (n)levels of the contours to get fewer 'circles').

Another way of doing your third image might be:

image(m1>m2)

this produces a matrix of TRUE/FALSE values which gets imaged as 0/1, so you have a two-colour image. Still not sure about your 'putting a line around' thing though...

Here's some code I wrote to do something similar. I wanted to highlight contiguous regions above a 0.95 threshold by drawing a box round them, so I got all the grid squares above 0.95 and did a clustering on them. Then do a bit of fiddling with the clustering output to get the rectangle coordinates of the regions:

computeHotspots = function(xyz, thresh, minsize=1, margin=1){
### given a list(x,y,z), return a data frame where each row
### is a (xmin,xmax,ymin,ymax) of bounding box of a contiguous area
### over the given threshhold.
### or approximately. lets use the clustering tools in R...

  overs <- which(xyz$z>thresh,arr.ind=T)

  if(length(overs)==0){
    ## found no hotspots
    return(NULL)
  }

  if(length(overs)==2){
    ## found one hotspot
    xRange <- cbind(xyz$x[overs[,1]],xyz$x[overs[,1]])
    yRange <- cbind(xyz$y[overs[,2]],xyz$y[overs[,2]])
  }else{

    oTree <- hclust(dist(overs),method="single")
    oCut <- cutree(oTree,h=10)

    oXYc <- data.frame(x=xyz$x[overs[,1]],y=xyz$y[overs[,2]],oCut)

    xRange <- do.call("rbind",tapply(oXYc[,1],oCut,range))
    yRange <- do.call("rbind",tapply(oXYc[,2],oCut,range))

  }

### add user-margins
 xRange[,1] <- xRange[,1]-margin
 xRange[,2] <- xRange[,2]+margin
 yRange[,1] <- yRange[,1]-margin
 yRange[,2] <- yRange[,2]+margin

## put it all together
 xr <- apply(xRange,1,diff)
 xm <- apply(xRange,1,mean)
 xRange[xr<minsize,1] <- xm[xr<minsize]-(minsize/2)
 xRange[xr<minsize,2] <- xm[xr<minsize]+(minsize/2)

 yr <- apply(yRange,1,diff)
 ym <- apply(yRange,1,mean)
 yRange[yr<minsize,1] <- ym[yr<minsize]-(minsize/2)
 yRange[yr<minsize,2] <- ym[yr<minsize]+(minsize/2)

  cbind(xRange,yRange)

}

Test code:

x=1:23
y=7:34
m1=list(x=x,y=y,z=outer(x,y,function(x,y){sin(x/3)*cos(y/3)}))
image(m1)
hs = computeHotspots(m1,0.95)

That should give you a matrix of rectangle coordinates:

> hs
  [,1] [,2] [,3] [,4]
1   13   15    8   11
2    3    6   17   20
3   22   24   18   20
4   13   16   27   30

Now you can draw them over the image with rect:

image(m1)
rect(hs[,1],hs[,3],hs[,2],hs[,4])

and to show they are where they should be:

image(list(x=m1$x,y=m1$y,z=m1$z>0.95))
rect(hs[,1],hs[,3],hs[,2],hs[,4])

You could of course adapt this to draw circles, but more complex shapes would be tricky. It works best when the regions of interest are fairly compact.

Barry