1 2 3 4.gDif <- function(x, y, type='polygons') { 5 xln <- length(x) 6 yln <- length(y) 7 if (xln==0 | yln==0) { 8 return(x) 9 } 10 rn <- row.names(x) 11 for (i in xln:1) { 12 z <- x[i,] 13 for (j in 1:yln) { 14 z <- rgeos::gDifference(z, y[j,], drop_lower_td=TRUE) 15 if (is.null(z)) { 16 break 17 } 18 } 19 if (is.null(z)) { 20 x <- x[-i,] 21 rn <- rn[-i] 22 } else { 23 if (type=='polygons') { 24 x@polygons[i] <- z@polygons 25 } else { 26 x@lines[i] <- z@lines 27 } 28 } 29 if (length(rn) > 0) { 30 row.names(x) <- rn 31 } 32 } 33 if ((type=='polygons') & (length(x) > 0)) { 34 35 w <- getOption('warn') 36 on.exit(options('warn' = w)) 37 options('warn'=-1) 38 39 j <- rgeos::gIsValid(x, byid=TRUE, reason=FALSE) 40 #j <- which(gArea(x, byid=TRUE) > 0) 41 if (!all(j)) { 42 bad <- which(!j) 43 for (i in bad) { 44 # it could be that a part of a polygon is a sliver, but that other parts are OK 45 a <- sp::disaggregate(x[i, ]) 46 if (length(a) > 1) { 47 jj <- which(rgeos::gIsValid(a, byid=TRUE, reason=FALSE)) 48 a <- a[jj, ] 49 if (length(a) > 0) { 50 x@polygons[i] <- aggregate(a)@polygons 51 j[i] <- TRUE 52 } 53 } 54 } 55 x <- x[j,] 56 rn <- rn[j] 57 } 58 59 if (length(rn) > 0) { 60 row.names(x) <- rn 61 } 62 } 63 x 64} 65 66 67setMethod(erase, signature(x='SpatialPolygons', y='SpatialPolygons'), 68 function(x, y, ...){ 69 70 valgeos <- .checkGEOS(); on.exit(rgeos::set_RGEOS_CheckValidity(valgeos)) 71 72 prj <- x@proj4string 73 if (is.na(prj)) prj <- y@proj4string 74 x@proj4string <- sp::CRS(as.character(NA)) 75 y@proj4string <- sp::CRS(as.character(NA)) 76 77 if (!.hasSlot(x, "data")) { 78 d <- data.frame(erase_dissolve_ID=1:length(x)) 79 rownames(d) <- row.names(x) 80 x <- sp::SpatialPolygonsDataFrame(x, data=d) 81 dropframe <- TRUE 82 } else { 83 dropframe <- FALSE 84 x$erase_dissolve_ID <- 1:nrow(x) 85 } 86 87 y <- aggregate(y) 88 89 int <- rgeos::gIntersects(x, y, byid=TRUE) 90 int1 <- apply(int, 2, any) 91 int2 <- apply(int, 1, any) 92 93 if (sum(int1) == 0) { # no intersections 94 return(x) 95 } 96 97 if (all(int1)) { 98 part1 <- NULL 99 } else { 100 part1 <- x[!int1,] 101 } 102 part2 <- .gDif(x[int1,], y[int2,]) 103 104 part2 <- sp::SpatialPolygonsDataFrame(part2, x@data[match(row.names(part2), rownames(x@data)), ,drop=FALSE]) 105 if (!is.null(part1)) { 106 part2 <- rbind(part1, part2) 107 } 108 109 if (length(part2@polygons) > 1) { 110 part2 <- aggregate(part2, colnames(part2@data)) 111 } 112 part2@proj4string <- prj 113 114 if (dropframe) { 115 return( as(part2, 'SpatialPolygons') ) 116 } else { 117 part2@data$erase_dissolve_ID <- NULL 118 return( part2 ) 119 } 120 } 121) 122 123setMethod(erase, signature(x='SpatialLines', y='SpatialPolygons'), 124 function(x, y, ...){ 125 126 valgeos <- .checkGEOS(); on.exit(rgeos::set_RGEOS_CheckValidity(valgeos)) 127 prj <- x@proj4string 128 if (is.na(prj)) prj <- y@proj4string 129 x@proj4string <- sp::CRS(as.character(NA)) 130 y@proj4string <- sp::CRS(as.character(NA)) 131 132 133 if (!.hasSlot(x, 'data')) { 134 d <- data.frame(ID=1:length(x)) 135 rownames(d) <- row.names(x) 136 x <- sp::SpatialLinesDataFrame(x, data=d) 137 dropframe <- TRUE 138 } else { 139 dropframe <- FALSE 140 } 141 142 y <- aggregate(y) 143 144 int <- rgeos::gIntersects(x, y, byid=TRUE) 145 int1 <- apply(int, 2, any) 146 int2 <- apply(int, 1, any) 147 148 if (sum(int1) == 0) { # no intersections 149 return(x) 150 } 151 152 if (all(int1)) { 153 part1 <- NULL 154 } else { 155 part1 <- x[!int1,] 156 } 157 part2 <- .gDif(x[int1,], y[int2,], 'lines') 158 159 part2 <- sp::SpatialLinesDataFrame(part2, x@data[match(row.names(part2), rownames(x@data)), ,drop=FALSE], match.ID = FALSE) 160 if (!is.null(part1)) { 161 part2 <- rbind(part1, part2) 162 } 163 164 if (length(part2@lines) > 1) { 165 part2 <- aggregate(part2, colnames(part2@data)) 166 } 167 168 part2@proj4string <- prj 169 if (dropframe) { 170 return( as(part2, 'SpatialLines') ) 171 } else { 172 return( part2 ) 173 } 174 175 } 176) 177 178