Skip to content

Commit

Permalink
zonal r-v table
Browse files Browse the repository at this point in the history
  • Loading branch information
rhijmans committed Nov 10, 2023
1 parent 9494a5e commit 40e1ef5
Show file tree
Hide file tree
Showing 6 changed files with 107 additions and 12 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: terra
Type: Package
Title: Spatial Data Analysis
Version: 1.7-59
Date: 2023-11-01
Version: 1.7-60
Date: 2023-11-10
Depends: R (>= 3.5.0)
Suggests: parallel, tinytest, ncdf4, sf (>= 0.9-8), deldir, XML, leaflet, htmlwidgets
LinkingTo: Rcpp
Expand Down
48 changes: 40 additions & 8 deletions R/zonal.R
Original file line number Diff line number Diff line change
Expand Up @@ -162,22 +162,54 @@ setMethod("zonal", signature(x="SpatRaster", z="SpatRaster"),


setMethod("zonal", signature(x="SpatRaster", z="SpatVector"),
function(x, z, fun="mean", na.rm=FALSE, w=NULL, weights=FALSE, exact=FALSE, touches=FALSE, as.raster=FALSE, as.polygons=FALSE, filename="", wopt=list()) {
function(x, z, fun="mean", na.rm=FALSE, w=NULL, weights=FALSE, exact=FALSE, touches=FALSE, as.raster=FALSE, as.polygons=FALSE, wide=FALSE, filename="", wopt=list()) {
opt <- spatOptions()
txtfun <- .makeTextFun(fun)
if (!inherits(txtfun, "character")) {
error("zonal", "this 'fun' is not supported. You can use extract instead")
} else {
if (is.null(w)) {
out <- x@cpp$zonal_poly(z@cpp, txtfun, weights[1], exact[1], touches[1], na.rm, opt)
if (txtfun == "table") {
if (!is.null(w)) {
error("cannot use 'w' when 'fun=table'")
}
v <- x@cpp$zonal_poly_table(z@cpp, weights[1], exact[1], touches[1], na.rm, opt)
messages(x, "zonal")

v <- lapply(v, function(i) if (length(i) == 0) NA else i)
v <- lapply(1:length(v), function(i) cbind(i, matrix(v[[i]], ncol=2)))
v <- do.call(rbind, v)
v <- as.data.frame(v)
colnames(v) <- c("polygon", "value", "count")
ff <- is.factor(x)[1]
if (ff) {
cg <- cats(x)[[1]]
i <- match(v$value, cg[,1])
act <- activeCat(x, 1) + 1
v$value <- cg[i, act]
}
if (as.polygons | wide) {
nms <- names(v)
v <- reshape(v, direction="wide", idvar=nms[1], timevar=nms[2])
names(v) <- gsub("count.", "", names(v))
v[is.na(v)] <- 0
}
if (as.polygons) {
values(z) <- v
return(z)
}
return(v)
} else {
if (txtfun != "mean") {
error("zonal", "fun must be 'mean' when using weights")
if (is.null(w)) {
out <- x@cpp$zonal_poly(z@cpp, txtfun, weights[1], exact[1], touches[1], na.rm, opt)
} else {
if (txtfun != "mean") {
error("zonal", "fun must be 'mean' when using weights")
}
out <- x@cpp$zonal_poly_weighted(z@cpp, w@cpp, weights[1], exact[1], touches[1], na.rm, opt)
}
out <- x@cpp$zonal_poly_weighted(z@cpp, w@cpp, weights[1], exact[1], touches[1], na.rm, opt)
messages(out, "zonal")
out <- .getSpatDF(out)
}
messages(out, "zonal")
out <- .getSpatDF(out)
}
if (as.raster) {
if (is.null(wopt$names)) {
Expand Down
5 changes: 3 additions & 2 deletions man/zonal.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ You can also summarize values of a SpatVector for each polygon (zone) defined by
as.raster=FALSE, filename="", overwrite=FALSE, wopt=list())

\S4method{zonal}{SpatRaster,SpatVector}(x, z, fun="mean", na.rm=FALSE, w=NULL, weights=FALSE,
exact=FALSE, touches=FALSE, as.raster=FALSE, as.vector=FALSE, filename="", wopt=list())
exact=FALSE, touches=FALSE, as.raster=FALSE, as.vector=FALSE, wide=TRUE,
filename="", wopt=list())

\S4method{zonal}{SpatVector,SpatVector}(x, z, fun=mean, ..., weighted=FALSE, as.polygons=FALSE)
}
Expand All @@ -32,7 +33,7 @@ You can also summarize values of a SpatVector for each polygon (zone) defined by
\item{fun}{function to be applied to summarize the values by zone. Either as character: "mean", "min", "max", "sum", "isNA", and "notNA" and, for relatively small SpatRasters, a proper function}
\item{...}{additional arguments passed to fun, such as \code{na.rm=TRUE}}
\item{w}{SpatRaster with weights. Should have a single-layer with non-negative values}
\item{wide}{logical. Should the values returned in a wide format? This only affects the results when \code{nlyr(z) == 2}}
\item{wide}{logical. Should the values returned in a wide format? For the \code{SpatRaster, SpatRaster} method this only affects the results when \code{nlyr(z) == 2}. For the \code{SpatRaster, SpatVector} method this only affects the results when \code{fun=table}}
\item{as.raster}{logical. If \code{TRUE}, a SpatRaster is returned with the zonal statistic for each zone}
\item{filename}{character. Output filename (ignored if \code{as.raster=FALSE}}
\item{overwrite}{logical. If \code{TRUE}, \code{filename} is overwritten}
Expand Down
1 change: 1 addition & 0 deletions src/RcppModule.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -963,6 +963,7 @@ RCPP_MODULE(spat){
.method("zonal", &SpatRaster::zonal)
.method("zonal_weighted", &SpatRaster::zonal_weighted)
.method("zonal_poly", &SpatRaster::zonal_poly)
.method("zonal_poly_table", &SpatRaster::zonal_poly_table)
.method("zonal_poly_weighted", &SpatRaster::zonal_poly_weighted)
// .method("zonal_old", &SpatRaster::zonal_old)
;
Expand Down
57 changes: 57 additions & 0 deletions src/raster_stats.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1373,6 +1373,63 @@ SpatDataFrame SpatRaster::zonal_poly(SpatVector x, std::string fun, bool weights
}


std::vector<double> tabfun(std::vector<double> x, std::vector<double> w) {
// if (w.size() == 0) {
std::map<double, long long unsigned> tab = table(x);
return vtable(tab);
// } else {

// }
}


std::vector<std::vector<double>> SpatRaster::zonal_poly_table(SpatVector x, bool weights, bool exact, bool touches,bool narm, SpatOptions &opt) {

std::vector<std::vector<double>> out;
std::string gtype = x.type();
if (gtype != "polygons") {
setError("SpatVector must have polygon geometry");
return out;
}

if (!hasValues()) {
setError("raster has no values");
return out;
}

unsigned nl = nlyr();
if (nl > 1) {
SpatOptions ops(opt);
SpatRaster r = subset({0}, ops);
out = r.zonal_poly_table(x, weights, exact, touches, narm, opt);
addWarning("only the first layer of the raster is used");
return out;
}

unsigned ng = x.size();
std::vector<std::vector<double>> zv(nl, std::vector<double>(ng));
out.resize(ng);
SpatRaster r = geometry(1);
for (size_t i=0; i<ng; i++) {
SpatGeom g = x.getGeom(i);
SpatVector p(g);
p.srs = x.srs;
std::vector<double> cell, wgt;
if (weights) {
rasterizeCellsWeights(cell, wgt, p, opt);
} else if (exact) {
rasterizeCellsExact(cell, wgt, p, opt);
} else {
cell = rasterizeCells(p, touches, opt);
}
std::vector<std::vector<double>> e = extractCell(cell);
out[i] = tabfun(e[0], wgt);
}

return out;
}


SpatDataFrame SpatRaster::zonal_poly_weighted(SpatVector x, SpatRaster w, bool weights, bool exact, bool touches, bool narm, SpatOptions &opt) {

SpatDataFrame out;
Expand Down
4 changes: 4 additions & 0 deletions src/spatRaster.h
Original file line number Diff line number Diff line change
Expand Up @@ -834,6 +834,10 @@ class SpatRaster {

SpatDataFrame zonal_poly(SpatVector x, std::string fun, bool weights, bool exact, bool touches, bool narm, SpatOptions &opt);
SpatDataFrame zonal_poly_weighted(SpatVector x, SpatRaster w, bool weights, bool exact, bool touches, bool narm, SpatOptions &opt);

std::vector<std::vector<double>> zonal_poly_table(SpatVector x, bool weights, bool exact, bool touches,bool narm, SpatOptions &opt);



// SpatDataFrame zonal_old(SpatRaster x, std::string fun, bool narm, SpatOptions &opt);
SpatRaster rgb2col(size_t r, size_t g, size_t b, SpatOptions &opt);
Expand Down

0 comments on commit 40e1ef5

Please sign in to comment.