1# Copyright 2000-2001 (c) Nicholas Lewin-Koh 2# modifications 2001-2005 Roger Bivand, 3# shape2poly based on code by Stephane Dray 4 5plot.polylist <- function(x, col, border=par("fg"), add=FALSE, 6 xlim=NULL, ylim=NULL, xlab="", ylab="", xpd=NULL, 7 density=NULL, angle=45, pbg=NULL, forcefill=TRUE, ...) { 8 .Deprecated("", package="sp", msg="objects other than Spatial objects defined in the sp package are deprecated") 9 if (!inherits(x, "polylist")) stop("Not a polygon list") 10 11 usrpoly <- function(x) { 12 p <- matrix(c(x[1], x[2], x[2], x[1], x[3], x[3], 13 x[4], x[4]), ncol=2) 14 p 15 } 16 if (!add) { 17 maplim <- attr(x, "maplim") 18 if (is.null(maplim)) 19 if (is.null(xlim) || is.null(ylim)) 20 stop("map limits missing") 21 if (is.null(xlim)) xlim <- maplim$x 22 if (is.null(ylim)) ylim <- maplim$y 23 plot(x=xlim, y=ylim, xlim=xlim, ylim=ylim, type="n", 24 asp=1, xlab=xlab, ylab=ylab, ...) 25 if (!forcefill && !is.null(pbg)) { 26 polygon(usrpoly(par("usr")), col = pbg, border = NA) 27 box() 28 } 29 } 30 pO <- attr(x, "plotOrder") 31 if (length(x) < 1L) stop("zero length polylist") 32 if (is.null(pO) || length(x) != length(pO)) pO <- 1:length(x) 33 pO <- as.integer(pO) 34 if (length(pO) != length(unique(pO))) stop("malformed plot order") 35 if (!identical(sort(pO), sort(1:length(x)))) stop("malformed plot order") 36 if (missing(col)) { 37 if (length(density) != length(x)) { 38 density <- rep(density, length(x), length(x)) 39 } 40 if (length(angle) != length(x)) { 41 angle <- rep(angle, length(x), length(x)) 42 } 43 for (j in pO) polygonholes(x[[j]], border=border, 44 xpd=xpd, density=density[j], angle=angle[j], 45 pbg=pbg, forcefill=forcefill) 46 } else { 47 if (length(col) != length(x)) { 48 col <- rep(col, length(x), length(x)) 49 } 50 for (j in pO) 51 polygonholes(x[[j]], col=col[j], border=border, 52 xpd=xpd, pbg=pbg, forcefill=forcefill) 53 } 54 if (!missing(col) | !is.null(density)) { 55 if (forcefill) 56 warning("From next major release, default fill behaviour will change") 57 } 58} 59 60polygonholes <- function(coords, col=NA, border=NULL, xpd=NULL, density=NULL, 61 angle=45, pbg=par("bg"), forcefill=TRUE) { 62 nParts <- attr(coords, "nParts") 63 if (is.null(nParts)) nParts <- 1 64 pFrom <- attr(coords, "pstart")$from 65 if (is.null(pFrom)) pFrom <- 1 66 pTo <- attr(coords, "pstart")$to 67 if (is.null(pTo)) pTo <- dim(coords)[1] 68 if (is.na(col)) hatch <- TRUE 69 else hatch <- FALSE 70 pO <- attr(coords, "plotOrder") 71 if (is.null(pO)) pO <- 1:nParts 72 for (i in pO) { 73 if (hatch) { 74 if (forcefill || attr(coords, "ringDir")[i] == 1) { 75 polygon(coords[pFrom[i]:pTo[i],], 76 border = border, xpd = xpd, 77 density = density, angle = angle) 78 } else { 79 polygon(coords[pFrom[i]:pTo[i],], 80 border = border, xpd = xpd, col=pbg, 81 density = NULL) 82 } 83 } else { 84 if (forcefill || attr(coords, "ringDir")[i] == 1) { 85 polygon(coords[pFrom[i]:pTo[i],], 86 border = border, xpd = xpd, col=col) 87 } else { 88 polygon(coords[pFrom[i]:pTo[i],], 89 border = border, xpd = xpd, col=pbg) 90 } 91 } 92 } 93} 94 95 96plotpolys <- function(pl, bb, col=NA, border=par("fg"), add=FALSE, 97 xlim=NULL, ylim=NULL, ...) { 98 .Deprecated("plot.polylist", package="maptools") 99 if (!inherits(pl, "polylist")) stop("Not a polygon list") 100 if (!add) { 101 maplim <- attr(pl, "maplim") 102 if (is.null(maplim) && missing(bb)) 103 if (is.null(xlim) || is.null(ylim)) 104 stop("map limits missing") 105 if (is.null(xlim)) { 106 if (is.null(maplim)) 107 xlim <- c(min(bb[,1]), max(bb[,3])) 108 else xlim <- maplim$x 109 } 110 if (is.null(ylim)) { 111 if (is.null(maplim)) 112 ylim <- c(min(bb[,2]), max(bb[,4])) 113 else ylim <- maplim$y 114 } 115 plot(x=xlim, y=ylim, xlim=xlim, ylim=ylim, type="n", 116 asp=1, xlab="", ylab="") 117 } 118 if (length(col) != length(pl)) { 119 col <- rep(col, length(pl), length(pl)) 120 } 121 for (j in 1:length(pl)) polygon(pl[[j]], col=col[j], border=border, ...) 122} 123 124shape2poly <- function(shape, region.id=NULL) { 125 if (is.null(shape$shp)) stop("No shp component in this list") 126 if (shape$shp$header$shape.type != 5) stop("Not a polygon shapefile") 127 nrecord <- length(shape$shp$shp) 128 res <- vector(mode="list", length=nrecord) 129 if (is.null(region.id) || length(region.id) != nrecord) { 130 attr(res, "region.id") <- as.character(1:nrecord) 131 } else { 132 attr(res, "region.id") <- as.character(region.id) 133 } 134 np <- integer(nrecord) 135 for (i in 1:nrecord) np[i] <- shape$shp$shp[[i]]$num.parts 136 for (i in 1:nrecord) { 137 if (np[i] > 1) { 138 res[[i]] <- .getMultishape(shape$shp$shp[[i]], np[i]) 139 } else { 140 res[[i]] <- as.matrix(shape$shp$shp[[i]]$points) 141 attr(res[[i]], "pstart") <- list(from=1, 142 to=shape$shp$shp[[i]]$num.points) 143 rownames(res[[i]]) <- NULL 144 colnames(res[[i]]) <- NULL 145 } 146 attr(res[[i]], "nParts") <- np[i] 147 rD <- integer(np[i]) 148 for (j in 1:np[i]) rD[j] <- ringDir(res[[i]], j) 149 attr(res[[i]], "ringDir") <- rD 150 attr(res[[i]], "bbox") <- as.vector(shape$shp$shp[[i]]$box) 151 } 152 153 class(res) <- "polylist" 154 attr(res, "maplim") <- shp2maplim(shape) 155 return(res) 156 157} 158 159shape2lines <- function(shape) { 160 if (is.null(shape$shp)) stop("No shp component in this list") 161 if (shape$shp$header$shape.type != 3) stop("maptype not line/arc") 162 n <- length(shape$shp$shp) 163 res <- vector(mode="list", length=n) 164 nParts <- integer(n) 165 for (i in 1:n) nParts[i] <- shape$shp$shp[[i]]$num.parts 166 for (i in 1:n) { 167 if (nParts[i] > 1) 168 res[[i]] <- .getMultishape(shape$shp$shp[[i]], 169 nParts[i]) 170 else { 171 res[[i]] <- as.matrix(shape$shp$shp[[i]]$points) 172 attr(res[[i]], "pstart") <- list(from=1, 173 to=shape$shp$shp[[i]]$num.points) 174 rownames(res[[i]]) <- NULL 175 colnames(res[[i]]) <- NULL 176 } 177 attr(res[[i]], "nParts") <- nParts[i] 178 } 179 class(res) <- "lineslist" 180 attr(res, "maplim") <- shp2maplim(shape) 181 res 182} 183 184shape2points <- function(shape) { 185 if (is.null(shape$shp)) stop("No shp component in this list") 186 if (shape$shp$header$shape.type != 1) 187 stop("maptype not points") 188# n <- length(shape$shp$shp) 189 res <- shape$shp$shp[,2:3] 190 class(res) <- "Mappoints" 191 attr(res, "maplim") <- shp2maplim(shape) 192 res 193} 194 195 196# based on SHPRingDir_2d, modified to use current ring only, and to strip 197# out last vertex if identical with first 198 199ringDir <- function(xy, ring) { 200 nParts <- attr(xy, "nParts") 201 if (ring > nParts) stop("ring too large") 202 from <- attr(xy, "pstart")$from 203 to <- attr(xy, "pstart")$to 204 a <- xy[from[ring]:to[ring],1] 205 b <- xy[from[ring]:to[ring],2] 206 nvx <- length(b) 207 208 if((a[1] == a[nvx]) && (b[1] == b[nvx])) { 209 a <- a[-nvx] 210 b <- b[-nvx] 211 nvx <- nvx - 1 212 } 213 214 tX <- 0.0 215 dfYMax <- max(b) 216 ti <- 1 217 for (i in 1:nvx) { 218 if (b[i] == dfYMax && a[i] > tX) ti <- i 219 } 220 if ( (ti > 1) & (ti < nvx) ) { 221 dx0 = a[ti-1] - a[ti] 222 dx1 = a[ti+1] - a[ti] 223 dy0 = b[ti-1] - b[ti] 224 dy1 = b[ti+1] - b[ti] 225 } else if (ti == nvx) { 226 dx0 = a[ti-1] - a[ti] 227 dx1 = a[1] - a[ti] 228 dy0 = b[ti-1] - b[ti] 229 dy1 = b[1] - b[ti] 230 } else { 231# /* if the tested vertex is at the origin then continue from 0 (1) */ 232 dx1 = a[2] - a[1] 233 dx0 = a[nvx] - a[1] 234 dy1 = b[2] - b[1] 235 dy0 = b[nvx] - b[1] 236 } 237 v3 = ( (dx0 * dy1) - (dx1 * dy0) ) 238 if ( v3 > 0 ) return (1) 239 else return (-1) 240} 241 242shp2maplim <- function(shape) { 243 if (is.null(shape$shp)) stop("No shp component in this list") 244 mapxlim<-c(shape$shp$header$xmin, shape$shp$header$xmax) 245 mapylim<-c(shape$shp$header$ymin, shape$shp$header$ymax) 246 list(x=mapxlim, y=mapylim) 247} 248 249.getMultishape <- function(shp, nParts) { 250 Pstart <- shp$parts 251 nVerts <- shp$num.points 252 from <- integer(nParts) 253 to <- integer(nParts) 254 from[1] <- 1 255 for (j in 1:nParts) { 256 if (j == nParts) to[j] <- nVerts 257 else { 258 to[j] <- Pstart[j+1] 259 from[j+1] <- to[j]+1 260 } 261 } 262 res <- as.matrix(shp$points[from[1]:to[1],]) 263 if (nParts > 1) { 264 for (j in 2:nParts) { 265 res <- rbind(res, c(NA, NA)) 266 res <- rbind(res, as.matrix(shp$points[from[j]:to[j],])) 267 } 268 } 269 rownames(res) <- NULL 270 colnames(res) <- NULL 271 for (j in 1:nParts) { 272 from[j] <- from[j] + (j-1) 273 to[j] <- to[j] + (j-1) 274 } 275 attr(res, "pstart") <- list(from=from, to=to) 276 res 277} 278 279shape2bbs <- function(shape) { 280 if (is.null(shape$shp)) stop("No shp component in this list") 281 if (shape$shp$header$shape.type != 5) stop("Not a polygon shapefile") 282 n <- length(shape$shp$shp) 283 res <- matrix(0, ncol=4, nrow=n) 284 for (i in 1:n) res[i,] <- as.vector(shape$shp$shp[[i]]$box) 285 res 286} 287 288Map2lines <- function(Map) { 289 .Deprecated("", package="sp", msg="objects other than Spatial objects defined in the sp package are deprecated") 290 if (!inherits(Map, "Map")) stop("not a Map") 291 if (attr(Map$Shapes,'shp.type') != 'arc') 292 stop("maptype not line/arc") 293 n <- attr(Map$Shapes,'nshps') 294 res <- vector(mode="list", length=n) 295 nParts <- integer(n) 296 for (i in 1:n) nParts[i] <- attr(Map$Shapes[[i]], "nParts") 297 for (i in 1:n) { 298 if (nParts[i] > 1) 299 res[[i]] <- .getMultiShp(Map$Shapes[[i]], nParts[i]) 300 else { 301 res[[i]] <- Map$Shapes[[i]]$verts 302 attr(res[[i]], "pstart") <- list(from=1, 303 to=attr(Map$Shapes[[i]], "nVerts")) 304 } 305 attr(res[[i]], "nParts") <- nParts[i] 306 shpID <- attr(Map$Shapes[[i]], "shpID") 307 attr(res[[i]], "shpID") <- ifelse (is.null(shpID), as.integer(NA), shpID) 308 } 309 class(res) <- "lineslist" 310 attr(res, "maplim") <- Map2maplim(Map) 311 res 312} 313 314Map2points <- function(Map) { 315 .Deprecated("", package="sp", msg="objects other than Spatial objects defined in the sp package are deprecated") 316 if (!inherits(Map, "Map")) stop("not a Map") 317# if (class(Map) != "Map") stop("not a Map") 318 if (attr(Map$Shapes,'shp.type') != 'point') 319 stop("maptype not points") 320 n <- attr(Map$Shapes,'nshps') 321 ncols <- 2 322 if (attr(Map$Shapes[[1]], "shp.type") == 11) ncols <- 3 323 res <- matrix(NA, nrow=n, ncol=ncols) 324 shpIDs <- integer(n) 325 for (i in 1:n) { 326 res[i,] <- Map$Shapes[[i]]$verts 327 shpID <- attr(Map$Shapes[[i]], "shpID") 328 shpIDs[i] <- ifelse (is.null(shpID), as.integer(NA), shpID) 329 } 330 class(res) <- "Mappoints" 331 attr(res, "maplim") <- Map2maplim(Map) 332 attr(res, "shpID") <- shpIDs 333 res 334} 335 336Map2poly <- function(Map, region.id=NULL, quiet=TRUE) { 337 .Deprecated("", package="sp", msg="objects other than Spatial objects defined in the sp package are deprecated") 338 if (!inherits(Map, "Map")) stop("not a Map") 339# if (class(Map) != "Map") stop("not a Map") 340 if (attr(Map$Shapes,'shp.type') != 'poly') 341 stop("maptype not poly") 342 if (!is.null(region.id)) 343 if (length(region.id) != length(unique(region.id))) 344 stop("region.id not unique") 345 res <- .get.polylist1(Map=Map, region.id=region.id, quiet=quiet) 346 attr(res, "maplim") <- Map2maplim(Map) 347 area <- sapply(res, function(x) attr(x, "area")) 348 pO <- order(area, decreasing=TRUE) 349 after <- as.integer(rep(NA, attr(Map$Shapes,'nshps'))) 350 351 attr(res, "after") <- after 352 attr(res, "plotOrder") <- pO 353 354 nD <- unique(sapply(res, function(x) dim(x)[2])) 355 if (length(nD) > 1L) stop("multiple dimension polylist components") 356 nD <- as.integer(nD) 357 attr(res, "nDims") <- nD 358 359 res 360} 361 362 363Map2poly1 <- function(Map, region.id=NULL, raw=TRUE) { 364 .Deprecated("", package="sp", msg="objects other than Spatial objects defined in the sp package are deprecated") 365 if (!inherits(Map, "Map")) stop("not a Map") 366# if (class(Map) != "Map") stop("not a Map") 367 if (attr(Map$Shapes,'shp.type') != 'poly') 368 stop("maptype not poly") 369 if (!is.null(region.id)) 370 if (length(region.id) != length(unique(region.id))) 371 stop("region.id not unique") 372 res <- .get.polylist(Map=Map, region.id=region.id, raw=raw) 373 attr(res, "maplim") <- Map2maplim(Map) 374 pO <- as.integer(1:attr(Map$Shapes,'nshps')) 375 after <- as.integer(rep(NA, attr(Map$Shapes,'nshps'))) 376 rD <- sapply(res, function(x) 377 attr(x, "ringDir")[attr(x, "plotOrder")[1]]) 378 r1 <- .mtInsiders(res, rD) 379 if (!all(sapply(r1, is.null))) { 380 lres <- .mtlbuild(.mtafters(r1), rD) 381 pO <- lres$pO 382 after <- lres$after 383 } 384# pO <- as.integer(1:attr(Map$Shapes,'nshps')) 385# after <- as.integer(rep(NA, attr(Map$Shapes,'nshps'))) 386# r1 <- .mtInsiders(res) 387# if (!all(sapply(r1, is.null))) { 388# after <- as.integer(sapply(r1, 389# function(x) ifelse(is.null(x), NA, max(x)))) 390# pO <- order(after, na.last=FALSE) 391# } 392 attr(res, "after") <- after 393 attr(res, "plotOrder") <- pO 394 if (!raw) { 395 rD <- sapply(res, function(x) attr(x, 396 "ringDir")[which(attr(x, "plotOrder") == 1)]) 397 if (any((rD == -1) & is.na(after))) { 398 oddCC <- which((rD == -1) & is.na(after)) 399 for (i in oddCC) { 400 tgt <- which(attr(res[[i]], "plotOrder") == 1) 401 nParts <- attr(res[[i]], "nParts") 402 tmp <- as.matrix(res[[i]]) 403 from <- attr(res[[i]], "pstart")$from[tgt] 404 to <- attr(res[[i]], "pstart")$to[tgt] 405 tmp[from:to,] <- res[[i]][to:from, ] 406 attributes(tmp) <- attributes(res[[i]]) 407 rD <- vector(length=nParts, mode="integer") 408 for (j in 1:nParts) rD[j] <- ringDir(tmp, j) 409 attr(tmp, "ringDir") <- rD 410 res[[i]] <- tmp 411 warning(paste("ring direction changed in polygon", i)) 412 } 413 } 414 } 415 if (raw) warning("From next release, default hole handling will change") 416 res 417} 418 419.mtInsiders <- function(pl, rD) { 420 bbox1 <- function(x) { 421 r1 <- range(x[,1], na.rm=TRUE) 422 r2 <- range(x[,2], na.rm=TRUE) 423 res <- c(r1[1], r2[1], r1[2], r2[2]) 424 res 425 } 426 427 n <- length(pl) 428 bbs <- matrix(0, nrow=n, ncol=4) 429 for (i in 1:n) bbs[i,] <- bbox1(pl[[i]]) 430 storage.mode(bbs) <- "double" 431 res <- .Call("mtInsiders", as.integer(n), bbs, 432 PACKAGE="maptools") 433 res1 <- vector(mode="list", length=n) 434 435 for (i in 1:n) { 436 if (!is.null(res[[i]])) { 437 ri <- res[[i]] 438 ixc <- pl[[i]][1,1] 439 iyc <- pl[[i]][1,2] 440 int <- logical(length(ri)) 441 for (j in 1:length(ri)) { 442 xj <- pl[[ri[j]]] 443 jxc <- na.omit(xj[,1]) 444 jyc <- na.omit(xj[,2]) 445 pip <- mt.point.in.polygon(ixc, 446 iyc, jxc, jyc) 447 int[j] <- ((pip == 1) | (pip > 1)) 448# int[j] <- ((pip == 1) | 449# ((pip > 1) & ((rD[i] == 1) & 450# (rD[ri[j]] == -1)))) 451 452 } 453 rj <- ri[int] 454 if (length(rj) > 0L) { 455 res1[[i]] <- as.integer(rj) 456 } 457 } 458 } 459 res1 460} 461 462.mtafters <- function(rl) { 463 464# argument is output from .insiders() - a list with either NULL components 465# (not included in any other polygon) or lists of polygons in which the polygon 466# in question is included; outputs a from:to matrix 467 468 n <- length(rl) 469 res <- NULL 470 for (i in 1:n) { 471 if (is.null(rl[[i]])) 472 res <- rbind(res, c(i, NA)) 473 else { 474 for (j in 1:length(rl[[i]])) { 475 res <- rbind(res, c(i, rl[[i]][j])) 476 } 477 } 478 } 479 res 480} 481 482.mtlbuild1 <- function(x) { 483 484# reverse list builder with from:to matrix as argument, used to try to find 485# circularities 486 487 lx <- vector(mode="list", length=length(unique(x[,1]))) 488 rle.x <- rle(x[,1]) 489 cs1.x <- cumsum(rle.x$lengths) 490 cs0.x <- c(1, cs1.x[1:(length(lx)-1)]+1) 491 ii <- 1 492 for (i in 1:length(lx)) { 493 if (rle.x$value[ii] == i) { 494 lx[[i]] <- as.integer(x[cs0.x[ii]:cs1.x[ii],2]) 495 ii <- ii+1 496 } 497 } 498 lx 499} 500 501.mtcircs <- function(x) { 502 503# try to find circularities from reverse list as argument (polygons reported 504# as being inside each other despite checking ring direction in .insiders); 505# only the first loop will be run in normal cases 506 507 res <- NULL 508 for (i in 1:length(x)) { 509 if (!is.na(match(i, unlist(x[x[[i]]])))) { 510 hits <- rep(FALSE, length(x[[i]])) 511 for (j in 1:length(hits)) { 512 jj <- x[[i]][j] 513 hits[j] <- (i %in% x[[jj]]) 514 } 515 if (length(which(hits)) > 1L) stop("multiple circulars") 516 pair <- c(i, x[[i]][hits]) 517 res <- rbind(res, pair) 518 } 519 } 520 res1 <- NULL 521 if (!is.null(res)) { 522 if (nrow(res) %% 2 != 0) stop("odd circulars") 523 gone <- rep(FALSE, nrow(res)) 524 for (i in 1:nrow(res)) { 525 if (!gone[i]) { 526 from <- res[i,1] 527 to <- res[i,2] 528 hit <- match(from, res[,2]) 529 if (!gone[hit]) { 530 if (res[hit,1] != to) 531 stop("mismatched circular") 532 res1 <- rbind(res1, c(from, to)) 533 gone[i] <- TRUE 534 } 535 } 536 } 537 } 538 res1 539} 540 541.mtlbuild <- function(x, rD) { 542 543# list analysis of matrix output from .afters combined with current ring 544# directions (which may be quite wrong) to generate a plot order and 545# vector of afters (NA for no dependency, 1 for dependency on being plotted 546# after another polygon) 547 548 ids <- x[,1] 549 ins <- x[,2] 550 n <- length(unique(ids)) 551 nas <- which(is.na(ins)) 552 ntop <- length(nas) 553 pO <- vector(length=n, mode="integer") 554 after <- rep(as.integer(NA), length=n) 555 gone <- rep(FALSE, n) 556 j <- 1 557 for (i in 1:ntop) { 558 ii <- ids[nas[i]] 559 if (!gone[ii]) { 560 gone[ii] <- TRUE 561 pO[j] <- ii 562#cat("level 1:", j, ii, "\n") 563 j <- j+1 564 } else warning(paste("level 1 circularity at", ii)) 565 ihits <- which(ins == ii) 566 567# for each top level (not inside any other) polygon, check to see if any 568# polygons are inside it, and insert orders to match; from outer to deepest in; 569# the gone vector is used to avoid multiple assignments to the plot 570# order list that can happen with circularity 571 572 if (length(ihits) > 0L) { 573 tihits <- ids[ihits] 574 rtihits <- rle(ids[ids %in%tihits]) 575 o <- order(rtihits$lengths) 576 for (jj in 1:length(rtihits$values)) { 577 jjj <- rtihits$values[o][jj] 578 if (!gone[jjj]) { 579 gone[jjj] <- TRUE 580 pO[j] <- jjj 581cat("level 2:", j, ii, "\n") 582 j <- j+1 583 } else warning(paste("level 2 circularity at", 584 jjj)) 585 after[jjj] <- as.integer(1) 586 } 587 } 588 } 589 xcircs <- .mtcircs(.mtlbuild1(x)) 590 591# Further attempts to trap circularities, possibly no longer needed, first 592# introduced before point-in-polygon test added to .insiders; TODO check 593# whether is.null(xcircs) is always TRUE 594 595 if (!is.null(xcircs)) { 596 for (i in 1:nrow(xcircs)) { 597 from <- xcircs[i,1] 598 to <- xcircs[i,2] 599 rDfrom <- rD[from] 600 rDto <- rD[to] 601 pOfrom <- which(pO == from) 602 pOto <- which(pO == to) 603 if (rDfrom == 1) { 604 if (pOfrom < pOto) { 605 pO[pOto] <- from 606 pO[pOfrom] <- to 607 } 608 } 609 if (rDto == 1) { 610 if (pOfrom > pOto) { 611 pO[pOto] <- from 612 pO[pOfrom] <- to 613 } 614 } 615 } 616 } 617 list(pO=pO, after=after) 618} 619 620.get.polylist <- function(Map, region.id=NULL, raw=TRUE) { 621 n <- attr(Map$Shapes,'nshps') 622 res <- vector(mode="list", length=n) 623 nParts <- integer(n) 624 for (i in 1:n) nParts[i] <- attr(Map$Shapes[[i]], "nParts") 625 for (i in 1:n) { 626 if (nParts[i] > 1) { 627 res[[i]] <- .getMultiShp(Map$Shapes[[i]], nParts[i], 628 raw=raw) 629 } else { 630 res[[i]] <- Map$Shapes[[i]]$verts 631 attr(res[[i]], "pstart") <- list(from=as.integer(1), 632 to=as.integer(nrow(Map$Shapes[[i]]$verts))) 633# attr(Map$Shapes[[i]], "nVerts")) 634 attr(res[[i]], "after") <- 1 635 attr(res[[i]], "plotOrder") <- 1 636 attr(res[[i]], "bbox") <- 637 as.vector(attr(Map$Shapes[[i]], "bbox")) 638 attr(res[[i]], "RingDir") <- 639 as.vector(attr(Map$Shapes[[i]], "RingDir")) 640 attr(res[[i]], "nParts") <- nParts[i] 641 attr(res[[i]], "ringDir") <- ringDir(res[[i]], 1) 642 } 643# attr(res[[i]], "shpID") <- attr(Map$Shapes[[i]], "shpID") 644 shpID <- attr(Map$Shapes[[i]], "shpID") 645 attr(res[[i]], "shpID") <- ifelse (is.null(shpID), as.integer(NA), shpID) } 646 if (is.null(region.id) || length(region.id) != n) { 647 attr(res, "region.id") <- as.character(1:n) 648 } else { 649 attr(res, "region.id") <- as.character(region.id) 650 } 651 class(res) <- "polylist" 652 invisible(res) 653} 654.get.polylist1 <- function(Map, region.id=NULL, quiet=TRUE) { 655 n <- attr(Map$Shapes,'nshps') 656 res <- vector(mode="list", length=n) 657 nParts <- integer(n) 658 for (i in 1:n) nParts[i] <- attr(Map$Shapes[[i]], "nParts") 659 for (i in 1:n) { 660 if (nParts[i] > 1) { 661 res[[i]] <- .getMultiShp1(Map$Shapes[[i]], nParts[i], 662 IID=i, quiet=quiet) 663 } else { 664 res[[i]] <- Map$Shapes[[i]]$verts 665 rD <- .ringDirxy(res[[i]]) 666 if (rD != 1) { 667 res[[i]] <- res[[i]][nrow(res[[i]]):1,] 668 if (!quiet) warning(paste( 669 "Ring direction changed in shape", i)) 670 } 671 attr(res[[i]], "pstart") <- list(from=as.integer(1), 672 to=as.integer(nrow(Map$Shapes[[i]]$verts))) 673# attr(Map$Shapes[[i]], "nVerts")) 674 attr(res[[i]], "after") <- 1 675 attr(res[[i]], "plotOrder") <- 1 676 attr(res[[i]], "bbox") <- 677 as.vector(attr(Map$Shapes[[i]], "bbox")) 678 attr(res[[i]], "RingDir") <- 679 as.vector(attr(Map$Shapes[[i]], "RingDir")) 680 attr(res[[i]], "nParts") <- nParts[i] 681 attr(res[[i]], "ringDir") <- .ringDirxy(res[[i]]) 682 cents <- .RingCentrd_2d(res[[i]]) 683 attr(res[[i]], "area") <- cents$area 684 attr(res[[i]], "centroid") <- list(x=cents$xc, 685 y=cents$yc) 686 } 687# attr(res[[i]], "shpID") <- attr(Map$Shapes[[i]], "shpID") 688 shpID <- attr(Map$Shapes[[i]], "shpID") 689 attr(res[[i]], "shpID") <- ifelse (is.null(shpID), as.integer(NA), shpID) } 690 if (is.null(region.id) || length(region.id) != n) { 691 attr(res, "region.id") <- as.character(1:n) 692 } else { 693 attr(res, "region.id") <- as.character(region.id) 694 } 695 class(res) <- "polylist" 696 invisible(res) 697} 698 699MapShapeIds <- function(Map) { 700 if (!inherits(Map, "Map")) stop("not a Map") 701# if (class(Map) != "Map") stop("not a Map") 702 sapply(Map$Shapes, function(x) attr(x, "shpID")) 703} 704.getMultiShp1 <- function(shp, nParts, IID, quiet=TRUE) { 705 Pstart <- shp$Pstart 706# nVerts <- attr(shp, "nVerts") 707 nVerts <- nrow(shp$verts) 708 from <- integer(nParts) 709 to <- integer(nParts) 710 area <- double(nParts) 711 xc <- double(nParts) 712 yc <- double(nParts) 713 from[1] <- 1 714 for (j in 1:nParts) { 715 if (j == nParts) to[j] <- nVerts 716 else { 717 to[j] <- Pstart[j+1] 718 from[j+1] <- to[j]+1 719 } 720 } 721 poly <- shp$verts[from[1]:to[1],] 722 rD <- .ringDirxy(poly) 723 if (rD != 1) { 724 poly <- poly[nrow(poly):1,] 725 if (!quiet) warning(paste( 726 "Ring direction changed in shape", IID, "ring 1")) 727 } 728 cents <- .RingCentrd_2d(poly) 729 area[1] <- cents$area 730 xc[1] <- cents$xc 731 yc[1] <- cents$yc 732 res <- poly 733 if (nParts > 1) { 734 for (j in 2:nParts) { 735 res <- rbind(res, c(NA, NA)) 736 poly <- shp$verts[from[j]:to[j],] 737 rD <- .ringDirxy(poly) 738 if (rD != 1) { 739 poly <- poly[nrow(poly):1,] 740 if (!quiet) warning(paste( 741 "Ring direction changed in shape", IID, "ring", j)) 742 } 743 cents <- .RingCentrd_2d(poly) 744 area[j] <- cents$area 745 xc[j] <- cents$xc 746 yc[j] <- cents$yc 747 res <- rbind(res, poly) 748 } 749 } 750 for (j in 1:nParts) { 751 from[j] <- from[j] + (j-1) 752 to[j] <- to[j] + (j-1) 753 } 754 attr(res, "nParts") <- nParts 755 attr(res, "pstart") <- list(from=as.integer(from), to=as.integer(to)) 756 attr(res, "bbox") <- as.vector(attr(shp, "bbox")) 757 attr(res, "RingDir") <- as.vector(attr(shp, "RingDir")) 758 rD <- integer(nParts) 759 for (j in 1:nParts) rD[j] <- ringDir(res, j) 760 attr(res, "ringDir") <- rD 761 pO <- order(area, decreasing=TRUE) 762 after <- as.integer(rep(NA, nParts)) 763 attr(res, "centroid") <- list(x=xc[pO[1]], y=yc[pO[1]]) 764 765 attr(res, "after") <- after 766 attr(res, "plotOrder") <- pO 767 attr(res, "area") <- sum(area) 768 769 res 770} 771 772.getMultiShp <- function(shp, nParts, raw=TRUE) { 773 Pstart <- shp$Pstart 774# nVerts <- attr(shp, "nVerts") 775 nVerts <- nrow(shp$verts) 776 from <- integer(nParts) 777 to <- integer(nParts) 778 from[1] <- 1 779 for (j in 1:nParts) { 780 if (j == nParts) to[j] <- nVerts 781 else { 782 to[j] <- Pstart[j+1] 783 from[j+1] <- to[j]+1 784 } 785 } 786 res <- shp$verts[from[1]:to[1],] 787 if (nParts > 1) { 788 for (j in 2:nParts) { 789 res <- rbind(res, c(NA, NA)) 790 res <- rbind(res, shp$verts[from[j]:to[j],]) 791 } 792 } 793 for (j in 1:nParts) { 794 from[j] <- from[j] + (j-1) 795 to[j] <- to[j] + (j-1) 796 } 797 attr(res, "nParts") <- nParts 798 attr(res, "pstart") <- list(from=as.integer(from), to=as.integer(to)) 799 attr(res, "bbox") <- as.vector(attr(shp, "bbox")) 800 attr(res, "RingDir") <- as.vector(attr(shp, "RingDir")) 801 rD <- integer(nParts) 802 for (j in 1:nParts) rD[j] <- ringDir(res, j) 803 attr(res, "ringDir") <- rD 804 pO <- as.integer(1:nParts) 805 after <- as.integer(rep(NA, nParts)) 806 res1 <- vector(mode="list", length=nParts) 807 for (i in 1:nParts) res1[[i]] <- res[from[i]:to[i],] 808 r1 <- .mtInsiders(res1, rD) 809 if (!all(sapply(r1, is.null))) { 810 lres <- .mtlbuild(.mtafters(r1), rD) 811 pO <- lres$pO 812 after <- lres$after 813 } 814# r1 <- .mtInsiders(res1) 815# if (!all(sapply(r1, is.null))) { 816# after <- as.integer(sapply(r1, 817# function(x) ifelse(is.null(x), NA, max(x)))) 818# pO <- order(after, na.last=FALSE) 819# } 820 821 attr(res, "after") <- after 822 attr(res, "plotOrder") <- pO 823 if (!raw) { 824 top <- which(pO == 1) 825 if (any((rD[-top] == -1) & is.na(after[-top]))) { 826 oddCC <- which((rD == -1) & is.na(after)) 827 for (i in oddCC) { 828 if (i != top) { 829# from1 <- from[i] 830# to1 <- to[i] 831 res[from[i]:to[i],] <- res[to[i]:from[i],] 832 attr(res, "ringDir")[i] <- ringDir(res, i) 833 warning(paste("ring direction changed in subpolygon")) 834 } 835 } 836 } 837 838 839 } 840 res 841} 842 843.get.polybbs <- function(Map) { 844 n <- length(Map$Shapes) 845 res <- matrix(0, ncol=4, nrow=n) 846 for (i in 1:n) res[i,] <- attr(Map$Shapes[[i]], "bbox") 847 res 848} 849 850Map2bbs <- function(Map) { 851 if (!inherits(Map, "Map")) stop("not a Map") 852# if (class(Map) != "Map") stop("not a Map") 853 if (attr(Map$Shapes,'shp.type') != 'poly') 854 stop("maptype not poly") 855 res <- .get.polybbs(Map) 856 res 857} 858 859Map2maplim <- function(Map) { 860 if (!inherits(Map, "Map")) stop("not a Map") 861# if (class(Map) != "Map") stop("not a Map") 862 mapxlim<-c(attr(Map$Shapes, 'minbb')[1], 863 attr(Map$Shapes, 'maxbb')[1]) 864 mapylim<-c(attr(Map$Shapes, 'minbb')[2], 865 attr(Map$Shapes, 'maxbb')[2]) 866 list(x=mapxlim, y=mapylim) 867} 868 869 870convert.pl <- function(pl) { 871 if (!inherits(pl, "multiparts")) stop("not a mulitpart polylist") 872 res <- vector(mode="list", length=length(pl)) 873 for (i in 1:length(pl)) { 874 lp <- length(pl[[i]]) 875 res[[i]] <- pl[[i]][[1]] 876 if (lp > 1) { 877 for (j in 2:lp) { 878 res[[i]] <- rbind(res[[i]], c(NA, NA)) 879 res[[i]] <- rbind(res[[i]], pl[[i]][[j]]) 880 } 881 } 882 } 883 if (!is.null(attr(pl, "region.id"))) 884 attr(res, "region.id") <- attr(pl, "region.id") 885 class(res) <- "polylist" 886 res 887} 888 889.RingCentrd_2d <- function(plmat) { 890 nVert <- nrow(plmat) 891 storage.mode(plmat) <- "double" 892 res <- .C("RFindCG", as.integer(nVert), plmat[,1], 893 plmat[,2], as.double(0), as.double(0), 894 as.double(0), PACKAGE="maptools") 895 896# x_base <- plmat[1,1] 897# y_base <- plmat[1,2] 898# Cy_accum <- 0.0 899# Cx_accum <- 0.0 900# Area <- 0.0 901# ppx <- plmat[2,1] - x_base 902# ppy <- plmat[2,2] - y_base 903# for (iv in 2:(nVert-2)) { 904# x = plmat[iv,1] - x_base 905# y = plmat[iv,2] - y_base 906# dx_Area <- ((x * ppy) - (y * ppx)) * 0.5 907# Area <- Area + dx_Area 908# Cx_accum <- Cx_accum + ( ppx + x ) * dx_Area 909# Cy_accum <- Cy_accum + ( ppy + y ) * dx_Area 910# ppx <- x 911# ppy <- y 912# } 913# xc <- (Cx_accum / (Area * 3)) + x_base 914# yc <- (Cy_accum / (Area * 3)) + y_base 915 list(xc=res[[4]], yc=res[[5]], area=abs(res[[6]])) 916} 917 918.ringDirxy <- function(xy) { 919 a <- xy[,1] 920 b <- xy[,2] 921 storage.mode(a) <- "double" 922 storage.mode(b) <- "double" 923 nvx <- length(b) 924 925# if((a[1] == a[nvx]) && (b[1] == b[nvx])) { 926# a <- a[-nvx] 927# b <- b[-nvx] 928# nvx <- nvx - 1 929# } 930# 931# tX <- 0.0 932# dfYMax <- max(b) 933# ti <- 1 934# for (i in 1:nvx) { 935# if (b[i] == dfYMax && a[i] > tX) ti <- i 936# } 937# if ( (ti > 1) & (ti < nvx) ) { 938# dx0 = a[ti-1] - a[ti] 939# dx1 = a[ti+1] - a[ti] 940# dy0 = b[ti-1] - b[ti] 941# dy1 = b[ti+1] - b[ti] 942# } else if (ti == nvx) { 943# dx0 = a[ti-1] - a[ti] 944# dx1 = a[1] - a[ti] 945# dy0 = b[ti-1] - b[ti] 946# dy1 = b[1] - b[ti] 947# } else { 948# /* if the tested vertex is at the origin then continue from 0 (1) */ 949# dx1 = a[2] - a[1] 950# dx0 = a[nvx] - a[1] 951# dy1 = b[2] - b[1] 952# dy0 = b[nvx] - b[1] 953# } 954# v3 = ( (dx0 * dy1) - (dx1 * dy0) ) 955 956 res <- .C("RFindCG", as.integer(nvx), a, b, 957 as.double(0), as.double(0), as.double(0), PACKAGE="maptools") 958 959 if ( res[[6]] > 0 ) return(as.integer(-1)) 960 else return(as.integer(1)) 961} 962 963