

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


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/
rsiggeobounces at stat.math.ethz.ch wrote on 20040714 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
>
> _______________________________________________
> RsigGeo mailing list
> RsigGeo at stat.math.ethz.ch
> https://www.stat.math.ethz.ch/mailman/listinfo/rsiggeo


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 superquick and uses nexttono 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).
rsiggeobounces at stat.math.ethz.ch wrote on 20040714 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 superquick and uses nexttono 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


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 superquick and uses nexttono 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

