[R-sig-Geo] from x,y plot to a raster with Raster package?

[Josep M Serra diaz]

Hi everyone,

I am wondering if I can speed up the process of transforming from  a 'plot'
(or matrix) to a raster.

The idea is that I have a matrix with x,y,z values. and I have a raster
with spatial data for x and for y, and I want to get a raster of z.

I would know how to do it on a cell-by-cell basis, but I was hoping there
was a more efficient way, since I have more than 14Milion Cells.

I left an simple example below to see the not very efficient way in which I
was dealing with it, hoping that someone will let me know that I have a
better chance trying with a function in the raster package that can do that
very easily. So far, I have tried with calc, but it is saying that 'I
cannot use this function'.




###### SIMPLE EXAMPLE

##dimensions of the matrix
mat <- list( x=c(23,25,27), y=c(33,34,35) , z=c(0.1,0.2,0.3))
mat2 <- matrix (data = mat$z,nrow=3, ncol=3,dimnames = mat[c(1,2)])

#input raster x axis
r1 <- raster ( nrows=10, ncols=10, crs=NA)
r1[] <- runif(n = 100,min=20,max=28)

#input raster y axis
r2 <- raster ( nrows=10, ncols=10, crs=NA)
r2[] <- runif(n = 100,min=31,max=36)


#AIM: I want a raster that has in the cell values the corresponding value
of the z dimension in mat.
#    That is r1 contains the values of the xaxis in mat, r2 contains inputs
values of the y axis in mat
#    I want the raster with the values of z in mat


#on a cell by cell basis it is 'easy'

cell.values <- lapply (1:ncell(r1),FUN = function (i){
  pos.x.axis <- if (r1[i] < range(mat$x)[1] | r1[i] > range(mat$x)[2]) {NA}
else{
    which(abs(r1[i]-mat$x) == min(abs(r1[i]-mat$x))) #it gives me the
closest postition in the y axis
    }

  pos.y.axis <- if (r2[i] < range(mat$y)[1] | r2[i] > range(mat$y)[2]) {NA}
else{
    which(abs(r2[i]-mat$x) == min(abs(r2[i]-mat$x))) #it gives me the
closest postition i nthe y axis
  }

  if (is.na (pos.x.axis) ==T | is.na(pos.y.axis)==T) {out.value <- NA} else{
    out.value <- mat2[pos.x.axis, pos.y.axis]
  }

  return (out.value)

})

raster.out <- r1 ; raster.out [] <- NA
raster.out[]<- unlist(cell.values)
plot (raster.out)


any hint on how to do it for 14Milion cells?

	[[alternative HTML version deleted]]