Hello,
I'm trying to use R to simulate biological invasions, and got stuck with the following: Is there a function that uses a matrix of occupied/empty values (1 or 0) as input, and producing an output matrix of the same size, with each cell containing the distance to the nearest occupied cell (1) of the input matrix ? For example, the input matrix: 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 Would produce an output matrix like: 1 0 0 1 1.4 1 0 1 2.2 1.4 1 1.4 2.8 2.2 2.0 2.2 For those familiar with AV Spatial Analyst, this function would be the equivalent of OutGrid = EucDistance (directionFN, allocationFN, maxDistance). Many thanks in advance for any help or pointers, Marius |
Marius Gilbert wrote:
> Hello, > > I'm trying to use R to simulate biological invasions, and got stuck with > the following: > > Is there a function that uses a matrix of occupied/empty values (1 or 0) > as input, and producing an output matrix of the same size, with each cell > containing the distance to the nearest occupied cell (1) of the input > matrix ? Try this, which uses an internal, undocumented Spatstat function 'exactPdt': library(spatstat) nmatdist <- function(aMat){ storage.mode(aMat) <- "logical" win <- owin(xrange=c(1,ncol(aMat)),yrange=c(1,nrow(aMat)),mask=aMat) m <- exactPdt(win) return(m$d) } > For example, the input matrix: > > 0 1 1 0 > 0 0 1 0 > 0 0 0 0 > 0 0 0 0 > > Would produce an output matrix like: > > 1 0 0 1 > 1.4 1 0 1 > 2.2 1.4 1 1.4 > 2.8 2.2 2.0 2.2 test: > m1 [,1] [,2] [,3] [,4] [,5] [1,] 0 0 0 0 0 [2,] 1 0 0 0 0 [3,] 1 1 0 0 0 [4,] 0 0 0 0 0 > nmatdist(m1) [,1] [,2] [,3] [,4] [,5] [1,] 1 1.414214 2.236068 2.828427 3.605551 [2,] 0 1.000000 1.414214 2.236068 3.162278 [3,] 0 0.000000 1.000000 2.000000 3.000000 [4,] 1 1.000000 1.414214 2.236068 3.162278 looks about right! Baz |
In reply to this post by Marius Gilbert
Here's another version that probably can be simplified: nmatdist <- function (m) { v <- as.vector (m) ind <- as.matrix (expand.grid (seq (nrow (m)), seq (ncol (m)))) d <- as.matrix (dist (ind, upper=TRUE)) ones <- as.numeric (dimnames(ind[v == 1,])[[1]]) matrix (sapply (seq(nrow(ind)), function (x) min (d[x, ones])), nrow=nrow(m)) } Denis White US EPA, 200 SW 35th St, Corvallis, Oregon, 97333 USA voice: 541.754.4476, email: white.denis at epa.gov web: www.epa.gov/wed/pages/staff/white/ r-sig-geo-bounces at stat.math.ethz.ch wrote on 2004-07-14 09:50:16: > Hello, > > I'm trying to use R to simulate biological invasions, and got stuck with > the following: > > Is there a function that uses a matrix of occupied/empty values (1 or 0) > as input, and producing an output matrix of the same size, with each cell > containing the distance to the nearest occupied cell (1) of the input > matrix ? > > For example, the input matrix: > > 0 1 1 0 > 0 0 1 0 > 0 0 0 0 > 0 0 0 0 > > Would produce an output matrix like: > > 1 0 0 1 > 1.4 1 0 1 > 2.2 1.4 1 1.4 > 2.8 2.2 2.0 2.2 > > For those familiar with AV Spatial Analyst, this function would be the > equivalent of OutGrid = EucDistance (directionFN, allocationFN, > maxDistance). > > Many thanks in advance for any help or pointers, > > Marius > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo at stat.math.ethz.ch > https://www.stat.math.ethz.ch/mailman/listinfo/r-sig-geo |
White.Denis at epamail.epa.gov wrote:
> > > > Here's another version that probably can be simplified: > > nmatdist <- function (m) > { > v <- as.vector (m) > ind <- as.matrix (expand.grid (seq (nrow (m)), seq (ncol (m)))) > d <- as.matrix (dist (ind, upper=TRUE)) > ones <- as.numeric (dimnames(ind[v == 1,])[[1]]) > matrix (sapply (seq(nrow(ind)), > function (x) min (d[x, ones])), nrow=nrow(m)) > } Eeek! Have you tried that on a 100x100 matrix? Darn near killed my machine! > mt=matrix(runif(10000)>.1,100,100) > nmatdist(mt) ... near swap death experience... The code used by Adrian Baddeley's spatstat routine uses a very neat method for working out the distances, which involves sweeping along rows and columns or something. He did explain it to me when I was in Perth but I can't recall it now! Anyway, its super-quick and uses next-to-no memory. Here's how long my function that calls the spatstat routine takes: > unix.time(nmatdist(mt)) [1] 0.02 0.01 0.02 0.00 0.00 it was so quick I thought I'd check I'd not done it on a small matrix by mistake: > dim(mt) [1] 100 100 Nope! Baz |
Yup. Space complexity is O((nrow(mat)*ncol(mat))^2). r-sig-geo-bounces at stat.math.ethz.ch wrote on 2004-07-14 10:49:14: > White.Denis at epamail.epa.gov wrote: > > > > > > > > Here's another version that probably can be simplified: > > > > nmatdist <- function (m) > > { > > v <- as.vector (m) > > ind <- as.matrix (expand.grid (seq (nrow (m)), seq (ncol (m)))) > > d <- as.matrix (dist (ind, upper=TRUE)) > > ones <- as.numeric (dimnames(ind[v == 1,])[[1]]) > > matrix (sapply (seq(nrow(ind)), > > function (x) min (d[x, ones])), nrow=nrow(m)) > > } > > Eeek! Have you tried that on a 100x100 matrix? Darn near killed my > machine! > > > mt=matrix(runif(10000)>.1,100,100) > > nmatdist(mt) > ... near swap death experience... > > The code used by Adrian Baddeley's spatstat routine uses a very neat > method for working out the distances, which involves sweeping along rows > and columns or something. He did explain it to me when I was in Perth > but I can't recall it now! > > Anyway, its super-quick and uses next-to-no memory. Here's how long my > function that calls the spatstat routine takes: > > > unix.time(nmatdist(mt)) > [1] 0.02 0.01 0.02 0.00 0.00 > > it was so quick I thought I'd check I'd not done it on a small matrix > by mistake: > > > dim(mt) > [1] 100 100 > > Nope! > > Baz > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo at stat.math.ethz.ch > https://www.stat.math.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Barry Rowlingson
Dear Barry,
> The code used by Adrian Baddeley's spatstat routine uses a very neat > method for working out the distances, which involves sweeping along rows > and columns or something. He did explain it to me when I was in Perth > but I can't recall it now! I knew there should be way to do it without estimating the full distance redistribution matrix O((nrow(mat)*ncol(mat))^2) because the equivalent AV Spatial analyst can do rather big raster (3000 x 3000 pixels) in a few seconds. > Anyway, its super-quick and uses next-to-no memory. Here's how long my > function that calls the spatstat routine takes: It is super quick indeed, I tested it on a 1600 x 250 matrix, ang got the output 1600x250 distance matrix in less than 5 seconds on a PIV 1700 Mhz, so fairly efficient indeed, Many thanks, there's no I could have found this myself, Cheers, Marius |
Free forum by Nabble | Edit this page |