I have a shapefile of 50 polygons. I've read it into R. Now I want to use R to export 50 shapefiles, one for each polygon. This seems like an intuitive application of apply. But I'm having a devil of a time doing it.

I tried a for loop (like Splitting vector layer in R).

STNCODES <- unique(ptsbuf_400@data$STNCODE)
for (i in 1:length(STNCODES)) {
  tmp <- data[ptsbuf_400@data$STNCODE == STNCODES[i], ]
  writeOGR(tmp, dsn=getwd(), STNCODES[i], driver="ESRI Shapefile",

I get the following error:

Error in data[ptsbuf_400@data$STNCODE == STNCODES[i], ] : object of type 'closure' is not subsettable

This seems to be the cause of the error, but I'm at a loss how to fix it:

data[ptsbuf_400@data$STNCODE == STNCODES[i], ]

Amtrak station data here

Example using SF works--likely there is some difference between the spatial object generated by readOGR/writeOGR and read_sf and write_sf--the latter seems to be a vector, while the sp object is not?

  • If you think this is a different question, and are stuck on your code, then just include a code snippet that illustrates what you have tried and where you are stuck.
    – PolyGeo
    Jan 19, 2020 at 21:41

1 Answer 1


First, you don't really need to @data to get columns in an sp object. Its not the root of your problem but it makes things neater:

STNCODES <- unique(ptsbuf_400$STNCODE)


  tmp <- data[ptsbuf_400@data$STNCODE == STNCODES[i], ]

should make a subset of ptsbuf, not data:

  tmp <- ptsbuf[ptsbuf_400$STNCODE == STNCODES[i], ]

data is an R function, so the error is telling you that you are trying to subset a function (an object of type 'closure') instead of the spatial object you want.

and for good measure here's how to put that in a function:

splitByAttributes = function(spdata, codename, dsn=getwd()){
CODES <- unique(spdata[[codename]])
for (i in 1:length(CODES)) {
  tmp <- spdata[spdata[[codename]] == CODES[i], ]
  writeOGR(tmp, dsn=dsn, CODES[i], driver="ESRI Shapefile",

which makes it easier to test and reuse.

> d = data.frame(x=runif(10),y=runif(10),g=paste0("C",runif(10)>.5))
> coordinates(d)=~x+y
> head(d)
              coordinates     g
1  (0.9877652, 0.1615483) CTRUE
2 (0.3598359, 0.06623409) CTRUE
3  (0.5372637, 0.0559466) CTRUE
4   (0.696683, 0.6759665) CTRUE
5  (0.2361368, 0.1350742) CTRUE
6  (0.9232246, 0.8288764) CTRUE
> splitByAttributes(d, "g")

That produces two shapefiles that I can then test by reading in:

> ctrue = readOGR(".","CTRUE")
OGR data source with driver: ESRI Shapefile 
Source: "/home/rowlings/Downloads/SO/vrt", layer: "CTRUE"
with 9 features
It has 1 fields
> cfalse = readOGR(".","CFALSE")
OGR data source with driver: ESRI Shapefile 
Source: "/home/rowlings/Downloads/SO/vrt", layer: "CFALSE"
with 1 features
It has 1 fields

9 points in one and only 1 in the other? seems unlikely. what was the original grouping?

> table(d$g)

     1      9 

Oh well seems its all working.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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