[R-sig-Geo] Subset SpatialLinesDataFrame with gIntersection and reattach data to new SpatialLinesDataFrame

[Roger Bivand]

On Sat, 14 Jul 2012, O'Hanlon, Simon J wrote:

> 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))

  SpatialLines SpatialPoints
          3303             1

shows that one returned object is non-conformant, so:

out_class <- sapply(out, class)
ri <- do.call("rbind", out[out_class == "SpatialLines"])

works.

Hope this helps,

Roger

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

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no