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 <[email protected]>: > > 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. > > _______________________________________________ > R-sig-Geo mailing list > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
