Using 'sp' method to subset SpatialPointsDataFrame with SpatialPolygon returns error in parLapply()

I have numerous very large data sets of ship tracking data (Automated Information System; AIS). These are originally in a file geodatabase and can be downloaded here (this is a smaller file for testing purposes). After reading in the file using readOGR, I subset the data by attributes and then into latitudinal bands like this (warning - file is large and takes a while to download and read in. You will also need ReadFGDB driver support in your version of rgdal):

  # Libraries
library(rgdal)
library(sp)
library(parallel)

unzip("temp_file.zip", exdir="temp")

# Remove unneeded fields

# Subset the data
broadcast <- broadcast[broadcast@data$SOG>0.5,] # remove all points w speed <0.5 broadcast <- broadcast[broadcast@data$Status==0 | broadcast@data$Status==7 | broadcast@data$Status==8,] # remove all points at anchor

# Split data set into 4 regions

chnks <- list()


Now chnks is a list with four SPtsDFs which I want to clip with a polygon, EEZ. (EEZ is the US West Coast Exclusive Economic zone). I'll substitute a simple polygon for EEZ here:

EEZ <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(-117, -122, -129, -129, -124, -124, -117),c(32, 30, 39, 47, 49, 39, 32)), hole=F)), "EEZ")))


If I use lapply to do spatial subset on chnks, it works fine:

chunked <- lapply(chnks, FUN = function(x){
# Clip the AIS points to EEZ
return(x[EEZ, ])
})


But if I try to use parLapply, it fails with an error:

  # create cluster
cluster <- makeCluster(4, outfile = "")
clusterExport(cl=cluster, varlist = "EEZ")

chunked <- parLapply(cl = cluster, chnks, fun = function(x){
# Clip the AIS points to EEZ
return(x[EEZ, ])
})


4 nodes produced errors; first error: object of type 'S4' is not subsettable

It seems like the sp method for [ (which references points[!is.na(over(points, geometry(polygons))),] according to this post) is not being loaded to the workers.

Does anyone have an idea how to avoid this error or explicitly export the sp method to the parallel clusters?

I also tried using plyr and llply() with .parallel = T and got the same failure message. I would also note that the sp package seems to load fine on all of the workers, which I would think should allow the [ subsetting method for spatial objects.

• Have you tried loading package methods on each of the nodes? – Edzer Pebesma Apr 12 '16 at 5:42

The problem you are experiencing is actually twofold.

• Firstly, and as pointed out by @EdzerPebesma, you are not loading the required packages on each of the 4 nodes separately. You have to use clusterEvalQ to tell each node which packages it is going to need to fulfill the required (spatial) tasks.

• Secondly, you need to assign a proj4string to the polygon you use for clipping.

Otherwise, you will sooner or later run into

Error: is(proj4string, "CRS") is not TRUE

even though the parallel cluster is now setup properly. Here is some sample code that should meet your demands.

## assign a crs to 'EEZ'
EEZ <- SpatialPolygons(
list(Polygons(
list(Polygon(cbind(c(-117, -122, -129, -129, -124, -124, -117),
c(32, 30, 39, 47, 49, 39, 32)), hole = FALSE)), "EEZ")),
)

### parallelization -----

## register parallel backend
cl <- makePSOCKcluster(detectCores() - 1)

## load required packages ('rgdal' automatically attaches 'sp')
clusterEvalQ(cl, library(rgdal))

## export required objects
clusterExport(cl, "EEZ")

### perform operation -----

chunked <- parLapply(cl, chnks, function(x) {
x[EEZ, ]
})

## deregister parallel backend
stopCluster(cl)


@TimSalabim once wrote a nice tutorial about reading raster data in parallel. It comes very close to the above code snippet, so maybe you want to have a look at it.

• Thanks @fdetsch! That was very helpful. I am doing this work on Ubuntu and had missed that I needed to use clusterEvalQ() in conjunction with parLapply() to load packages. It still seems odd that sp looked to be loading on the workers, but not the methods functionality. Is that because the sp methods are stored in package methods as @Edzer Pebesma mentioned in the comment above, not in sp? – Cotton.Rockwood Apr 12 '16 at 23:13
• Also, do I need to export chnks to the workers even though it is the list being passed to parLapply? It would seem counterproductive to have to send the whole list object to all workers when they only need the element they will process. – Cotton.Rockwood Apr 13 '16 at 15:31
• I confirmed that I do not need to export chnks since it is passed to parLapply. Doing so only takes up memory on the workers unnecessarily. The solution to my issue was explicitly loading the package rgdal (and dependency sp) to the workers with clusterEvalQ(cl, library(rgdal)). Thanks! – Cotton.Rockwood Apr 13 '16 at 18:37
• A classic win-win situation. I didn't know that exporting chnks is not necessary, thanks for that! – fdetsch Apr 14 '16 at 6:24