Subset SpatialLinesDataFrame with gIntersection and reattach data to new SpatialLinesDataFrame

Dear list,
I have a shapefile (the publically available AQUASTAT Rivers of Africa), which I read into a SpatialLinesDataFrame. I would like to crop this to the extent of a bounding box as I only want to focus on a portion of the data. I would also like to reattach the original data from the data.frame after the cropping has been done (by preserving the IDs of the SpatialLines).

In response to a previous question I had on this list (http://r-sig-geo.2731867.n2.nabble.com/Crop-a-SpatialPolygonsDataFrame-by-a-bounding-box-anotehr-polygon-td7580278.html), Roger Bivand directed me to some lovely code here: https://stat.ethz.ch/pipermail/r-sig-geo/2012-June/015337.html which achieves the same end result for SpatialPolyogns, but I can't seem to get it to work for SpatialLines.


Here is a reproducible example. If you would like to try this out please feel free to run the following which will download the shapefile from my Dropbox. Be warned that the gIntersection command can take a while as it has to run through ~69,000 geometries in the rivers shapefile.

#             Download and unzip the AQUASTAT Rivers of Africa dataset
oldwd <- getwd()
tmp <- tempdir()
setwd(tmp)
url <- "http://dl.dropbox.com/s/pvpi8y9550prxyd/Data.zip"
f <- file.path(tmp,"Data.zip")
download.file(url,f) # File is ~8Mb
unzip(f)

#             Read data into SpatialLinesDataFrame
rivers <- readShapeLines(paste(tmp , "/rivers_africa" , sep = "" ) , proj4string = CRS("+proj=latlong +datum=WGS84") )

#             Create polygon for extent of study area
xmin <- -18;  xmax <- 3; ymin <- 4; ymax <- 15     # Coordinates for bounding box
bb <- cbind(x=c(xmin, xmin, xmax, xmax, xmin), y=c(ymin, ymax, ymax, ymin, ymin))    #Create bounding box
SP <- SpatialPolygons(list(Polygons(list(Polygon(bb)), "1")), proj4string=CRS(proj4string(rivers)))              # Create polygon from bounding box

#             This is Rogers code to first determine for each SpatialLine if it is within the bounding box
gI <- gIntersects(rivers, SP, byid=TRUE)
out <- vector(mode="list", length=length(which(gI)))
ii <- 1

#             And then create a list of only those SpatialLines for which gI = TRUE
for (i in seq(along=gI)) if (gI[i]) {out[[ii]] <-
gIntersection(rivers[i,], SP); row.names(out[[ii]]) <-
row.names(rivers)[i]; ii <- ii+1}


However I t hen run into an error I can't solve. I now want to take this 
list of objects of SpatialLines and return a single SpatialLines object. 
However the rbind method which works for SpatialPolygons does not work 
for SpatialLines (in this instance at least).

table(sapply(out, class))

ri <- do.call("rbind", out)
Error in rbind(<S4 object of class "SpatialLines">, <S4 object of class "SpatialLines">,  :
no method for coercing this S4 class to a vector

methods(rbind)
[1] rbind.data.frame               rbind.fill
[3] rbind.fill.matrix              rbind.gtable*
[5] rbind.SpatialLines             rbind.SpatialLinesDataFrame
[7] rbind.SpatialPixels            rbind.SpatialPixelsDataFrame
[9] rbind.SpatialPoints            rbind.SpatialPointsDataFrame
[11] rbind.SpatialPolygons          rbind.SpatialPolygonsDataFrame

Okay  so maybe try explicitly calling the "rbind.SpatialLines" method:

rsl <- do.call("rbind.SpatialLines", out)
Error in slot(x, "lines") :
 no slot of name "lines" for this object of class "SpatialPoints"

slotNames(out[[1]])
[1] "lines"       "bbox"        "proj4string"

Does anyone know why the function thinks I am giving it SpatialPoints. 
Does anyone know what I am doing wrong or why this wouldn't work? I am 
hoping that If I can get the rbind to work I can remake my 
SpatialLinesDataFrame with the following:

rid <- sapply(slot(rsl, "lines"), function(x) slot(x, "ID"))
small_riv <- SpatialLinesDataFrame(rsl , data = rivers at data[ c(rid) , ] )

Thanks you for any help,

Simon



	[[alternative HTML version deleted]]