R equivalent of ArcMap split by Attributes?

Question

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",
           overwrite_layer=TRUE)
}

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?

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)

Second,

  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",
           overwrite_layer=TRUE)
  }
}

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)

CFALSE  CTRUE 
     1      9 

Oh well seems its all working.