Extracting the parent polygon from a nested SpatialPolygonsDataFrame or 'Dissolving' holes from a parent ploygon

Question

EDIT After more research and still no solution, I am adding a substantial edit as well as a link to the .shp file.

The shape file is included here

I have a SpatialPolygonsDataFrame that contains 9 polygons, each of which also contains multiple nested polygons - the 'holes'. The summary of the data are here.

> summary(data)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min       max
x  483298.9  643204.4
y 4782172.1 4997248.3
Is projected: TRUE 
proj4string :
[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
       Id          IndID  
 Min.   :0   BHS_011_A:1  
 1st Qu.:0   BHS_015_A:1  
 Median :0   BHS_083_A:1  
 Mean   :0   BHS_089_A:1  
 3rd Qu.:0   BHS_091_A:1  
 Max.   :0   BHS_129_A:1  
             (Other)  :3  

A sample of the structure of the data is below.

Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 9 obs. of  2 variables:
  .. ..$ Id   : int [1:9] 0 0 0 0 0 0 0 0 0
  .. ..$ IndID: Factor w/ 9 levels "BHS_011_A","BHS_015_A",..: 1 2 3 4 5 6 7 8 9
  ..@ polygons   :List of 9
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 5
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 513497 4986246
  .. .. .. .. .. .. ..@ area   : num 76614017
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:287, 1:2] 509244 507384 507214 507010 506899 ...
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 509678 4979511
  .. .. .. .. .. .. ..@ area   : num 1462398
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 509301 509269 509194 509007 509412 ...
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 515572 4988493
  .. .. .. .. .. .. ..@ area   : num 1579348
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:10, 1:2] 514520 514570 514684 516501 515996 ...
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"

As seen in the example below (one of nine) the parent Polygon has multiple holes.

data <- readOGR(".", "IndLineBuff")
plot(data[data$IndID == "MTG_005_A",])

This is my first foray into the sp(), rgdal(), rgeos(), and other spatial packages and I have found a number of helpful posts regarding the use of operators to extract area, etc, but nonetheless, questions remain. While this post offers solutions that are close, I cannot seem to adapt the code to fit the needs described herein.

I want to obtain a SpatialPolygonsDataFrame that contains only the parent (biggest) polygon from each set of sublists (i.e. data@polygon). It seems like I should be able to either extract only the parent polygon or 'dissolve' the holes.

The final out come would be 9 polygons, each the parent of the 9 lists, that I could export as an ESRI shapefile.

Any suggestions would be appreciated.

Answer 1

The following hack might work:

data@polygons = lapply(data@polygons, function(x) { x@Polygons = x@Polygons[1]; x})
plot(data[data$IndID == "MTG_005_A",])

In your case, it seems you've buffered linear features or so, and want to get rid of the holes in polygons resulting. Each of the 9 polygons seems to consist of the main polygon followed by a set of wholes. The function selects the first of this list, returning the list-of-lists required. To understand the class structure, study ?"SpatialPolygons-class", ?"Polygons-class" and ?"Polygon-class" carefully.

Answer 2

When I was first learning about sp and related packages I found Barry Rowlingson’s site very useful.

It includes:

• A cheatsheet
  Really useful reference sheet with hints about how to do things using R packages\ http://www.maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/cheatsheet.html

• Introduction to Spatial
  An explanation of the spatialDataframe structures created by package sp. Also explains CRS (coordinate reference systems) and projections (the proj4string and proj4.Project). Finally, also covers use and manipulation of raster images, and extracting data from them using R\ http://www.maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/static/talk1.pdf

Getting to your question, it is messy because you are dealing with lists of lists of things, but it can be done.

The following code seems to work (I've tested it (a bit) on a shapefile I have already, rather than downloaded yours). If you replace NSWACTA with the name of your shapefile after you have read it in, it should work for your case. No promises though :-).

The first part is introductory, to show what's going on with all the lists of lists of classes. Hope this helps you figure out what's going on. It is the dropHole function and it's methods which actually do what you want.

class(NSWACTA)
# NSWACTA has class sp::SpatialPolygonsDataFrame
slotNames(NSWACTA)
# The SpatialPolygonsDataFrame has a slot called @polygons
class(NSWACTA@polygons)
# The @polgyons slot contains a list
class(NSWACTA@polygons[[1]])
# Elements of that list are of class sp::Polygons 
slotNames(NSWACTA@polygons[[1]])
# An sp::Polygons class has several slots, one of which is @Polygons
class(NSWACTA@polygons[[1]]@Polygons)
# The @Polygons slot contains a list
class(NSWACTA@polygons[[1]]@Polygons[[1]])
# The entries in the @Polygons list are of class sp::Polygon
slotNames(NSWACTA@polygons[[1]]@Polygons[[1]])
# An sp::Polygon has several slots. 
# @coords holds the coordinates which define the boundary.
# @hole is a logical field indicating whether or not the sp::Polygon is a hole 

setGeneric('dropHole',def=function(poly, ...){
  standardGeneric('dropHole')
})

# Drop a single sp::Polygon if it is holey
setMethod('dropHole',signature=signature('Polygon'),
          def=function(poly) {
            #return only Polygons which are not holes
            if (poly@hole) NULL else poly
          }
          )

# Drop holey sp::Polygon entries in the @Polygons list of an sp::Polygons class
setMethod('dropHole', signature = signature('Polygons'),
          def = function(poly) {
            noHoles <- lapply(poly@Polygons, dropHole)
            #Remove the NULL entries from the list
            noHoles <- Filter(Negate(is.null), noHoles)
            # Turn back into a (single) Polygons
            # The generator function (sp::Polygons) fills in the other slots!
            # return the new sp:Polygons object
            sp::Polygons(noHoles, ID = poly@ID)
          }
          )

# Drop holey parts of sp::Polygons in the @polygons list 
# of an sp::SpatialPolygonsDataFrame
setMethod('dropHole', signature = signature('SpatialPolygonsDataFrame'),
          def = function(poly) {
            noHoles <- lapply(poly@polygons, dropHole)
            # Put the un holey Polygons list back into the @polygons slot 
            poly@polygons <- noHoles
            #return the modified SpatialPolygonsDataFrame 
            poly
          }
          )

newmap <- dropHole(NSWACTA)

Here's the output from the first part. Take care to note the differences between upper and lower case in names, and between singular and plurals! (ie polygons, Polygons, and Polygon)

d> class(NSWACTA)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
d> # NSWACTA has class sp::SpatialPolygonsDataFrame
d> slotNames(NSWACTA)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
d> # The SpatialPolygonsDataFrame has a slot called @polygons
d> class(NSWACTA@polygons)
[1] "list"
d> # The @polgyons slot contains a list
d> class(NSWACTA@polygons[[1]])
[1] "Polygons"
attr(,"package")
[1] "sp"
d> # Elements of that list are of class sp::Polygons 
d> slotNames(NSWACTA@polygons[[1]])
[1] "Polygons"  "plotOrder" "labpt"     "ID"        "area"     
d> # An sp::Polygons class has several slots, one of which is @Polygons
d> class(NSWACTA@polygons[[1]]@Polygons)
[1] "list"
d> # The @Polygons slot contains a list
d> class(NSWACTA@polygons[[1]]@Polygons[[1]])
[1] "Polygon"
attr(,"package")
[1] "sp"
d> # The entries in the @Polygons list are of class sp::Polygon
d> slotNames(NSWACTA@polygons[[1]]@Polygons[[1]])
[1] "labpt"   "area"    "hole"    "ringDir" "coords" 
d> # An sp::Polygon has several slots. 
d> # @coords holds the coordinates which define the boundary.
d> # @hole is a logical field indicating whether or not the sp::Polygon is a hole