rgeos – Update

So I have been remiss in posting updates of progress on rgeos as my summer has gotten busier but the project is progressing smoothly and everything is on schedule. We have just recently passed an important personal milestone, parsing and running of all the JTS unit tests are now working and the rgeos code passes all but a handful of tests. (Failing tests appear mostly to do with issues surrounding dimension collapse which, I do not believe, is currently handled by geos). We are nearly feature complete and will now be focusing on testing and writing up documentation / examples.

Below are a couple of examples that came up as answers to questions posted to the R-sig-geo mailing list. Note that the naming convention of the geos functions has changed to use a single g (for geos) as a prefix instead of using the awkward RGEOS per Hadley’s earlier suggestion. Old function names are retained for the time being but will throw a depricated warning.

Pruning a Delaunay tessalation using a polygon:

library(rgeos)
library(spdep)

#bounding polygon
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"))) 

# points that form the basis of the Delaunay tessalation
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.nb = tri2nb(dat) 


filternb = function(nb, coords, bound) { 
        sl = nb2lines(nb,as.list(rep(NA,length(nb))), coords) 
        
        # any line that crosses the polygon should be removed
        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) 

par(mfrow=c(1,2))

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

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

IDing polygons containing points:

This can be done with existing functions in sp (see the overlay function) but this is how you could do it with rgeos. (I just realized that this is pretty much exactly the same as the cities example from my last post)

library(rgeos) 
gt <- GridTopology(c(0.05,0.05), c(.1,.1), c(10,10)) 
grd <- SpatialGrid(gt) 
spi <- as(grd, "SpatialPixels") 
spol <- as(spi, "SpatialPolygons") 

set.seed(1) 
x=runif(25) 
y=runif(25) 

pts = SpatialPoints(cbind(x,y)) 

ct = gContains( spol,pts,byid=c(TRUE,TRUE) ) 
colsub = apply(ct, 2,any) 

plot(spol) 
plot(spol[colsub,],add=T,col='green') 
plot(pts,add=T,pch=16)