Friday, April 27, 2012

Create polygons from a matrix


The following function matrix.poly allows for the addition of polygons to a plot based on a matrix and defined matrix positions. I have used this function on occasion to highlight specific matrix locations (e.g. in the above figure). You can do the same by overlaying another image (left in above plot) but with this function you will have all other polygon plotting possibilities (e.g. borders etc.).
Another useful application is in the creation of polygons from a matrix representing a data field on a regular grid. The generated polygons can then be easily drawn onto maps of various projection.



matrix.poly function...
matrix.poly <- function(x, y, z=mat, n=NULL){
 if(missing(z)) stop("Must define matrix 'z'")
 if(missing(n)) stop("Must define at least 1 grid location 'n'")
 if(missing(x)) x <- seq(0,1,,dim(z)[1])
 if(missing(y)) y <- seq(0,1,,dim(z)[2])
 poly <- vector(mode="list", length(n))
 for(i in seq(length(n))){
  ROW <- ((n[i]-1) %% dim(z)[1]) +1
  COL <- ((n[i]-1) %/% dim(z)[1]) +1
 
  dist.left <- (x[ROW]-x[ROW-1])/2
  dist.right <- (x[ROW+1]-x[ROW])/2
  if(ROW==1) dist.left <- dist.right
  if(ROW==dim(z)[1]) dist.right <- dist.left
 
  dist.down <- (y[COL]-y[COL-1])/2
  dist.up <- (y[COL+1]-y[COL])/2
  if(COL==1) dist.down <- dist.up
  if(COL==dim(z)[2]) dist.up <- dist.down
 
  xs <- c(x[ROW]-dist.left, x[ROW]-dist.left, x[ROW]+dist.right, x[ROW]+dist.right)
  ys <- c(y[COL]-dist.down, y[COL]+dist.up, y[COL]+dist.up, y[COL]-dist.down)
  poly[[i]] <- data.frame(x=xs, y=ys)
 }
 return(poly)
}
Created by Pretty R at inside-R.org


to reproduce the above figures...
###
#required packages
library(maps)
library(mapproj)
 
#required functions (from "www.menugget.blogspot.com")
source("matrix.poly.R")
source("val2col.R")
 
png("matrix.poly.ex.png", width=9, height=3, units="in", res=200)
par(mfcol=c(1,3), mar=c(5,5,0.5,0.5))
set.seed(1)
 
#Ex 1 - add another image layer
m=8
n=10
x <- seq(m)
y <- seq(n)
z <- matrix(runif(m*n), nrow=m, ncol=n)
image(x,y,z, col=grey.colors(20))
N <- sample(1:(m*n),20)
z2 <- NaN*z
z2[N] <- 1
image(x,y,z2, col=rgb(0,0,1,0.4), add=TRUE)
box()
 
#Ex 2 - add polygons
m=8
n=10
x <- seq(m)
y <- seq(n)
z <- matrix(runif(m*n), nrow=m, ncol=n)
image(x,y,z, col=grey.colors(20))
N <- sample(1:(m*n),20)
poly <- matrix.poly(x,y,z=z,n=N)
sapply(poly, function(X){polygon(X, col=rgb(1,1,0,0.4), border=1)})
box()
 
#Ex 2 - add polygons to unequal grid
x <- cumsum(round(runif(m, min=1, max=10)))
y <- cumsum(round(runif(n, min=1, max=10)))
z <- matrix(runif(m*n), nrow=m, ncol=n)
image(x,y,z, col=grey.colors(20))
N <- sample(1:(m*n),20)
poly <- matrix.poly(x,y,z=z,n=N)
sapply(poly, function(X){polygon(X, col=rgb(1,0,0,0.4), border=1)})
box()
 
dev.off()
 
####
png("matrix.poly.ex2.png", width=6, height=3, units="in", res=400)
par(mfcol=c(1,2), mar=c(5,5,0.5,0.5))
set.seed(1)
 
#Ex 1 - simple image
m=10
n=10
x <- seq(from=-60, to=60,, m)
y <- seq(from=30, to=85,, n)
z <- matrix(runif(m*n), nrow=m, ncol=n)
image(x,y,z, col=rainbow(20))
 
#Ex.2 project polygons
par(mar=c(0,0,0,0))
map("world", proj="orthographic", orient=c(60,0,20), par=NULL, t="n")
poly <- matrix.poly(x,y,z=z,n=seq(length(z)))
COL <- val2col(z,col=rainbow(20))
for(i in seq(poly)){
 polygon(mapproject(poly[[i]]), col=COL[i], border=COL[i], lwd=0.3)
}
map("world", proj="orthographic", orient=c(60,0,20), par=NULL, add=T)
map.grid(col=rgb(0,0,0,0.5), labels=F)
 
dev.off()
Created by Pretty R at inside-R.org




2 comments:

  1. This is a nice function. I encountered an issue when a matrix defines a mask with only a single cell:
    testmat <- matrix(c(0,0,0, 0,1,0, 0,0,0),nrow=3) matrix.poly(c(1,2,3), c(1,2,3), testmat, which(testmat==1))
    # Error in if (ROW == 1) dist.left <- dist.right : missing value where TRUE/FALSE needed.

    I don't seem to be able to figure out how to avoid the error. Any hints?

    ReplyDelete
    Replies
    1. Thanks for finding this error. The problem was in the definition of indices in the for loop - It was "for(i in seq(n)){}", but should have been "for(i in seq(length(n))){".

      Delete