sr

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")
 

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)}
 }
 

src/extract.cpp

@@ -26,28 +26,40 @@
 
 
 
-std::vector<double> circ_dist(double xres, double yres, double d, std::vector<size_t> &dim) {
-	size_t nx = 1 + 2 * floor(d/xres);
-	size_t ny = 1 + 2 * floor(d/yres);
-
+std::vector<double> circ_dist(double xres, double yres, double d, std::vector<size_t> &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<double> out{1};
 		return out;
 	} 
 	dim = {ny, nx};
 
-	SpatRaster x({ny, nx, 1}, {0., nx * xres, 0., ny * yres}, "+proj=utm +zone=1 +datum=WGS84");
+	SpatRaster x({ny, nx, 1}, {0., nx * xres, 0., ny * yres}, crs);
 	std::vector<double> 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<double> out;
 	x.getValuesSource(0, out);				
 	out[out.size()/2] = 1;
@@ -79,7 +91,13 @@ std::vector<std::vector<double>> SpatRaster::extractBuffer(const std::vector<dou
 
 	std::vector<size_t> dim;
 	std::vector<double> cd;
-	std::vector<double> cb = circ_dist(xres(), yres(), b, dim);
+	double ymean = 0;
+	bool lonlat = is_lonlat();
+	if (lonlat) {
+		ymean = vmean(y, true);
+	}
+
+	std::vector<double> cb = circ_dist(xres(), yres(), b, dim, lonlat, ymean);
 	bool docb = false;
 	std::vector<bool> adj(cb.size(), false);
 	if (cb.size() > 1) {