-
Notifications
You must be signed in to change notification settings - Fork 96
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
8 changed files
with
308 additions
and
3 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,239 @@ | ||
| # Author: Robert Hijmans | ||
| # December 2024 | ||
| # License GPL3 | ||
|
|
||
|
|
||
| getv <- function(x, a) { | ||
| cumsum(values(aggregate(x, a, sum, na.rm=TRUE)) / 1000) | ||
| } | ||
|
|
||
| splitNS3 <- function(x) { | ||
| v <- getv(x, c(1, ncol(x))) | ||
| m1 <- which.min(abs(v - max(v)/3)) | ||
| m2 <- which.min(abs(v - max(v) * 2/3)) | ||
| list(n=x[1:m1, ,drop=FALSE], m=x[(m1+1):m2, ,drop=FALSE], s=x[(m2+1):nrow(x), ,drop=FALSE]) | ||
| } | ||
|
|
||
|
|
||
| splitNS <- function(x) { | ||
| v <- getv(x, c(1, ncol(x))) | ||
| m <- which.min(abs(v - max(v)/2)) | ||
| list(n=x[1:m, ,drop=FALSE], s=x[(m+1):nrow(x), ,drop=FALSE]) | ||
| } | ||
|
|
||
| splitWE3 <- function(x) { | ||
| v <- getv(x, c(nrow(x), 1)) | ||
| m1 <- which.min(abs(v - max(v)/3)) | ||
| m2 <- which.min(abs(v - max(v) * 2/3)) | ||
| list(w=x[, 1:m1, drop=FALSE], m=x[, (m1+1):m2, drop=FALSE], e=x[, (m2+1):ncol(x), drop=FALSE]) | ||
| } | ||
|
|
||
|
|
||
| splitWE <- function(x) { | ||
| v <- getv(x, c(nrow(x), 1)) | ||
| m <- which.min(abs(v - max(v)/2)) | ||
| list(w=x[, 1:m, drop=FALSE], e=x[, (m+1):ncol(x), drop=FALSE]) | ||
| } | ||
|
|
||
|
|
||
| setMethod("divide", signature(x="SpatRaster"), | ||
| function(x, n=2, start="ns", as.raster=FALSE, na.rm=TRUE) { | ||
|
|
||
| # if (!is.null(border)) stopifnot(inherits(border, "SpatVector")) | ||
|
|
||
| if (nlyr(x) > 1) { | ||
| warn("divide", "only the first layer is used") | ||
| x <- x[[1]] | ||
| } | ||
|
|
||
| n <- round(n) | ||
| if (length(n) > 1) { | ||
| if (!all(n %in% c(-2,-1,1,2))) { | ||
| error("divide", "if (length(n) > 1), values can only be -3, -2 for WE and 2, 3 for NS") | ||
| } | ||
| } else { | ||
| if (!isTRUE(n > 0)) { | ||
| error("divide", "n must be > 0") | ||
| } | ||
| start <- match.arg(tolower(start), c("ns", "ew")) | ||
| north <- start == "ns" | ||
| } | ||
|
|
||
| if (!hasValues(x)) { | ||
| out <- x | ||
| if (length(n) > 1) { | ||
| nrow(out) <- sum(abs(n[n<0]) + 1) | ||
| ncol(out) <- sum(abs(n[n<0]) + 1) | ||
| } else { | ||
| if (north) { | ||
| rows <- ceiling(n / 2) | ||
| cols <- n - rows | ||
| } else { | ||
| cols <- ceiling(n / 2) | ||
| rows <- n - cols | ||
| } | ||
| nrow(out) <- rows * 2 | ||
| ncol(out) <- cols * 2 | ||
| } | ||
| out <- as.polygons(out) | ||
| } else { | ||
|
|
||
| x <- classify(trim(x), cbind(NA, 0)) | ||
| out <- list(x) | ||
| if (length(n) > 1) { | ||
| for (i in 1:length(n)) { | ||
| if (n[i] == 1) { | ||
| out <- unlist(lapply(out, \(i) splitNS(i))) | ||
| } else if (n[i] == -1) { | ||
| out <- unlist(lapply(out, \(i) splitWE(i))) | ||
| } else if (n[i] == 2) { | ||
| out <- unlist(lapply(out, \(i) splitNS3(i))) | ||
| } else { # if (n[i] == -2) { | ||
| out <- unlist(lapply(out, \(i) splitWE3(i))) | ||
| } | ||
| } | ||
| } else { | ||
| for (i in 1:n) { | ||
| if (north) { | ||
| out <- unlist(lapply(out, \(s) splitNS(s))) | ||
| } else { | ||
| out <- unlist(lapply(out, \(s) splitWE(s))) | ||
| } | ||
| north <- !north | ||
| } | ||
| } | ||
| out <- lapply(out, \(i) as.polygons(ext(i))) |> vect() | ||
| crs(out) <- crs(out) | ||
| } | ||
| out$zones <- 1:nrow(out) | ||
|
|
||
| if (na.rm) { | ||
| border <- as.polygons(not.na(x), TRUE) | ||
| out <- crop(out, border) | ||
| } | ||
| if (isTRUE(as.raster) || is.na(as.raster)) { | ||
| r <- rasterize(out, x, "zones") | ||
| if (is.na(as.raster)) { | ||
| return(list(r=r, v=out)) | ||
| } else { | ||
| return(r) | ||
| } | ||
| } else { | ||
| out | ||
| } | ||
| } | ||
| ) | ||
|
|
||
|
|
||
|
|
||
| check_frac <- function(f) { | ||
| if (is.null(f)) return(f) | ||
| if ((any(f <= 0)) || (any(f >= 1))) { | ||
| stop("f values must be > 0 and < 1") | ||
| } | ||
| if (length(f) == 1) { | ||
| f <- seq(f, 1, f) | ||
| f[length(f)] <- 1 | ||
| } else { | ||
| if (length(unique(f)) < length(f)) { | ||
| stop("f values must be unique") | ||
| } | ||
| ord <- order(f) | ||
| if (!all(ord == 1:length(ord))) { | ||
| stop("f values must be in ascending order") | ||
| } | ||
| if (!isTRUE(all.equal(sum(f), 1, tolerance=0.0001))) { | ||
| stop("f values must sum to 1") | ||
| } | ||
| } | ||
| f | ||
| } | ||
|
|
||
|
|
||
|
|
||
| strip_polygon <- function(x, vertical, horizontal) { | ||
|
|
||
| totalArea <- expanse(x, transform=FALSE, unit="km") | ||
| e <- ext(x) | ||
| ex <- data.frame(t(as.vector(e + 1))) | ||
| e <- data.frame(t(as.vector(e))) | ||
|
|
||
| if (!is.null(vertical)) { | ||
| edges <- sapply(vertical, function(fraction){ | ||
| target <- totalArea * fraction | ||
| target_fun <- function(xm){ | ||
| expanse(crop(x, ext(ex$xmin, xm, ex$ymin, ex$ymax)), transform=FALSE, unit="km") - target | ||
| } | ||
| stats::uniroot(target_fun, lower=e$xmin+0.0000001, upper=e$xmax)$root | ||
| }) | ||
| xbnds <- matrix(c(ex$xmin, rep(edges,rep(2,length(edges))), ex$xmax), ncol=2, byrow=TRUE) | ||
| xbnds <- cbind(xbnds, ex$ymin, ex$ymax) | ||
| xbnds <- do.call(rbind, apply(xbnds, 1, \(i) as.polygons(ext(i)))) | ||
| xbnds$vid <- 1:nrow(xbnds) | ||
| } | ||
| if (!is.null(horizontal)) { | ||
| edges <- sapply(horizontal, function(fraction){ | ||
| target <- totalArea * fraction | ||
| target_fun <- function(ym){ | ||
| expanse(crop(x, ext(ex$xmin, ex$xmax, ex$ymin, ym)), transform=FALSE, unit="km") - target | ||
| } | ||
| stats::uniroot(target_fun, lower=e$ymin+0.0000001, upper=e$ymax)$root | ||
| }) | ||
| ybnds <- matrix(c(ex$ymin, rep(edges,rep(2,length(edges))), ex$ymax), ncol=2, byrow=TRUE) | ||
| ybnds <- cbind(ex$xmin, ex$xmax, ybnds) | ||
| ybnds <- do.call(rbind, apply(ybnds, 1, \(i) as.polygons(ext(i)))) | ||
| ybnds$hid <- 1:nrow(ybnds) | ||
|
|
||
| if (!is.null(vertical)) { | ||
| bnds <- union(xbnds, ybnds) | ||
| } else { | ||
| bnds <- ybnds | ||
| } | ||
| } | ||
| intersect(x, bnds) | ||
| } | ||
|
|
||
|
|
||
| divide_polygon <- function(x, n, w, alpha, ...) { | ||
| s <- terra::spatSample(x, max(n*4, 1000, log(n) * 100), "regular") | ||
| xy <- terra::crds(s) | ||
| if (!is.null(w)) { | ||
| e <- extract(w, s, ID=FALSE) | ||
| alpha <- rep_len(alpha, 2) | ||
| xy[,1] <- xy[,1] * alpha[1] | ||
| xy[,2] <- xy[,2] * alpha[2] | ||
| xy <- na.omit(cbind(xy, e)) | ||
| } | ||
| k <- stats::kmeans(xy, centers = n, ...) | ||
| ctrs <- k$centers[, 1:2] | ||
| if (!is.null(w)) { | ||
| ctrs[,1] <- ctrs[,1] / alpha[1] | ||
| ctrs[,2] <- ctrs[,2] / alpha[2] | ||
| } | ||
| v <- terra::voronoi(vect(ctrs, crs=xcrs), bnd=x) | ||
| terra::crop(v, x) | ||
| } | ||
|
|
||
|
|
||
| setMethod("divide", signature(x="SpatVector"), | ||
| function(x, n, w=NULL, alpha=1, ...) { | ||
| if (geomtype(x) == "polygons") { | ||
| error("divide", "the geometry type must be polgyons") | ||
| } | ||
| if (is.list(n)) { | ||
| vertical <- check_frac(n$v) | ||
| horizontal <- check_frac(n$h) | ||
| if (is.null(vertical) && is.null(horizontal)) return(x) | ||
| out <- lapply(1:nrow(x), function(i) strip_polygon(x[i], vertical, horizontal)) | ||
| } else { | ||
| n <- round(n) | ||
| stopifnot(n > 0) | ||
| if (n == 1) return(deepcopy(x)) | ||
| xcrs <- crs(x) | ||
| crs(x) <- "+proj=utm +zone=1" | ||
| out <- lapply(1:nrow(x), function(i) divide_polygon(x[i], n, w, alpha, ...)) | ||
| } | ||
| return(do.call(rbind, out)) | ||
| } | ||
| ) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,62 @@ | ||
| \name{divide} | ||
|
|
||
| \alias{divide} | ||
|
|
||
| \title{ | ||
| Subdivide a raster or polygons | ||
| } | ||
|
|
||
| \description{ | ||
| Divide a \code{SpatRaster} into \code{n} parts with approximately the same sum of weights (cell values). | ||
|
|
||
| Divides a \code{SpatVector} of polygons into \code{n} compact and approximately equal area parts. The results are not deterministic so you should use set.seed to be able to reproduce your results. If you get a warning about non-convergence, you can increase the number of iterations used with additional argument iter.max | ||
| } | ||
|
|
||
| \usage{ | ||
| \S4method{divide}{SpatRaster}(x, n=2, start="ns", as.raster=FALSE, na.rm=TRUE) | ||
|
|
||
| \S4method{divide}{SpatVector}(x, n=5, w=NULL, alpha=1) | ||
| } | ||
|
|
||
| \arguments{ | ||
| \item{x}{SpatRaster or SpatVector of polygons} | ||
| \item{n}{numeric. Can be a single positive integer to indicate the number of parts (SpatVector) or the number of splits (SpatRaster). | ||
|
|
||
| If \code{x} is a SpatRaster, it can also be a vector with values -2, -1, 1, or 2. Where 1 means one split and 2 means two splits, and the negative sign indicates an East-West (vertical) split as opposed to a North-South split. | ||
|
|
||
| If \code{x} is a SpatVector it can be a list with at least one of these elements: \code{horizontal} and \code{vertical} that specify the proportions of the area that splits should cover. This can either be a single fraction such as 1/3, or a sequence of increasing fractions such as \code{c(1/4, 1/2, 1)}} | ||
|
|
||
| \item{start}{character. To indicate the initial direction of splitting the raster. "ns" for North-South (horizontal) or "ew" for East-West (vertical)} | ||
| \item{as.raster}{logical. If \code{FALSE} a SpatVector is returned. If \code{FALSE}, a SpatRaster is returned. If \code{NA} a list with a SpatRaster and a SpatVector is returned} | ||
| \item{na.rm}{logical. If \code{TRUE} cells in \code{x} that are \code{NA} are not included in the output} | ||
|
|
||
| \item{w}{SpatRaster with, for example, environmental data} | ||
| \item{alpha}{numeric. One or two numbers that act as weights for the x and y coordinates} | ||
| \item{...}{additional arguments such as \code{iter.max} passed on to \code{\link{kmeans}}} | ||
|
|
||
| \item{f}{numeric vector of fractions. These must be > 0 and < 1, and in ascending order} | ||
| \item{vertical}{logical. If \code{TRUE} the strips are vertical} | ||
|
|
||
| } | ||
|
|
||
|
|
||
| \value{ | ||
| SpatVector or SpatRaster, or a list with bith | ||
| } | ||
|
|
||
|
|
||
| \examples{ | ||
| f <- system.file("ex/elev.tif", package="terra") | ||
| r <- rast(f) | ||
| x <- divide(r, 3) | ||
| # plot(r); lines(x) | ||
|
|
||
|
|
||
| f <- system.file("ex/lux.shp", package="terra") | ||
| v <- vect(f) | ||
| d <- divide(v, 3) | ||
| dv <- divide(v, list(h=.5)) | ||
| } | ||
|
|
||
| \keyword{spatial} | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters