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

[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>:
>
> 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
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>