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

  download.file("https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2014/01/Zone9_2014_01.zip", "temp_file.zip")
  unzip("temp_file.zip", exdir="temp")
  broadcast <- readOGR("temp/Zone9_2014_01.gdb", layer = "Zone9_2014_01_Broadcast")

  # Remove unneeded fields
  broadcast@data <- broadcast@data[,c("SOG", "BaseDateTime", "Status", "VoyageID", "MMSI")]

  # 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
  coords <- coordinates(broadcast)[,2]

  chnks <- list()
  chnks[[1]] <- broadcast[which(coords<=37),]
  chnks[[2]] <- broadcast[which(coords>37 & coords<=40),]
  chnks[[3]] <- broadcast[which(coords>40 & coords<=44),]
  chnks[[4]] <- broadcast[which(coords>44),]

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.

  • 2
    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(Polygon(cbind(c(-117, -122, -129, -129, -124, -124, -117), 
                       c(32, 30, 39, 47, 49, 39, 32)), hole = FALSE)), "EEZ")), 
  proj4string = CRS(proj4string(broadcast))

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

@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
  • 1
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.