Clip SpatialPolygonsDataFrame and preserve data in R

Question

I'm trying to clip a SpatialPolygonsDataFrame and keep it as a SpatialPolygonsDataFrame. However, when I use gIntersect, it converts it to a SpatialPolygons and I lose all the data. I'm pretty new at working with spatial data in R, so I'm probably missing something simple.

As an exercise, I'm trying to cut off the outlying islands of Hawaii.

library(tigris)
st<-tigris::states(cb=T)
hawaii<-subset(st,STATEFP=="15")
class(hawaii)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

hawaii2<-gIntersection(as(extent(-178.3347+18, -154.8068, 18.91036, 28.40212), "SpatialPolygons"), hawaii, byid = TRUE, drop_lower_td = T)
class(hawaii2)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"

Answer 1

You can just coerce the data back into a SpatialPolygonsDataFrame and add the attributes back to the data. You just have to index the original data.frame rows so you know what to relate back to the data.

One issue with your example data is that it is a multipart geometry (multiple polygons collapsed into a single feature and attribute). I am going to explode the geometry so that each polygon is represented individually. The example will still work if this step is not done but it provides a better example in showing how attribute tracking is happening.

First, add packages and data.

library(tigris)
library(sp)
library(rgeos)
library(raster)
library(sf)

st <- tigris::states(cb=T)
hawaii <- subset(st,STATEFP=="15")
  dim(hawaii)

You will see that the dimensions currently indicate that there is only one feature with nine attribute columns. Here, we are going to explode the geometry so that we have the actual number of polygons (n=28). We use sf to do this but coerce back to an sp object at the end.

( x <- as(hawaii, "sf") )
( x <- sf::st_cast(x, substring("MULTIPOLYGON", 6, last = 1000000L)) )
hawaii <- as(x, "Spatial")
  dim(hawaii)

Now we can go back to the example at hand by creating a bounding polygon and then clipping to it. Note that I am using the id = rownames(hawaii@data) argument in gIntersection to ensure that the polygon ids match the rownames in the sp.df data.frame object that we set aside. Looking at the names of the resulting hawaii SpatialPolygons object will show that they match the original rownames in the @data slot of the original data and those of sp.df.

e <- as(extent(-178.3347+18, -154.8068, 18.91036, 28.40212), 
        "SpatialPolygons")
      proj4string(e) <- proj4string(hawaii) 
    sp.df <- hawaii@data   
hawaii <- gIntersection(e, hawaii, byid = TRUE, 
                        id = rownames(hawaii@data), 
                        drop_lower_td = TRUE)
      names(hawaii)

We can now coerce the data back to a SpatialPolygonsDataFrame and use the sp.df data.frame in the data argument. We can use which and %in% to match rownames(sp.df) and names(hawaii). Looking at the object dimensions again we see that there are now 11 features (polygons) with 9 attribute columns.

hawaii <- SpatialPolygonsDataFrame(hawaii, 
    data = sp.df[which(rownames(sp.df) %in% names(hawaii)),])
dim(hawaii)  

Answer 2

Just to offer you an alternative, you could use raster::intersect() function to accomplish the same result. From its help file:

If x is a Spatial* object, a new Spatial* object is returned. If x or y has a data.frame, these are also returned (after merging if necessary) as part of a Spatial*DataFrame, and this is how intersect is different from rgeos::gIntersection on which it depends.

Tweaking a little bit your code:

library(tigris)

#Defining general map projection
crs <- "+proj=longlat +datum=WGS84 +no_defs"

st <- tigris::states(cb=T)
hawaii <- subset(st, STATEFP=="15")

#Assigning general projection to object 'hawaii'
raster::projection(hawaii) <- crs

#Creating your zomming extent and assigning it the general projection so intersect does not issues warnings
window <- as(spatstat::as.extent(c(-178.3347+18, -154.8068, 18.91036, 28.40212)), "SpatialPolygons") 
raster::projection(window) <- crs

#Intersection of both objects
hawaii2<-raster::intersect(hawaii, window)

The object hawaii is:

> str(hawaii@data)
'data.frame':   1 obs. of  9 variables:
 $ STATEFP : chr "15"
 $ STATENS : chr "01779782"
 $ AFFGEOID: chr "0400000US15"
 $ GEOID   : chr "15"
 $ STUSPS  : chr "HI"
 $ NAME    : chr "Hawaii"
 $ LSAD    : chr "00"
 $ ALAND   : chr "16633990195"
 $ AWATER  : chr "11777809026"

The object hawaii2 is now:

> class(hawaii2)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

> str(hawaii2@data)
'data.frame':   1 obs. of  9 variables:
 $ STATEFP : chr "15"
 $ STATENS : chr "01779782"
 $ AFFGEOID: chr "0400000US15"
 $ GEOID   : chr "15"
 $ STUSPS  : chr "HI"
 $ NAME    : chr "Hawaii"
 $ LSAD    : chr "00"
 $ ALAND   : chr "16633990195"
 $ AWATER  : chr "11777809026"

Hope this helps