1# PBSmapping utilities 2 3SpatialPolygons2PolySet <- function(SpP) { 4 # require(PBSmapping) 5 if (!requireNamespace("PBSmapping", quietly = TRUE)) 6 stop("package PBSmapping required for SpatialPolygons2PolySet") 7 pls <- slot(SpP, "polygons") 8 n <- length(pls) 9 PID <- NULL 10 SID <- NULL 11 POS <- NULL 12 X <- NULL 13 Y <- NULL 14 for (i in 1:n) { 15 srs <- slot(pls[[i]], "Polygons") 16 m <- length(srs) 17 for (j in 1:m) { 18 crds <- slot(srs[[j]], "coords") 19 k <- nrow(crds) 20 PID <- c(PID, rep(i, k)) 21 SID <- c(SID, rep(j, k)) 22# POS <- c(POS, 1:k) 23# Daniel Rodolphe Schlaepfer 140806 24 if(slot(srs[[j]], "hole")) POS <- c(POS, k:1) 25 else POS <- c(POS, 1:k) 26 X <- c(X, crds[,1]) 27 Y <- c(Y, crds[,2]) 28 } 29 } 30 PID <- as.integer(PID) 31 SID <- as.integer(SID) 32 POS <- as.integer(POS) 33 storage.mode(X) <- "double" 34 storage.mode(Y) <- "double" 35 pj <- .pbsproj(SpP) 36 zn <- NULL 37 if (pj == "UTM") { 38 zn <- attr(pj, "zone") 39 attr(pj, "zone") <- NULL 40 } 41 res <- PBSmapping::as.PolySet(data.frame(PID=PID, SID=SID, POS=POS, 42 X=X, Y=Y), projection=pj, zone=zn) 43 res 44} 45 46SpatialLines2PolySet <- function(SL) { 47# require(maps) 48 pls <- slot(SL, "lines") 49 n <- length(pls) 50 PID <- NULL 51 SID <- NULL 52 POS <- NULL 53 X <- NULL 54 Y <- NULL 55 for (i in 1:n) { 56 srs <- slot(pls[[i]], "Lines") 57 m <- length(srs) 58 for (j in 1:m) { 59 crds <- coordinates(srs[[j]]) 60 k <- nrow(crds) 61 PID <- c(PID, rep(i, k)) 62 SID <- c(SID, rep(j, k)) 63 POS <- c(POS, 1:k) 64 X <- c(X, crds[,1]) 65 Y <- c(Y, crds[,2]) 66 } 67 } 68 PID <- as.integer(PID) 69 SID <- as.integer(SID) 70 POS <- as.integer(POS) 71 storage.mode(X) <- "double" 72 storage.mode(Y) <- "double" 73 # require(PBSmapping) 74 if (!requireNamespace("PBSmapping", quietly = TRUE)) 75 stop("package PBSmapping required for SpatialPolygons2PolySet") 76 pj <- .pbsproj(SL) 77 zn <- NULL 78 if (pj == "UTM") { 79 zn <- attr(pj, "zone") 80 attr(pj, "zone") <- NULL 81 } 82 res <- PBSmapping::as.PolySet(data.frame(PID=PID, SID=SID, POS=POS, 83 X=X, Y=Y), projection=pj, zone=zn) 84 res 85} 86 87.pbsproj <- function(Sobj) { 88 p4str <- proj4string(Sobj) 89 if (is.na(p4str)) return("1") 90 res <- grep("longlat", p4str, fixed=TRUE) 91 if (length(res) > 0L) return("LL") 92 res <- regexpr("utm", p4str, fixed=TRUE) 93 if (res > 0) { 94 val <- "UTM" 95 res <- regexpr("+zone=", p4str, fixed=TRUE) 96 sres <- substring(p4str, res+attr(res, "match.length")) 97 zn0 <- regexpr("[[:digit:]]+", sres) 98 attr(val, "zone") <- as.integer(substring(sres, zn0, 99 zn0+attr(zn0, "match.length"))) 100 } else val <- "1" 101 val 102} 103 104PolySet2SpatialPolygons <- function(PS, close_polys=TRUE) { 105 if (!inherits(PS, "PolySet")) stop("not a PolySet object") 106 prj <- attr(PS, "projection") 107 if (is.null(prj)) stop("unknown coordinate reference system") 108 if (prj == "LL") p4s <- "+proj=longlat +ellps=WGS84" 109 else if (prj == "UTM") { 110# apparent change in PBS object attributes 111 zn <- attr(prj, "zone") 112 if (is.null(zn)) zn <- attr(PS, "zone") 113 if (is.null(zn)) stop("no valid zone attribute") 114 p4s <- paste("+proj=utm +ellps=WGS84 +zone=", zn, sep="") 115 } else { 116 p4s <- as.character(NA) 117 warning("unknown coordinate reference system") 118 } 119# stop("unknown coordinate reference system") 110310 120 hasPID <- "PID" %in% names(PS) 121 if (!hasPID) stop("object does not have PID column") 122 res0 <- split(PS, PS$PID) 123 hasSID <- "SID" %in% names(PS) 124 outPolygons <- vector(mode="list", length=length(res0)) 125 if (hasSID) { 126 res1 <- lapply(res0, function(x) split(x, x$SID)) 127 if (close_polys) res1 <- lapply(res1, 128 function(i) lapply(i, function(x) { 129 n <- nrow(x) 130 if (!isTRUE(identical(x$X[1], x$X[n])) || 131 !isTRUE(identical(x$Y[1], x$Y[n]))) rbind(x, x[1,]) 132 else x 133 }) 134 ) 135# extra level added to fix bug found by A Lobos 080413 136 for (i in seq(along=outPolygons)) { 137 outPolygons[[i]] <- Polygons(lapply(res1[[i]], function(x) 138 Polygon(cbind(x$X, x$Y))), ID=names(res1)[i]) 139 } 140# PIDs added as IDs 080413 141 } else { 142 if (close_polys) res0 <- lapply(res0, function(x) { 143 n <- nrow(x) 144 if (!isTRUE(identical(x$X[1], x$X[n])) || 145 !isTRUE(identical(x$Y[1], x$Y[n]))) rbind(x, x[1,]) 146 else x 147 }) 148 for (i in seq(along=outPolygons)) { 149 outPolygons[[i]] <- Polygons(list(Polygon(cbind(res0[[i]]$X, 150 res0[[i]]$Y))), ID=as.character(i)) 151 } 152 } 153 outSP <- SpatialPolygons(outPolygons, proj4string=CRS(p4s)) 154 outSP 155} 156 157PolySet2SpatialLines <- function(PS) { 158 if (!inherits(PS, "PolySet")) stop("not a PolySet object") 159 prj <- attr(PS, "projection") 160 prj <- attr(PS, "projection") 161 if (is.null(prj)) stop("unknown coordinate reference system") 162 if (prj == "LL") p4s <- "+proj=longlat +ellps=WGS84" 163 else if (prj == "UTM") { 164# apparent change in PBS object attributes 165 zn <- attr(prj, "zone") 166 if (is.null(zn)) zn <- attr(PS, "zone") 167 if (is.null(zn)) stop("no valid zone attribute") 168 p4s <- paste("+proj=utm +ellps=WGS84 +zone=", zn, sep="") 169 } else { 170 p4s <- as.character(NA) 171 warning("unknown coordinate reference system") 172 } 173# stop("unknown coordinate reference system") 110310 174 hasPID <- "PID" %in% names(PS) 175 if (!hasPID) stop("object does not have PID column") 176 res0 <- split(PS, PS$PID) 177 hasSID <- "SID" %in% names(PS) 178 outLines <- vector(mode="list", length=length(res0)) 179 if (hasSID) { 180 res1 <- lapply(res0, function(x) split(x, x$SID)) 181 for (i in seq(along=outLines)) { 182 outLines[[i]] <- Lines(lapply(res1[[i]], function(x) 183 Line(cbind(x$X, x$Y))), ID=as.character(i)) 184 } 185 } else { 186 for (i in seq(along=outLines)) { 187 outLines[[i]] <- Lines(lapply(res0[[i]], function(x) 188 Line(cbind(res0[[i]]$X, res0[[i]]$Y))), ID=as.character(i)) 189 } 190 } 191 outSP <- SpatialLines(outLines, proj4string=CRS(p4s)) 192 outSP 193} 194