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