From 4c3ae951d649e7484479685b0413c0e12b95772b Mon Sep 17 00:00:00 2001 From: rhijmans Date: Mon, 30 Dec 2024 22:21:40 -0800 Subject: [PATCH] sr --- R/extract.R | 5 +++-- man/extract.Rd | 2 +- src/extract.cpp | 38 ++++++++++++++++++++++++++++---------- 3 files changed, 32 insertions(+), 13 deletions(-) diff --git a/R/extract.R b/R/extract.R index a8b0b0077..0d949f059 100644 --- a/R/extract.R +++ b/R/extract.R @@ -231,8 +231,8 @@ function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weight if (geo == "points") { if (search_radius > 0) { - xy <- crds(y) - e <- x@pntr$extractBuffer(xy[,1], xy[,2], searchradius) + pts <- crds(y) + e <- x@pntr$extractBuffer(pts[,1], pts[,2], search_radius) messages(x) e <- do.call(cbind, e) colnames(e) <- c(names(x)[1], "distance", "cell") @@ -280,6 +280,7 @@ function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weight } + opt <- spatOptions() e <- x@pntr$extractVectorFlat(y@pntr, "", FALSE, touches[1], small[1], method, isTRUE(cells[1]), isTRUE(xy[1]), isTRUE(weights[1]), isTRUE(exact[1]), opt) x <- messages(x, "extract") diff --git a/man/extract.Rd b/man/extract.Rd index 9d8d49e55..4012724ba 100644 --- a/man/extract.Rd +++ b/man/extract.Rd @@ -55,7 +55,7 @@ Alternatively, you can use \code{\link{zonal}} after using \code{\link{rasterize \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} -\item{search_radius}{positive number. A search-radius that is used when \code{y} has point geometry. If this value is larger than zero, it is the maximum distance used to find the a cell with a value that is nearest to the cell that the point falls in if that cell that has a missing (\code{NA}) value. The value of this nearest cell, the distance to the original cell, and the new cell number are returned} +\item{search_radius}{positive number. A search-radius that is used when \code{y} has point geometry. If this value is larger than zero, it is the maximum distance used to find the a cell with a value that is nearest to the cell that the point falls in if that cell that has a missing (\code{NA}) value. The value of this nearest cell, the distance to the original cell, and the new cell number are returned. The radius should be expressed in m if the data have lon/lat coordinates or in the distance unit of the crs in other cases (typically also m). For lon/lat data, the mean latitude of the points is used to compute the distances, so this may be inprecise for data with a large latitudinal range} \item{...}{additional arguments to \code{fun} if \code{y} is a SpatVector. For example \code{na.rm=TRUE}. Or arguments passed to the \code{SpatRaster,SpatVector} method if \code{y} is a matrix (such as the \code{method} and \code{cells} arguments)} } diff --git a/src/extract.cpp b/src/extract.cpp index 1ae5b3289..08f3bdfec 100644 --- a/src/extract.cpp +++ b/src/extract.cpp @@ -26,10 +26,26 @@ -std::vector circ_dist(double xres, double yres, double d, std::vector &dim) { - size_t nx = 1 + 2 * floor(d/xres); - size_t ny = 1 + 2 * floor(d/yres); - +std::vector circ_dist(double xres, double yres, double d, std::vector &dim, bool lonlat, double ymean) { + + size_t nx, ny; + std::string crs; + if (lonlat) { + deg2rad(ymean); + double xr = xres; + double yr = yres; + deg2rad(xr); + deg2rad(yr); + double dx = distance_cos(0, ymean, xr, ymean); + double dy = distance_cos(0, ymean-0.5*yr, 0, ymean+0.5*yr); + nx = 1 + 2 * floor(d/dx); + ny = 1 + 2 * floor(d/dy); + crs = "+proj=longlat"; + } else { + nx = 1 + 2 * floor(d/xres); + ny = 1 + 2 * floor(d/yres); + crs = "+proj=utm +zone=1"; + } if ((nx == 1) || (ny == 1)) { dim = {1, 1}; std::vector out{1}; @@ -37,17 +53,13 @@ std::vector circ_dist(double xres, double yres, double d, std::vector v(nx*ny, NAN); v[v.size()/2] = 1; SpatOptions opt; x.setValues(v, opt); - x = x.distance(NAN, NAN, false, "m", false, "cosine", opt); -// SpatRaster y = x.arith(d, "<=", false, true, opt); -// x = x.mask(y, false, NAN, NAN, opt); - std::vector out; x.getValuesSource(0, out); out[out.size()/2] = 1; @@ -79,7 +91,13 @@ std::vector> SpatRaster::extractBuffer(const std::vector dim; std::vector cd; - std::vector cb = circ_dist(xres(), yres(), b, dim); + double ymean = 0; + bool lonlat = is_lonlat(); + if (lonlat) { + ymean = vmean(y, true); + } + + std::vector cb = circ_dist(xres(), yres(), b, dim, lonlat, ymean); bool docb = false; std::vector adj(cb.size(), false); if (cb.size() > 1) {