Skip to content

Pruning a delaunay triangulation with a shape

4 messages · Pierre Roudier, Colin Rundel

Original
Pierre Roudier
# 
Dear researchers,

I'm scratching my head on a problem involving neighbourhood
definition: I am using the Delaunay triangulation to generate a
neighbourhood relationship object on a set of georeferenced points.

I am using the tripack package for that (tri.mesh() for the
triangulation and neighbours() to create the neighbourhood).

However, I would like to prune this neighbourhood using a polygon, so
that some neighbourhood relationships won't be taken into account.
Here's an illustrating example:

dat <- data.frame(
x=c(0.34433527, 0.08592805, 0.55564179, 0.03938242, 0.98044051,
0.19835405, 0.94186612, 0.56208017, 0.31093811, 0.54341230,
0.93508703, 0.38160473, 0.89435383, 0.55457428, 0.22406338),
y=c(0.997803213, 0.094993232, 0.509774611, 0.615526709, 0.004314438,
0.676662461, 0.026060349, 0.165807011, 0.596449068, 0.469553704,
0.888788079, 0.163129754, 0.340335312, 0.621845245, 0.019412254)
)

dat.tr <- tri.mesh(dat)

bound <- data.frame(
x=c(0.34063939, 0.56754974, 0.95361248, 0.96464284, 0.60694389,
0.58173163, 0.58330740, 0.91421832, 1.00403700, 0.96464284,
0.50294332, 0.39263968, 0.22560845, 0.02548613, 0.25397225,
-0.01233226),
y=c(1.02171515, 0.70486571, 0.92207697, 0.84834471, 0.62116963,
0.53747356, 0.47569788, 0.35812482, 0.02134774, -0.07231216,
0.15287015, 0.14489909, -0.03444965, 0.09707276, 0.56138672,
0.58729265)
)

plot(dat)
polygon(bound, col='red')
plot(dat.tr, add=TRUE, col='blue')

In that example, I would like to eliminate all the (blue)
neighbourhood relationships that are not in the red polygon.

Would anybody has a simple way to do this?

Thanks,

Pierre
Colin Rundel
# 
Depending on what you want to do downstream this may or may not be the best
solution for this problem (say if you need to use the tri class at some
later point) but if you just want the neighbors then you could use something
like the following which makes use of rgeos to find which relationships
(lines) are within the polygon. I written a quick function to translate a
tri object into a set of spatial lines, the IDs of each Lines object has the
index of the two points that it connects.

tri.asSpLines = function(x) {

    if (!inherits(x, "tri")) 
        stop("x must be of class \"tri\"")
    tnabor <- integer(x$tlnew)
    nnabs <- integer(x$n)
    nptr <- integer(x$n)
    nptr1 <- integer(x$n)
    nbnos <- integer(x$n)
    ans <- .Fortran("troutq", as.integer(x$nc), as.integer(x$lc), 
        as.integer(x$n), as.double(x$x), as.double(x$y),
as.integer(x$tlist), 
        as.integer(x$tlptr), as.integer(x$tlend), as.integer(6), 
        nnabs = as.integer(nnabs), nptr = as.integer(nptr), nptr1 =
as.integer(nptr1), 
        tnabor = as.integer(tnabor), nbnos = as.integer(nbnos), 
        na = as.integer(0), nb = as.integer(0), nt = as.integer(0), 
        PACKAGE = "tripack")

	LinesList = list()
	k=0
    for (i in 1:x$n) {
        inb <- ans$tnabor[ans$nptr[i]:ans$nptr1[i]]
        for (j in inb) {
        	LinesList[k] = Lines(list(Line(cbind(x=c(x$x[i], x$x[j]),
y=c(x$y[i], x$y[j])))), ID=paste(i,j) )
        	k=k+1
        }
    }
	
	return(SpatialLines(LinesList))
}

dat <- data.frame( 
x=c(0.34433527, 0.08592805, 0.55564179, 0.03938242, 0.98044051, 
0.19835405, 0.94186612, 0.56208017, 0.31093811, 0.54341230, 
0.93508703, 0.38160473, 0.89435383, 0.55457428, 0.22406338), 
y=c(0.997803213, 0.094993232, 0.509774611, 0.615526709, 0.004314438, 
0.676662461, 0.026060349, 0.165807011, 0.596449068, 0.469553704, 
0.888788079, 0.163129754, 0.340335312, 0.621845245, 0.019412254) 
) 

dat.tr <- tri.mesh(dat) 
dat.sl <- tri.asSpLines(dat.tr)

#only change to bound is to add the first point to the end as sp expects
polygons to be closed
bound <- SpatialPolygons(list(Polygons(list(Polygon(cbind( 
x=c(0.34063939, 0.56754974, 0.95361248, 0.96464284, 0.60694389, 
0.58173163, 0.58330740, 0.91421832, 1.00403700, 0.96464284, 
0.50294332, 0.39263968, 0.22560845, 0.02548613, 0.25397225, 
-0.01233226, 0.34063939), 
y=c(1.02171515, 0.70486571, 0.92207697, 0.84834471, 0.62116963, 
0.53747356, 0.47569788, 0.35812482, 0.02134774, -0.07231216, 
0.15287015, 0.14489909, -0.03444965, 0.09707276, 0.56138672, 
0.58729265, 1.02171515) ) )),ID="bound")))

#this is the only rgeos function, it returns true for any geometry that
crosses the bound polygon.
#byid parameter treats each Lines object as a separate geometry instead of a
single multiline geometry
(out=gCrosses(dat.sl,bound,byid=TRUE))

plot(bound,col='red')
plot(dat.sl[!out,],col='blue',add=T)
row.names(dat.sl[!out,]

Hope this helps,
-Colin
Pierre Roudier
# 
Thanks Colin,

Yes, the code is great, that is exactly what I want to do. rgeos will
really be a great tool for us!

However, in the best of the worlds, I would indeed need to get back to
a tri object (or a nb object). Do you know any straightforward way to
do this?

Cheers,

Pierre

2010/7/23 rundel <rundel at gmail.com>:
Colin Rundel
# 
I'm not very familiar with either the tri or nb classes so I'm not entirely
sure what is possible. From a cursory look it seems like this is not
possible with tri as it seems to only internally store the point locations
and recalculates the triangulation whenever the object is printed or
plotted. I have a function below that I think works with nb classes, but I
haven't tested it.



bound <- SpatialPolygons(list(Polygons(list(Polygon(cbind( 
x=c(0.34063939, 0.56754974, 0.95361248, 0.96464284, 0.60694389, 
0.58173163, 0.58330740, 0.91421832, 1.00403700, 0.96464284, 
0.50294332, 0.39263968, 0.22560845, 0.02548613, 0.25397225, 
-0.01233226, 0.34063939), 
y=c(1.02171515, 0.70486571, 0.92207697, 0.84834471, 0.62116963, 
0.53747356, 0.47569788, 0.35812482, 0.02134774, -0.07231216, 
0.15287015, 0.14489909, -0.03444965, 0.09707276, 0.56138672, 
0.58729265, 1.02171515) ) )),ID="bound")))


dat <- data.frame( 
x=c(0.34433527, 0.08592805, 0.55564179, 0.03938242, 0.98044051, 
0.19835405, 0.94186612, 0.56208017, 0.31093811, 0.54341230, 
0.93508703, 0.38160473, 0.89435383, 0.55457428, 0.22406338), 
y=c(0.997803213, 0.094993232, 0.509774611, 0.615526709, 0.004314438, 
0.676662461, 0.026060349, 0.165807011, 0.596449068, 0.469553704, 
0.888788079, 0.163129754, 0.340335312, 0.621845245, 0.019412254) 
) 

dat.tr <- tri.mesh(dat) 
dat.sl <- tri.asSpLines(dat.tr)

dat.nb = tri2nb(dat)


filternb = function(nb, coords, bound) {
	sl = nb2lines(nb,as.list(rep(NA,length(nb))), coords)
	
	out=!gCrosses(sl,bound,byid=TRUE)
	
	k=1
	for(i in 1:length(nb)) {
		l = length(nb[[i]])
		nb[[i]] = nb[[i]][ out[k:(k+l-1)] ]
		k=k+l
	}
	
	return(nb)
}

dat.nb.new = filternb(dat.nb, dat, bound)

plot(bound,col='red')
plot(dat.nb.new,dat,add=T)

-Colin