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

Question

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)

  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?

Addition

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.

Answer 1

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")), 
  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
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.