I have a spatialPolygonsDataFrame consisting of 3 polygons. The third polygon has the same shape as the first but has a hole where the second polygon is located.
I built the hole in the using the answer from another question (How to add a hole to a polygon within a SpatialPolygonsDataFrame?).
library(raster)
library(sp)
# create rasters and store them in a list
r1 <- raster(xmn=1, xmx=5, ymn=1, ymx=5, nrows=4, ncols=4)
r1[] <- 1:length(r1)
# create SpatialPolygonsDataFrame
Sr1 = Polygon(cbind(c(1,5,4,1,1),c(1,2,5,4,1)))
Sr2 = Polygon(cbind(c(2,4,3,2),c(3,2,4,3)))
SpP = SpatialPolygons(list(Polygons(list(Sr1), "s1"), Polygons(list(Sr2), "s2")),
1:2)
dat = data.frame(ID = c("s1", "s2"), value = c("a", "b"))
row.names(dat) <- c("s1", "s2")
p <- SpatialPolygonsDataFrame(SpP, data = dat,
match.ID = TRUE)
AddHoleToPolygon <-function(poly,hole){
# invert the coordinates for Polygons to flag it as a hole
coordsHole <- hole@polygons[[1]]@Polygons[[1]]@coords
newHole <- Polygon(coordsHole,hole=TRUE)
# punch the hole in the main poly
listPol <- poly@polygons[[1]]@Polygons
listPol[[length(listPol)+1]] <- newHole
punch <- Polygons(listPol,poly@polygons[[1]]@ID)
# make the polygon a SpatialPolygonsDataFrame as the entry
new <- SpatialPolygons(list(punch),proj4string=poly@proj4string)
new <- SpatialPolygonsDataFrame(new,data=as(poly,"data.frame"))
return(new)
}
punchedPoly <-AddHoleToPolygon(p[1,],p[2,])
p1 <- rbind(p, punchedPoly, makeUniqueIDs = TRUE)
p1 <- p1[2:3,]
When I use mask() to "crop" the raster r1, then the hole is created, although the triangular polygon has a value and indeed is not a real hole. But it gets "overridden" by the third polygon with the hole:
masked_hole <- mask(r1, p1)
plot(masked_hole)
When I change the order of the polygons, then no hole is created:
m3 <- mask(r1, p1[c(2,1),])
plot(m3)
The function rasterize is affected in the same manner:
r2 <- rasterize(p1, r1, field = "value")
plot(r2)
r3 <- rasterize(p1[c(2,1),], r1, field = "value")
plot(r3)
In my real data I have holes where there are no "filling" polygons and those ones I want to keep as holes.
How can I fix the spatialPolygonsDataFrame for polygons that are creating holes where there are none?
How can I fix this issue without reordering but "transform" the hole-creating polygons?
It was a bug in the raster package which has been fixed meanwhile (see https://github.com/rspatial/raster/issues/60).