[R-sig-Geo] Pruning a delaunay triangulation with a shape

[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
 
-- 
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Pruning-a-delaunay-triangulation-with-a-shape-tp5324088p5327334.html
Sent from the R-sig-geo mailing list archive at Nabble.com.