I am working with a dataset consisting of a list of coordinates for a set of polygons. Some of these polygons are holes. However, the dataset does not follow the convention of listing the coordinates of holes in counter-clockwise direction. Therefore, I have to check for each polygon whether they are a hole, and use the result as argument to the Polygon
function. Here is a simple dummy dataset:
library(sp)
# Coordinates for two dummy polygons (squares), one of which is fully contained within the other
x1 <- c(5,9,9,5,5)
y1 <- c(5,5,1,1,5)
x2 <- c(6,8,8,6,6)
y2 <- c(4,4,2,2,4)
tables <- list() # Will contain the coordinates of each object, each in a data.frame
tables[[1]] <- data.frame(x1,y1)
tables[[2]] <- data.frame(x2,y2)
nobjs <- 2 # Number of objects
Next, I implemented an ugly but functional way of determining whether a polygon is a hole (i.e. fully contained within another polygon), by checking whether all of the nodes of a certain polygon are contained in any other polygon of the dataset:
holes <- logical(nobjs)
# Check for polygons that are holes
for(i in 1:nobjs){
vectTF <- vector("logical",length(tables[[i]]$V1))
for(j in 1:length(tables[[i]][,1])){
vectTFj <- vector("logical",length(tables))
for(k in 1:length(tables)){
if(k==i){
vectTFj[k] <- FALSE
} else {
vectTFj[k] <- point.in.polygon(tables[[i]][j,1],tables[[i]][j,2],tables[[k]][,1],tables[[k]][,2])
}
}
vectTF[j] <- as.logical(sum(vectTFj))
}
holes[i] <- as.logical(sum(vectTF))
}
In this example, the vector holes
correctly tells me that my second object is contained within another polygon, while the first one is not:
holes
[1] FALSE TRUE
Next, I try to build a Polygons object using the vector holes
as the hole
argument to the Polygon
function:
objs <- list()
id <- integer(nobjs)
for(i in 1:nobjs){
objs[[i]] <- Polygons(list(Polygon(tables[[i]][,1:2],hole=holes[i])),as.character(i))
id[i] <- i
}
Now, when I check the structure of the new Polygons
object, I see that none of my polygons are marked as being holes:
str(objs)
List of 2
$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. ..@ Polygons :List of 1
.. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. ..@ labpt : num [1:2] 7 3
.. .. .. .. ..@ area : num 16
.. .. .. .. ..@ hole : logi FALSE
.. .. .. .. ..@ ringDir: int 1
.. .. .. .. ..@ coords : num [1:5, 1:2] 5 9 9 5 5 5 5 1 1 5
.. .. .. .. .. ..- attr(, "dimnames")=List of 2
.. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..$ : chr [1:2] "x1" "y1"
.. ..@ plotOrder: int 1
.. ..@ labpt : num [1:2] 7 3
.. ..@ ID : chr "1"
.. ..@ area : num 16
$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. ..@ Polygons :List of 1
.. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. ..@ labpt : num [1:2] 7 3
.. .. .. .. ..@ area : num 4
.. .. .. .. ..@ hole : logi FALSE
.. .. .. .. ..@ ringDir: int 1
.. .. .. .. ..@ coords : num [1:5, 1:2] 6 8 8 6 6 4 4 2 2 4
.. .. .. .. .. ..- attr(, "dimnames")=List of 2
.. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..$ : chr [1:2] "x2" "y2"
.. ..@ plotOrder: int 1
.. ..@ labpt : num [1:2] 7 3
.. ..@ ID : chr "2"
.. ..@ area : num 4
Can anyone tell me what I'm doing wrong?