A belated thanks to Stephen Kaluzny and Nick Ellis for
pointing out that the check.islands function in S+Spatial
Stats will not do what I was trying to do. A special thanks
to Nick Ellis for providing me a function that is able
detect spatial units that have no neighbors. Their reply's
are given below but first here is the original posting:
> My problem is that I would like to detect spatial units
that
> have no neighbors using the
> "check.islands" function that comes with the
> S+SpatialStats module. For example a matrix
> x <- matrix(1,ncol=5,nrow=4)
> x[3,] <- x[,4] <- 0
> x
> [,1] [,2] [,3] [,4] [,5]
> [1,] 1 1 1 0 1
> [2,] 1 1 1 0 1
> [3,] 0 0 0 0 0
> [4,] 1 1 1 0 1
> has four "islands" (1-land , 0-water). In general, how can
> one calculate the number of islands
> for any configuration stored in a matrix using the
> S+SpatialStats module?
Stephen wrote:
The check.islands function requires a spatial.neighbor
object. If you
have a matrix of neighbor weights then you can construct a
spatial.neighbor object using the function spatial.neighbor.
A matrix of
neighbor weights must be a square matrix such that if
element [i,j] is
non-zero, then spatial regions j is considered a neighbor of
region i.
The neighbor relation is not considered symmetric by
default; [i,j] in
the neighbor weights matrix being non-zero means that j is a
neighbor of
i but unless [j,i] is non-zero in the neighbor weights
matrix, i is not a
neighbor of j.
The check.islands function shows those regions that have no
neighbors
AND are not neighbors of any other region. (The summary
method for
objects of class spatial.neighbor also shows this
information).
Your matrix is not square so I do not know what the value
represent.
If your x was:
> x <- matrix(1,ncol=5,nrow=5) # a square matrix
> x[3,] <- x[,4] <- 0
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 0 1
[2,] 1 1 1 0 1
[3,] 0 0 0 0 0
[4,] 1 1 1 0 1
[5,] 1 1 1 0 1
Then you can create a spatial.neighbor object:
> snx <- spatial.neighbor(neighbor.matrix=x)
But this does not have any islands:
> summary(snx)
Matrix was NOT defined as symmetric
Number of Regions: 5
Average Number of Connections: 4
Average Weight: 1
Least Number of Connections: 4 for Regions with
Indices:
[1] 1 2 4 5
Maximum Number of Connections: 4 for Regions with
Indices:
[1] 1 2 4 5
Missing Row Indices:
[1] 3
Missing Column Indices:
[1] 4
Indices of Regions with No Connections (islands):
[1] "none"
This shows that location 3 has no neighbors yet it is the
neighbor for
some locations and location 4 is not the neighbor for any
location
but it has neighbors.
Hope this explanation helps.
Stephen Kaluzny (spk@insightful.com)
Insightful Corp.
---------------------------------------------------
and Nick Ellis wrote:
An 'island' in the sense of the check.islands function is a
single point
without neighbours, so that, by that definition your example
would have only
1 island.
Below is a way of finding the islands (let's call them
continents) in your
example.
[....]
==========================================================================
x <- matrix(1,ncol=5,nrow=4)
x[3,] <- x[,4] <- 0 # matrix with 4 'continents'
x
xpad <- matrix(0,6,7) # padded version (for indexing below)
xpad[2:5,2:6] <- x
i <- matrix(1:20,4,5) # index for each of the 4x5 points
ipad <- matrix(0,6,7) # padded version
ipad[2:5,2:6] <- i
ipad
right<-ncol(xpad);left<-1;top<-1;bottom<-nrow(xpad)
# find all neighbours: a neighbour is point to immediate N S
E W NE SE SW or
NW
nhbr <- rbind(
cbind(ipad[,-right][(xpad[,-right]==1)
&(xpad[,-right]==xpad[,-left])],ipad[,-left][(xpad[,-left]==1)
&(xpad[,-right]==xpad[,-left])]),
cbind(ipad[-bottom,][(xpad[-bottom,]==1)
&(xpad[-bottom,]==xpad[-top,])],ipad[-top,][(xpad[-top,]==1)
&(xpad[-bottom,]==xpad[-top,])]),
cbind(ipad[-bottom,-right][(xpad[-bottom,-right]==1)
&(xpad[-bottom,-right]==xpad[-top,-left])],ipad[-top,-left][(xpad[-top,-left
]==1) &(xpad[-bottom,-right]==xpad[-top,-left])]),
cbind(ipad[-bottom,-left][(xpad[-bottom,-left]==1)
&(xpad[-bottom,-left]==xpad[-top,-right])],ipad[-top,-right][(xpad[-top,-rig
ht]==1) &(xpad[-bottom,-left]==xpad[-top,-right])]),
)
nhbr <- rbind(nhbr,nhbr[,2:1]) # if a neighbours b then b
neighbours a
nhbr <- nhbr[order(nhbr[,1]),]
nhbr
xbig <- matrix(0,20,20) # 20 x 20 neighbourhood matrix for
the 20 points in
the 4x5 matrix x
xbig[nhbr] <- 1
diag(xbig) <- 1 # a neighbours a
# Note this means the following 'sea' points are neighbours
of themselves
3,7,11,13:16,19
xprod <- xbig # xprod will find all points in same
'continent'
prev <- xprod*0
while (any(prev!=xprod)) {
prev <- xprod
xprod <- (xprod%*%xbig > 0)+0
}
xprod <-
xprod[!duplicated(apply(xprod,1,paste,collapse=".") ),] #
xprod
has unique rows showing constituents of each continent
lapply(1:nrow(xprod),function(i) (1:20)[xprod[i,]>0]) # a
list of the
'continents' (contiguous land masses, includes 'sea')
sea <- c(3,7,11,13:16,19)
lst <- lapply(1:nrow(xprod),function(i)
(1:20)[-sea][xprod[i,-sea]>0])
lst[sapply(lst,length)>0] # a list of the 'continents'
(contiguous land
masses, excludes 'sea')
# The number of contintents is the length(lst)
==========================================================================
Nick Ellis
CSIRO Marine Research mailto:Nick.Ellis@marine.csiro.au
PO Box 120 ph +61 (07) 3826 7260
Cleveland QLD 4163 fax +61 (07) 3826 7222
Australia http://www.marine.csiro.au
Again, many thanks Nick.
Cheers,
Mrigesh
--
----------------------------------------------------------
Mrigesh Kshatriya
Department of Zoology & Entomology,
University of Pretoria, Pretoria, 0002, South Africa
voice: (27-12) 420-4282 fax: (27-12) 420-3210
home: (27-12) 420-3677
email: mkshatriya@zoology.up.ac.za
----------------------------------------------------------
|