fixes #1419

rspatial · Feb 6, 2024 · c3a80f2 · c3a80f2
1 parent d3b5220
commit c3a80f2
Show file tree

Hide file tree

Showing 12 changed files with 61 additions and 55 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,8 @@
 # version 1.7-72
 
+## enhancements
+
+- `extract` has new argument "small=TRUE" to allow for strict use of "touches=FALSE" [#1419](https://github.com/rspatial/terra/issues/1419) by Floris Vanderhaeghe.
 
 
 # version 1.7-71

diff --git a/R/cells.R b/R/cells.R
@@ -29,10 +29,10 @@ setMethod("cells", signature(x="SpatRaster", y="numeric"),
 
 
 setMethod("cells", signature("SpatRaster", "SpatVector"),
-	function(x, y, method="simple", weights=FALSE, exact=FALSE, touches=is.lines(y)) {
+	function(x, y, method="simple", weights=FALSE, exact=FALSE, touches=is.lines(y), small=TRUE) {
 		method = match.arg(tolower(method), c("simple", "bilinear"))
 		opt <- spatOptions()
-		d <- x@ptr$vectCells(y@ptr, touches[1], method[1], weights[1], exact[1], opt)
+		d <- x@ptr$vectCells(y@ptr, touches[1], small[1], method[1], weights[1], exact[1], opt)
 		if (geomtype(y) == "points") {
 			d <- matrix(d, nrow=nrow(y), byrow=TRUE)
 			d <- cbind(1:nrow(y), d)

diff --git a/R/extract.R b/R/extract.R
@@ -80,7 +80,7 @@ use_layer <- function(e, y, layer, nl) {
 }
 
 
-extract_table <- function(x, y, ID=FALSE, weights=FALSE, exact=FALSE, touches=FALSE, na.rm=FALSE) {
+extract_table <- function(x, y, ID=FALSE, weights=FALSE, exact=FALSE, touches=FALSE, small=TRUE, na.rm=FALSE) {
 
 	if (weights && exact) {
 		exact = FALSE
@@ -107,7 +107,7 @@ extract_table <- function(x, y, ID=FALSE, weights=FALSE, exact=FALSE, touches=FA
 			)
 		}
 
-		e <- x@ptr$extractVector(y@ptr, touches[1], "simple", FALSE, FALSE, 
+		e <- x@ptr$extractVector(y@ptr, touches[1], small[1], "simple", FALSE, FALSE, 
 			isTRUE(weights[1]), isTRUE(exact[1]), opt)
 		x <- messages(x, "extract")
 		e <- lapply(e, wtable, na.rm=na.rm)
@@ -131,7 +131,7 @@ extract_table <- function(x, y, ID=FALSE, weights=FALSE, exact=FALSE, touches=FA
 		}
 		if (nlyr(x) == 1) return(out[[1]]) else return(out)
 	} else {
-		e <- x@ptr$extractVectorFlat(y@ptr, "", FALSE, touches[1], "", FALSE, FALSE, FALSE, FALSE, opt)
+		e <- x@ptr$extractVectorFlat(y@ptr, "", FALSE, touches[1], small[1], "", FALSE, FALSE, FALSE, FALSE, opt)
 		x <- messages(x, "extract")
 		e <- data.frame(matrix(e, ncol=nlyr(x)+1, byrow=TRUE))
 		colnames(e) <- c("ID", names(x))
@@ -169,7 +169,7 @@ extract_fun <- function(x, y, fun, ID=TRUE, weights=FALSE, exact=FALSE, touches=
 	}
 
 	opt <- spatOptions()
-	e <- x@ptr$extractVectorFlat(y@ptr, fun, na.rm, touches[1], "", FALSE, FALSE, weights, exact, opt)
+	e <- x@ptr$extractVectorFlat(y@ptr, fun, na.rm, touches[1], small[1], "", FALSE, FALSE, weights, exact, opt)
 	x <- messages(x, "extract")
 	e <- data.frame(matrix(e, ncol=nl*nf, byrow=TRUE))
 	if (nf == 1) {
@@ -220,7 +220,7 @@ do_fun <- function(x, e, fun, ...) {
 
 
 setMethod("extract", signature(x="SpatRaster", y="SpatVector"),
-function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weights=FALSE, exact=FALSE, touches=is.lines(y), layer=NULL, bind=FALSE, raw=FALSE, ...) {
+function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weights=FALSE, exact=FALSE, touches=is.lines(y), small=TRUE, layer=NULL, bind=FALSE, raw=FALSE, ...) {
 
 	geo <- geomtype(y)
 	if (geo == "points") {		
@@ -240,9 +240,9 @@ function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weight
 				if (!is.null(layer)) {
 					warn("extract", "argument 'layer' is ignored when 'fun=table'")
 				}
-				e <- extract_table(x, y, ID=ID, weights=weights, exact=exact, touches=touches, ...)
+				e <- extract_table(x, y, ID=ID, weights=weights, exact=exact, touches=touches, small=small, ...)
 			} else {
-				e <- extract_fun(x, y, txtfun, ID=ID, weights=weights, exact=exact, touches=touches, bind=bind, layer=layer, ...)
+				e <- extract_fun(x, y, txtfun, ID=ID, weights=weights, exact=exact, touches=touches, small=small, bind=bind, layer=layer, ...)
 			}
 			return(e)
 		} else if (weights || exact) {
@@ -253,7 +253,7 @@ function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weight
 	} 
 
 	opt <- spatOptions()
-	e <- x@ptr$extractVectorFlat(y@ptr, "", FALSE, touches[1], method, isTRUE(cells[1]), isTRUE(xy[1]), isTRUE(weights[1]), isTRUE(exact[1]), opt)
+	e <- x@ptr$extractVectorFlat(y@ptr, "", FALSE, touches[1], small[1], method, isTRUE(cells[1]), isTRUE(xy[1]), isTRUE(weights[1]), isTRUE(exact[1]), opt)
 	x <- messages(x, "extract")
 
 	cn <- c("ID", names(x))

diff --git a/R/zonal.R b/R/zonal.R
@@ -162,7 +162,7 @@ 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, wide=TRUE, filename="", wopt=list())  {
+	function(x, z, fun="mean", na.rm=FALSE, w=NULL, weights=FALSE, exact=FALSE, touches=FALSE, small=TRUE, as.raster=FALSE, as.polygons=FALSE, wide=TRUE, filename="", wopt=list())  {
 		opt <- spatOptions()
 		txtfun <- .makeTextFun(fun)
 		if (!inherits(txtfun, "character")) {
@@ -172,7 +172,7 @@ setMethod("zonal", signature(x="SpatRaster", z="SpatVector"),
 				if (!is.null(w)) {
 					error("cannot use 'w' when 'fun=table'")
 				}
-				v <- x@ptr$zonal_poly_table(z@ptr, weights[1], exact[1], touches[1], na.rm, opt)
+				v <- x@ptr$zonal_poly_table(z@ptr, weights[1], exact[1], touches[1], small[1], na.rm, opt)
 				messages(x, "zonal")
 
 				v <- lapply(v, function(i) if (length(i) == 0) NA else i)
@@ -201,12 +201,12 @@ setMethod("zonal", signature(x="SpatRaster", z="SpatVector"),
 				return(v)
 			} else {
 				if (is.null(w)) {
-					out <- x@ptr$zonal_poly(z@ptr, txtfun, weights[1], exact[1], touches[1], na.rm, opt)
+					out <- x@ptr$zonal_poly(z@ptr, txtfun, weights[1], exact[1], touches[1], small[1], na.rm, opt)
 				} else {
 					if (txtfun != "mean") {
 						error("zonal", "fun must be 'mean' when using weights")
 					}
-					out <- x@ptr$zonal_poly_weighted(z@ptr, w@ptr, weights[1], exact[1], touches[1], na.rm, opt)
+					out <- x@ptr$zonal_poly_weighted(z@ptr, w@ptr, weights[1], exact[1], touches[1], small[1], na.rm, opt)
 				}
 				messages(out, "zonal")
 				out <- .getSpatDF(out)

diff --git a/man/cells.Rd b/man/cells.Rd
@@ -33,6 +33,7 @@ Get the cell numbers covered by a SpatVector or SpatExtent. Or that match values
   \item{weights}{logical. If \code{TRUE} and \code{y} has polygons, the approximate fraction of each cell that is covered is returned as well}
   \item{exact}{logical. If \code{TRUE} and \code{y} has polygons, the exact fraction of each cell that is covered is returned as well}
   \item{touches}{logical. If \code{TRUE}, values for all cells touched by lines or polygons are extracted, not just those on the line render path, or whose center point is within the polygon. Not relevant for points}
+  \item{small}{logical. If \code{TRUE}, values for all cells in touched polygons are extracted if none of the cells center points is within the polygon; even if \code{touches=FALSE}}
 }
 
 \value{

diff --git a/man/extract.Rd b/man/extract.Rd
@@ -28,7 +28,7 @@ Alternatively, you can use \code{\link{zonal}} after using \code{\link{rasterize
 
 \usage{
 \S4method{extract}{SpatRaster,SpatVector}(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE,
-    ID=TRUE, weights=FALSE, exact=FALSE, touches=is.lines(y),
+    ID=TRUE, weights=FALSE, exact=FALSE, touches=is.lines(y), small=TRUE,
 	layer=NULL, bind=FALSE, raw=FALSE, ...)
 
 \S4method{extract}{SpatRaster,SpatExtent}(x, y, cells=FALSE, xy=FALSE)
@@ -51,6 +51,7 @@ Alternatively, you can use \code{\link{zonal}} after using \code{\link{rasterize
 \item{weights}{logical. If \code{TRUE} and \code{y} has polygons, the approximate fraction of each cell that is covered is returned as well, for example to compute a weighted mean}
 \item{exact}{logical. If \code{TRUE} and \code{y} has polygons, the exact fraction of each cell that is covered is returned as well, for example to compute a weighted mean}
 \item{touches}{logical. If \code{TRUE}, values for all cells touched by lines or polygons are extracted, not just those on the line render path, or whose center point is within the polygon. Not relevant for points; and always considered \code{TRUE} when \code{weights=TRUE} or \code{exact=TRUE}}
+\item{small}{logical. If \code{TRUE}, values for all cells in touched polygons are extracted if none of the cells center points is within the polygon; even if \code{touches=FALSE}}
 \item{layer}{character or numeric to select the layer to extract from for each geometry. If \code{layer} is a character it can be a name in \code{y} or a vector of layer names. If it is numeric, it must be integer values between \code{1} and \code{nlyr(x)}}
 \item{bind}{logical. If \code{TRUE}, a SpatVector is returned consisting of the input SpatVector \code{y} and the \code{cbind}-ed extracted values}
 \item{raw}{logical. If \code{TRUE}, a matrix is returned with the "raw" numeric cell values. If \code{FALSE}, a data.frame is returned and the cell values are transformed to factor, logical, or integer values, where appropriate}

diff --git a/man/zonal.Rd b/man/zonal.Rd
@@ -21,7 +21,7 @@ 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.polygons=FALSE, wide=TRUE, 
+		exact=FALSE, touches=FALSE, small=TRUE, as.raster=FALSE, as.polygons=FALSE, wide=TRUE, 
 		filename="", wopt=list())
 
 \S4method{zonal}{SpatVector,SpatVector}(x, z, fun=mean, ..., weighted=FALSE, as.polygons=FALSE) 
@@ -42,7 +42,7 @@ You can also summarize values of a SpatVector for each polygon (zone) defined by
   \item{weights}{logical. If \code{TRUE} and \code{y} has polygons, the approximate fraction of each cell that is covered is returned as well, for example to compute a weighted mean}
   \item{exact}{logical. If \code{TRUE} and \code{y} has polygons, the exact fraction of each cell that is covered is returned as well, for example to compute a weighted mean}
   \item{touches}{logical. If \code{TRUE}, values for all cells touched by lines or polygons are extracted, not just those on the line render path, or whose center point is within the polygon. Not relevant for points; and always considered \code{TRUE} when \code{weights=TRUE} or \code{exact=TRUE}}
-
+  \item{small}{logical. If \code{TRUE}, values for all cells in touched polygons are extracted if none of the cells center points is within the polygon; even if \code{touches=FALSE}}
   \item{weighted}{logical. If \code{TRUE}, a weighted.mean is computed and \code{fun} is ignored. Weights are based on the length of the lines or the area of the polygons in \code{x} that intersect with \code{z}. This argument is ignored of \code{x} is a SpatVector or points}  
   \item{as.polygons}{logical. Should the zonal statistics be combined with the geometry of \code{z}?}
   \item{na.rm}{logical. If \code{TRUE}, \code{NA}s are removed}

diff --git a/src/extract.cpp b/src/extract.cpp
@@ -645,7 +645,7 @@ std::vector<double> SpatRaster::extractXYFlat(const std::vector<double> &x, cons
 
 
 // <geom<layer<values>>>
-std::vector<std::vector<std::vector<double>>> SpatRaster::extractVector(SpatVector v, bool touches, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
+std::vector<std::vector<std::vector<double>>> SpatRaster::extractVector(SpatVector v, bool touches, bool small, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
 
 	if (!source[0].srs.is_same(v.srs, true)) {
 		v = v.project(getSRS("wkt"), false);
@@ -736,7 +736,7 @@ std::vector<std::vector<std::vector<double>>> SpatRaster::extractVector(SpatVect
 					rasterizeCellsExact(cell, wgt, p, opt);
 				}
 			} else {
-				cell = rasterizeCells(p, touches, opt);
+				cell = rasterizeCells(p, touches, small, opt);
             }
 			srcout = extractCell(cell);
             for (size_t j=0; j<nl; j++) {
@@ -898,7 +898,7 @@ std::vector<double> SpatRaster::extractVectorFlat(SpatVector v, std::string fun,
 */
 
 
-std::vector<double> SpatRaster::extractVectorFlat(SpatVector v, std::vector<std::string> funs, bool narm, bool touches, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
+std::vector<double> SpatRaster::extractVectorFlat(SpatVector v, std::vector<std::string> funs, bool narm, bool touches, bool small, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
 
 	if (!source[0].srs.is_same(v.srs, true)) {
 		v = v.project(getSRS("wkt"), false);
@@ -1027,7 +1027,7 @@ std::vector<double> SpatRaster::extractVectorFlat(SpatVector v, std::vector<std:
 				rasterizeCellsExact(cell, wgt, p, opt);
 			}
 		} else {
-			cell = rasterizeCells(p, touches, opt);
+			cell = rasterizeCells(p, touches, small, opt);
 		}
 
 		if (havefun) {
@@ -1240,7 +1240,7 @@ std::vector<double> SpatRaster::extractCellFlat(std::vector<double> &cell) {
 }
 
 
-std::vector<double> SpatRaster::vectCells(SpatVector v, bool touches, std::string method, bool weights, bool exact, SpatOptions &opt) {
+std::vector<double> SpatRaster::vectCells(SpatVector v, bool touches, bool small, std::string method, bool weights, bool exact, SpatOptions &opt) {
 
 	std::string gtype = v.type();
 	if (gtype != "polygons") weights = false;
@@ -1279,7 +1279,7 @@ std::vector<double> SpatRaster::vectCells(SpatVector v, bool touches, std::strin
 				cells.insert(cells.end(), cnr.begin(), cnr.end());
 				wghts.insert(wghts.end(), wght.begin(), wght.end());
 			} else {
-				std::vector<double> geomc = rasterizeCells(p, touches, opt);
+				std::vector<double> geomc = rasterizeCells(p, touches, small, opt);
 				std::vector<double> id(geomc.size(), i);
 				out.insert(out.end(), id.begin(), id.end());
 				cells.insert(cells.end(), geomc.begin(), geomc.end());
@@ -1350,12 +1350,12 @@ std::vector<std::vector<std::vector<double>>> SpatRasterStack::extractCell(std::
 
 
 // this is rather inefficient (repeated rasterization)
-std::vector<std::vector<std::vector<std::vector<double>>>> SpatRasterStack::extractVector(SpatVector v, bool touches, std::string method, SpatOptions &opt) {
+std::vector<std::vector<std::vector<std::vector<double>>>> SpatRasterStack::extractVector(SpatVector v, bool touches, bool small, std::string method, SpatOptions &opt) {
 	unsigned ns = nsds();
 	std::vector<std::vector<std::vector<std::vector<double>>>> out(ns);
 	for (size_t i=0; i<ns; i++) {
 		SpatRaster r = getsds(i);
-		out[i] = r.extractVector(v, touches, method, false, false, false, false, opt);
+		out[i] = r.extractVector(v, touches, small, method, false, false, false, false, opt);
 	}
 	return out;
 }

diff --git a/src/raster_stats.cpp b/src/raster_stats.cpp
@@ -1283,7 +1283,7 @@ SpatDataFrame SpatRaster::zonal_weighted(SpatRaster z, SpatRaster w, bool narm,
 }
 
 
-SpatDataFrame SpatRaster::zonal_poly(SpatVector x, std::string fun, bool weights, bool exact, bool touches,bool narm, SpatOptions &opt) {
+SpatDataFrame SpatRaster::zonal_poly(SpatVector x, std::string fun, bool weights, bool exact, bool touches,bool small, bool narm, SpatOptions &opt) {
 
 	SpatDataFrame out;
 	std::string gtype = x.type();
@@ -1330,7 +1330,7 @@ SpatDataFrame SpatRaster::zonal_poly(SpatVector x, std::string fun, bool weights
 		} else if (exact) {
 			rasterizeCellsExact(cell, wgt, p, opt);
 		} else {
-			cell = rasterizeCells(p, touches, opt);
+			cell = rasterizeCells(p, touches, small, opt);
         }
 
 		std::vector<std::vector<double>> e = extractCell(cell);
@@ -1383,7 +1383,7 @@ std::vector<double> tabfun(std::vector<double> x, std::vector<double> w) {
 }
 
 
-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>> SpatRaster::zonal_poly_table(SpatVector x, bool weights, bool exact, bool touches, bool small, bool narm, SpatOptions &opt) {
 
 	std::vector<std::vector<double>> out;
 	std::string gtype = x.type();
@@ -1401,7 +1401,7 @@ std::vector<std::vector<double>> SpatRaster::zonal_poly_table(SpatVector x, bool
 	if (nl > 1) {
 		SpatOptions ops(opt);
 		SpatRaster r = subset({0}, ops);
-		out = r.zonal_poly_table(x, weights, exact, touches, narm, opt);
+		out = r.zonal_poly_table(x, weights, exact, touches, small, narm, opt);
 		addWarning("only the first layer of the raster is used");		
 		return out;
 	}
@@ -1420,7 +1420,7 @@ std::vector<std::vector<double>> SpatRaster::zonal_poly_table(SpatVector x, bool
 		} else if (exact) {
 			rasterizeCellsExact(cell, wgt, p, opt);
 		} else {
-			cell = rasterizeCells(p, touches, opt);
+			cell = rasterizeCells(p, touches, small, opt);
         }	
 		std::vector<std::vector<double>> e = extractCell(cell);
 		out[i] = tabfun(e[0], wgt);
@@ -1430,7 +1430,7 @@ std::vector<std::vector<double>> SpatRaster::zonal_poly_table(SpatVector x, bool
 }
 
 
-SpatDataFrame SpatRaster::zonal_poly_weighted(SpatVector x, SpatRaster w, bool weights, bool exact, bool touches, bool narm, SpatOptions &opt) {
+SpatDataFrame SpatRaster::zonal_poly_weighted(SpatVector x, SpatRaster w, bool weights, bool exact, bool touches, bool small, bool narm, SpatOptions &opt) {
 
 	SpatDataFrame out;
 	std::string gtype = x.type();
@@ -1468,7 +1468,7 @@ SpatDataFrame SpatRaster::zonal_poly_weighted(SpatVector x, SpatRaster w, bool w
 		} else if (exact) {
 			rasterizeCellsExact(cell, wgt, p, opt);
 		} else {
-			cell = rasterizeCells(p, touches, opt);
+			cell = rasterizeCells(p, touches, small, opt);
         }
 
 		std::vector<std::vector<double>> e = extractCell(cell);

diff --git a/src/rasterize.cpp b/src/rasterize.cpp
@@ -762,7 +762,7 @@ SpatRaster SpatRaster::rasterize(SpatVector x, std::string field, std::vector<do
 }
 
 
-std::vector<double> SpatRaster::rasterizeCells(SpatVector &v, bool touches, SpatOptions &opt) {
+std::vector<double> SpatRaster::rasterizeCells(SpatVector &v, bool touches, bool small, SpatOptions &opt) {
 // note that this is only for lines and polygons
     SpatOptions ropt(opt);
 	SpatRaster r = geometry(1);
@@ -791,17 +791,21 @@ std::vector<double> SpatRaster::rasterizeCells(SpatVector &v, bool touches, Spat
 	SpatVector pts = rcr.as_points(false, true, false, ropt);
 	std::vector<double> cells;
 	if (pts.empty()) {
-		pts = v.as_points(false, true);
-		SpatDataFrame vd = pts.getGeometryDF();
-		std::vector<double> x = vd.getD(0);
-		std::vector<double> y = vd.getD(1);
-		cells = r.cellFromXY(x, y);	
-		cells.erase(std::remove_if(cells.begin(), cells.end(),
-                    [](const double& value) { return std::isnan(value); }), cells.end());
-		std::sort( cells.begin(), cells.end() );
-		cells.erase(std::unique(cells.begin(), cells.end()), cells.end());
-		if (cells.empty()) {
-			cells.resize(1, NAN);
+		if (small) {
+			pts = v.as_points(false, true);
+			SpatDataFrame vd = pts.getGeometryDF();
+			std::vector<double> x = vd.getD(0);
+			std::vector<double> y = vd.getD(1);
+			cells = r.cellFromXY(x, y);	
+			cells.erase(std::remove_if(cells.begin(), cells.end(),
+						[](const double& value) { return std::isnan(value); }), cells.end());
+			std::sort( cells.begin(), cells.end() );
+			cells.erase(std::unique(cells.begin(), cells.end()), cells.end());
+			if (cells.empty()) {
+				cells.resize(1, NAN);
+			}
+		} else {
+			cells.resize(1, NAN);			
 		}
 	} else {
 		SpatDataFrame vd = pts.getGeometryDF();