1# Copyright 2001-2021 by Roger Bivand 2# 3# 4# Modified by Micah Altman 2010 5 6 7 8poly2nb <- function(pl, row.names=NULL, snap=sqrt(.Machine$double.eps), 9 queen=TRUE, useC=TRUE, foundInBox=NULL) { 10 verbose <- get("verbose", envir = .spdepOptions) 11 .ptime_start <- proc.time() 12 sf <- NULL 13 if (extends(class(pl), "SpatialPolygons")) { 14 sf <- FALSE 15 } else { 16 if (inherits(pl, "sf")) { 17 row.names <- row.names(pl) 18 regid <- NULL 19 pl <- sf::st_geometry(pl) 20 } 21 if (inherits(pl, "sfc")) { 22 if (length(grep("POLYGON", class(pl)[1])) == 0L) 23 pl <- try(st_cast(pl, "MULTIPOLYGON"), silent=TRUE) 24 if (inherits(pl, "try-error")) 25 stop("Polygon geometries required") 26 if (attr(pl, "n_empty") > 0L) 27 stop("Empty geometries found") 28 sf <- TRUE 29 } 30 } 31 if (is.null(sf)) stop("Not a polygon object") 32 33 if (sf) { 34 n <- length(pl) 35 } else { 36 n <- length(slot(pl, "polygons")) 37 } 38 if (n < 1) stop("non-positive number of entities") 39 if (is.null(row.names)) regid <- row.names(pl) 40 else regid <- NULL 41 if (is.null(regid)) { 42 if(is.null(row.names)) regid <- as.character(1:n) 43 else { 44 if(length(row.names) != n) 45 stop("row.names wrong length") 46 else if (length(unique(row.names)) != length(row.names)) 47 stop("non-unique row.names given") 48 else regid <- row.names 49 } 50 } 51 if (snap < 0) snap <- abs(snap) 52# if (snap < .Machine$double.eps) { 53# bsnap <- .Machine$double.eps 54# } else { 55# bsnap <- snap 56# } 57 vbsnap <- c(-snap, snap) 58 if (verbose) cat("handle IDs:", (proc.time() - .ptime_start)[3], "\n") 59 .ptime_start <- proc.time() 60 61 if (sf) { 62 xpl0 <- as.data.frame(sf::st_coordinates(pl)) 63 xpl <- unname(split(xpl0[,1:2], xpl0[, length(xpl0)])) 64# https://github.com/r-spatial/spdep/issues/50 65 xxpl <- lapply(xpl, function(x) do.call("cbind", x[-1,])) 66 } else { 67 xpl <- slot(pl, "polygons") 68 xxpl <- vector(mode="list", length=length(xpl)) 69 for (i in 1:length(xpl)) { 70 xpli <- slot(xpl[[i]], "Polygons") 71 zz <- lapply(xpli, function(j) slot(j, "coords")[-1,]) 72 xxpl[[i]] <- do.call("rbind", zz) 73 } 74 } 75 nrs <- sapply(xxpl, nrow) 76 bb <- t(sapply(xxpl, function(x) { 77 rx <- range(x[,1]) + vbsnap 78 ry <- range(x[,2]) + vbsnap 79 c(rbind(rx, ry)) 80 })) 81 if (verbose) 82 cat("massage polygons:", (proc.time() - .ptime_start)[3], "\n") 83 84 .ptime_start <- proc.time() 85# dbsnap <- as.double(bsnap) 86 dsnap <- as.double(snap) 87 if (is.null(foundInBox)) { 88 if (!sf) { 89 pl0 <- try(st_as_sfc(pl), silent=TRUE) 90 if (inherits(pl0, "try-error")) { 91 warning("poly2nb: spatial indexing abandoned,\n", 92 deparse(substitute(pl)), 93 " could not be coerced to \"sfc\":\n", 94 attr(pl0, "condition")$message) 95 genBBIndex <- function(bb) { 96 n <- nrow(bb) 97 bxv <- as.vector(bb[,c(1,3)]) 98 byv <- as.vector(bb[,c(2,4)]) 99 obxv <- order(bxv) 100 rbxv <- c(1:(n*2))[obxv] 101 mbxv <- match(1:(n*2),obxv) 102 obyv <- order(byv) 103 rbyv <- c(1:(n*2))[obyv] 104 mbyv <- match(1:(n*2),obyv) 105 return(list(bb=bb, bxv=bxv, byv=byv, obxv=obxv, 106 obyv=obyv, mbxv=mbxv, mbyv=mbyv, rbyv=rbyv, 107 rbxv=rbxv)) 108 } 109 BBindex <- genBBIndex(bb) 110 if (verbose) cat("size of BBindex:", object.size(BBindex), 111 "\n") 112 foundInBox <- lapply(1:(n-1), function(i) 113 findInBox(i, BBindex)) 114 } else { 115 pl <- pl0 116 } 117 } 118 if (is.null(foundInBox)) { 119# https://github.com/r-spatial/spdep/issues/65 120 if (!is.na(st_is_longlat(pl)) && 121 st_is_longlat(pl) && sf_use_s2()) { 122 if (packageVersion("sf") < "1.0.4") { 123 fB1 <- st_intersects(pl, sparse=TRUE) 124 } else { 125 fB1 <- st_intersects(pl, sparse=TRUE, model="closed") 126 } 127 } else { 128 cdsnap <- as.double(c(-snap, -snap, snap, snap)) 129 cbb <- t(apply(bb, 1, function(x) x+cdsnap)) 130 envs <- apply(cbb, 1, function(x) 131 st_polygon(list(rbind(x[c(1,2)], x[c(3,2)], x[c(3,4)], 132 x[c(1,4)], x[c(1,2)])))) 133 envs_sfc <- st_as_sfc(envs, crs=st_crs(pl)) 134 fB1 <- st_intersects(envs_sfc, sparse=TRUE) 135 rm(envs_sfc) 136 rm(envs) 137 rm(cbb) 138 } 139 fB1a <- lapply(seq_along(fB1), function(i) 140 {fB1[[i]][fB1[[i]] > i]}) 141 foundInBox <- fB1a[-length(fB1a)] 142 rm(fB1) 143 rm(fB1a) 144 } 145 if (verbose) cat("findInBox:", (proc.time() - .ptime_start)[3]) 146 } 147 stopifnot(is.list(foundInBox)) 148 stopifnot(length(foundInBox) == (n-1L)) 149 stopifnot(all(unlist(sapply(foundInBox, 150 function(x) {if(!is.null(x)) is.integer(x)})))) 151 nfIBB <- sum(sapply(foundInBox, length)) 152 if (verbose) { 153 cat(" list size", nfIBB, "\n") 154 cat("generate foundInBox:", (proc.time() - .ptime_start)[3], "\n") 155 } 156 .ptime_start <- proc.time() 157 158 criterion <- ifelse(queen, 0, 1) 159 if (useC) { 160 ans <- .Call("poly_loop2", as.integer(n), foundInBox, bb, xxpl, 161 as.integer(nrs), as.double(dsnap), as.integer(criterion), 162 as.integer(nfIBB), PACKAGE="spdep") 163 } else { 164 polypoly2 <- function(poly1, nrs1, poly2, nrs2, snap) { 165 if (any(nrs1 == 0 || nrs2 == 0)) return(0L) 166 res <- .Call("polypoly", poly1, nrs1, poly2, 167 nrs2, snap, PACKAGE="spdep") 168 res 169 } 170 171 ans <- vector(mode="list", length=n) 172 for (i in 1:n) ans[[i]] <- integer(0) 173 for (i in 1:(n-1)) { 174 for (j in foundInBox[[i]]) { 175 jhit <- .Call("spOverlap", bb[i,], 176 bb[j,], PACKAGE="spdep") 177 if (jhit > 0) { 178 khit <- 0 179 khit <- polypoly2(xxpl[[i]], nrs[i], xxpl[[j]], 180 nrs[j], dsnap) 181 182 if (khit > criterion) { 183 ans[[i]] <- c(ans[[i]], j) 184 ans[[j]] <- c(ans[[j]], i) 185 } 186 } 187 } 188 } 189 for (i in 1:n) { 190 if (length(ans[[i]]) == 0L) ans[[i]] <- 0L 191 if (length(ans[[i]]) > 1L) ans[[i]] <- sort(ans[[i]]) 192 } 193 } 194 if (verbose) cat("work loop:", (proc.time() - .ptime_start)[3], "\n") 195 .ptime_start <- proc.time() 196 197 class(ans) <- "nb" 198 attr(ans, "region.id") <- regid 199 attr(ans, "call") <- match.call() 200 if (queen) attr(ans, "type") <- "queen" 201 else attr(ans, "type") <- "rook" 202 ans <- sym.attr.nb(ans) 203 if (verbose) cat("done:", (proc.time() - .ptime_start)[3], "\n") 204 .ptime_start <- proc.time() 205 ans 206} 207 208 209# faster findInBox 210 211qintersect<-function(x,y) { 212 # streamlined intersect function for unique vectors 213 as.integer(y[match(x, y, 0L)]) 214} 215 216findInBox<-function(i, sp, bigger=TRUE) { 217 n <- dim(sp$bb)[1] 218 219# use index structure to identify which other BB's fall in i's BB 220# by getting id's of polygons with BBmin_j < BBmax_i, BBmax_j > BBmin_i for x and y 221# then taking the intersection of these four lists of id's 222 223 tmp<-vector(mode="list", length=4) 224 # ! i1 > j3 --> i1 <= j3 225 tmp[[1]] <- sp$rbxv[sp$mbxv[i]:(n*2)] 226 tmp[[1]]<- tmp[[1]][which(tmp[[1]]>n)] - n 227 # ! i2 > j4 --> i2 <= bj4 228 tmp[[2]] <- sp$rbyv[sp$mbyv[i]:(n*2)] 229 tmp[[2]]<- tmp[[2]][which(tmp[[2]]>n)] - n 230 # ! i3 < j1 -> i3 >= j1 231 tmp[[3]] <- sp$rbxv[1:sp$mbxv[i+n]] 232 tmp[[3]] <- tmp[[3]][which(tmp[[3]]<=n)] 233 # ! i4 < j2 -> i4 >= j2 234 tmp[[4]] <- sp$rbyv[1:sp$mbyv[i+n]] 235 tmp[[4]]<- tmp[[4]][which(tmp[[4]]<=n)] 236 237 # for performance, order the comparison of the lists 238 239 lentmp <- order(sapply(tmp,length)) 240 241 # use qintersect, since these are already vectors and unique 242 result <- qintersect(tmp[[lentmp[2]]],tmp[[lentmp[1]]]) 243 result <- qintersect(tmp[[lentmp[3]]],result) 244 result <- qintersect(tmp[[lentmp[4]]],result) 245 246 if (bigger) { 247 result<-result[which(result>i)] 248 } 249 return(sort(result)) 250} 251 252 253 254