[R-sig-Geo] speed problem with %over% function

[Colin Rundel]

For this example at least you are severely handicapping rgeos as gIntersect needs to create all new points that result from the intersections which then need to be translated back to R from the geos representation. If you use something like the following you should see much better results:

> system.time(grille.plac4 <- gIntersects(shape,grille.plac,byid=TRUE))
   user  system elapsed 
  0.365   0.008   0.371 
> plot(grille.plac[c(grille.plac4),])

In this case I am using a spatial predicate instead of a topology operator, also note that I've switched the order of the polygon and the points since some of the predicates can make use of prepared geometries to improve performance and rgeos assumes the first geometry should be the prepared geometry (think in terms of the fact that you are comparing a bunch of points to one polygon so prepare the polygon (create a spatial index) to improve performance).  

-Colin

-----

Colin Rundel
Visiting Instructor
Duke University, Department of Statistical Science
colin.rundel at stat.duke.edu
www.stat.duke.edu/~cr173/

On Feb 29, 2012, at 11:19 AM, <Bastien.Ferland-Raymond at mrnf.gouv.qc.ca> <Bastien.Ferland-Raymond at mrnf.gouv.qc.ca> wrote:

> Thanks to Francis and Edzer for their responses.
> 
> I been working on my problem for the last 2 days and I managed to get some very interesting information.  I also prepared a reproducible example so people can try it.
> 
> I've found 3 functions that do the thing I want:
> 1) %over%
> 2) overlay()  (which is deprecated)
> 3) gIntersection()  (from the rgeos package)
> 
> So I've tried them all on my data to compare the speed between each of them and ArcGIS.  The punchline of my analysis is that %over% is the slowest (and by a lot) of all functions including ArcGIS.  overlay() is the fastest.
> 
> To get these results, here is a reproducible example:
> 
> ### first, download the shapes from:
> ### http://www.nceas.ucsb.edu/files/scicomp/demo/read-write-shapefiles.zip
> ### This is not my files, I took them from the example at:
> ### http://www.nceas.ucsb.edu/scicomp/usecases/ReadWriteESRIShapeFiles
> 
> ###  Use the rgdal library to read the files
>> library(rgdal)
>> shape <- readOGR("C:\\Documents and Settings\\ferba1\\Bureau\\shape","nw-counties")
> OGR data source with driver: ESRI Shapefile
> Source: "C:\Documents and Settings\ferba1\Bureau\shape", layer: "nw-counties"
> with 208 features and 8 fields
> Feature type: wkbPolygon with 2 dimensions
> 
> # choose only one territory to make the analysis easier
>> shape <- shape[shape at data$COUNTYNAME=="Lemhi",]
>> plot(shape, axes=T)
> 
>> etendu <- bbox(shape)
> 
> # this code is used to transfer km in degre-decimals. It's not necessary to
> # understand it for this problem.  The only important part is that you can
> # change "separation" to increase the density of the grid. "separation"
> # represent the distance between the grid point, in that example 0.5km.
>> separation <- 0.5
>> circ.minlat <- 2*pi*6378.16*cos(pi*etendu[2,1]/180)
>> pas.degre.lon <- circ.minlat/360
>> nb.degre.lon <- separation/pas.degre.lon
>> pas.degre.lat <- 111.3199
>> nb.degre.lat <- separation/pas.degre.lat
> ###
> 
> # make the grid
>> grille.plac <- GridTopology(etendu[,1],
> +                             cellsize=c(nb.degre.lon,nb.degre.lat),
> +                             cells.dim=(ceiling((etendu[,2]-etendu[,1])/c(nb.degre.lon,nb.degre.lat))+2))
> 
>> grille.plac <- coordinates(grille.plac)
> 
>> grille.plac <- SpatialPointsDataFrame( grille.plac,data.frame(plac_ID=1:nrow(grille.plac)), proj4string = CRS(proj4string(shape)))
> 
> 
>> points(grille.plac,col="blue", pch=".")
> 
> ### The overlay() function
>> system.time(plac.dans.contour1 <- overlay(grille.plac, shape) )
>   user  system elapsed
>   2.33    0.00    2.34
>> grille.plac1 <- grille.plac[!is.na(plac.dans.contour1),]
>> points(grille.plac1,col="red", pch=".")
> 
> ###  The %over% function
>> system.time(plac.dans.contour2<-grille.plac%over%shape)
>   user  system elapsed
> 953.72    2.33  960.80
>> grille.plac2 <- grille.plac[!is.na(plac.dans.contour2)[,1],]
>> points(grille.plac2,col="yellow", pch=".")
> 
> ###  The gIntersection() function
> library(rgeos)
>> system.time(grille.plac3 <- gIntersection(grille.plac, shape))
>   user  system elapsed
>  37.83    0.05   37.99
>> points(grille.plac3,col="green", pch=".")
> 
> ###  end of code
> 
> As you can see, the %over% function is the slowest with 961 seconds to run.  The overlay() function was the best with 2.34 seconds to run.  The gIntersection was in the middle with 38 seconds to run.  Finally, the intersection tool of ArcGIS took 20 seconds to do the same thing.  Those times seem to be context specific, as in my problem, gIntersection() was the slowest with 799 seconds, while %over% did the task in 250 seconds and overlay() in 119.  overlay() was still the fastest function.
> 
> So, it seems there clearly is room for amelioration of the %over% function.  And my next question is, why is overlay() deprecated? In the help, it says: "This function is deprecated, as it has some inconsistences." but I did not find any reference on what are those "inconsistences".  Going from overlay() to %over% seems more like a downgrade than a upgrade to me.  I implemented overlay() in my code instead of %over%, saving an incredible amount of time. However I'm a little stressed about it as it may cause problems or I may lose the function in a future version of the sp package.
> 
> So I guess anybody motivated to ameliorate the %over% function could start by understanding the difference between it and the overlay() function.  Sadly, I'm no programmer so I can' help here.
> 
> Best regards,
> Bastien Ferland-Raymond
> 
> 
> Date: Tue, 28 Feb 2012 11:38:04 -0500
> From: Francis Markham <francis.markham at anu.edu.au>
> 
> I've had examples in the past where using %over% from the sp package takes
> all available RAM (7GB) and several hours, while ArcGIS takes about 300MB
> and 5 minutes, so I would agree that there is plenty of room for
> improvement here. I'll try to give a reproducable example in the coming
> weeks when I return from travel.
> 
> This is a critical issue for me insofar as spatial joins are a routine
> procedure for me and but without reasonable speed for producedures such as
> this I cannot perform all my analysis in R, making for a more complex,
> error prone work flow and scuttling the possibility of "reproducable
> research."
> 
> On a related note, does the 'sp' package have an accessible bug tracker?
> I'd like to be able to contribute to improving this very useful package if
> at all possible, but I'm not sure where to begin.
> 
> Warm regards,
> 
> Francis Markham
> Research Associate
> Fenner School of Environment and Society
> Australian National University