diff --git a/R/distance.R b/R/distance.R index fb0fc944c..13a24c0a6 100644 --- a/R/distance.R +++ b/R/distance.R @@ -39,8 +39,9 @@ setMethod("buffer", signature(x="SpatRaster"), # #) -setMethod("distance", signature(x="SpatRaster", y="missing"), - function(x, y, target=NA, exclude=NULL, unit="m", method="haversine", values=FALSE, filename="", ...) { + +setMethod("distance", signature(x="SpatRaster", y="missing"), + function(x, y, target=NA, exclude=NULL, unit="m", method="haversine", maxdist=NA, values=FALSE, filename="", ...) { if (!(method %in% c("geo", "haversine", "cosine"))) { error("distance", "not a valid method. Should be one of: 'geo', 'haversine', 'cosine'") @@ -59,7 +60,7 @@ setMethod("distance", signature(x="SpatRaster", y="missing"), } else { exclude <- NA } - x@pntr <- x@pntr$rastDistance(target, exclude, keepNA, tolower(unit), TRUE, method, isTRUE(values), opt) + x@pntr <- x@pntr$rastDistance(target, exclude, keepNA, tolower(unit), TRUE, method, isTRUE(values), maxdist, opt) messages(x, "distance") } ) diff --git a/man/distance.Rd b/man/distance.Rd index 126cacc9e..4f569eb01 100644 --- a/man/distance.Rd +++ b/man/distance.Rd @@ -43,7 +43,8 @@ If \code{y} is missing, the distance between each points in \code{x} with all ot \usage{ -\S4method{distance}{SpatRaster,missing}(x, y, target=NA, exclude=NULL, unit="m", method="haversine", values=FALSE, filename="", ...) +\S4method{distance}{SpatRaster,missing}(x, y, target=NA, exclude=NULL, unit="m", method="haversine", + maxdist=NA, values=FALSE, filename="", ...) \S4method{distance}{SpatRaster,SpatVector}(x, y, unit="m", rasterize=FALSE, method="cosine", filename="", ...) @@ -54,7 +55,7 @@ If \code{y} is missing, the distance between each points in \code{x} with all ot \S4method{distance}{matrix,matrix}(x, y, lonlat, pairwise=FALSE, unit="m", method="geo") \S4method{distance}{matrix,missing}(x, y, lonlat, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, - unit="m", method="geo") + unit="m", method="geo") } \arguments{ @@ -64,7 +65,8 @@ If \code{y} is missing, the distance between each points in \code{x} with all ot \item{exclude}{numeric. The value of the cells that should not be considered for computing distances} \item{unit}{character. Can be either "m" or "km"} \item{method}{character. One of "geo", "cosine" or "haversine" (the latter cannot be used for distances to lines or polygons). With "geo" the most precise but slower method of Karney (2003) is used. The other two methods are faster but less precise} - \item{values}{logical, If \code{TRUE}, the value of the nearest non-target cell is returned instead of the distance to that cell} + \item{maxdist}{numeric. Distance above this values are not set to \code{NA}} + \item{values}{logical. If \code{TRUE}, the value of the nearest non-target cell is returned instead of the distance to that cell} \item{rasterize}{logical. If \code{TRUE} distance is computed from the cells covered by the geometries after rasterization. This can be much faster in some cases} \item{filename}{character. Output filename} \item{...}{additional arguments for writing files as in \code{\link{writeRaster}}} diff --git a/src/distRaster.cpp b/src/distRaster.cpp index 082debfb0..c639cd3df 100644 --- a/src/distRaster.cpp +++ b/src/distRaster.cpp @@ -244,7 +244,7 @@ void dist_only(std::vector &d, const std::vector& vx, const std: } -SpatRaster SpatRaster::distance_crds(std::vector& x, std::vector& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt) { +SpatRaster SpatRaster::distance_crds(std::vector& x, std::vector& y, const std::string& method, bool skip, bool setNA, std::string unit, double max_dist, SpatOptions &opt) { SpatRaster out = geometry(); if (x.empty()) { @@ -264,7 +264,6 @@ SpatRaster SpatRaster::distance_crds(std::vector& x, std::vector return(out); } - unsigned nc = ncol(); if (nrow() > 1000) { opt.steps = std::max(opt.steps, (size_t) 4); @@ -317,6 +316,9 @@ SpatRaster SpatRaster::distance_crds(std::vector& x, std::vector if (m != 1) { for (double &v : d) v *= m; } + if (max_dist > 0) { + for (double &v : d) v = v > max_dist ? NAN : v; + } if (!out.writeBlock(d, i)) return out; } readStop(); @@ -338,6 +340,9 @@ SpatRaster SpatRaster::distance_crds(std::vector& x, std::vector if (m != 1) { for (double &v : d) v *= m; } + if (max_dist > 0) { + for (double &v : d) v = v > max_dist ? NAN : v; + } if (!out.writeBlock(d, i)) return out; } } @@ -387,7 +392,7 @@ SpatRaster SpatRaster::distance_vector(SpatVector p, bool rasterize, std::string x = out.rasterize(p, "", {1}, NAN, false, "", false, false, false, ops); if (!lonlat) { - return x.distance(NAN, 0, false, unit, false, method, false, opt); + return x.distance(NAN, 0, false, unit, false, method, false, -1, opt); } if (poly) { @@ -408,7 +413,7 @@ SpatRaster SpatRaster::distance_vector(SpatVector p, bool rasterize, std::string bool setNA = false; - out = x.distance_crds(pxy[0], pxy[1], method, poly, setNA, unit, opt); + out = x.distance_crds(pxy[0], pxy[1], method, poly, setNA, unit, -1, opt); } else { @@ -459,7 +464,7 @@ SpatRaster SpatRaster::distance_vector(SpatVector p, bool rasterize, std::string std::vector> pxy = p.coordinates(); SpatOptions ops(opt); bool setNA = false; - out = distance_crds(pxy[0], pxy[1], method, false, setNA, unit, opt); + out = distance_crds(pxy[0], pxy[1], method, false, setNA, unit, -1, opt); } } return out; @@ -611,7 +616,7 @@ SpatRaster SpatRaster::distance_vector(SpatVector p, std::string unit, SpatOptio */ -SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, const std::string method, bool values, SpatOptions &opt) { +SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, const std::string method, bool values, double threshold, SpatOptions &opt) { SpatRaster out = geometry(1); if (!hasValues()) { @@ -631,7 +636,7 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std: std::vector lyr = {i}; SpatRaster r = subset(lyr, ops); ops.names = {nms[i]}; - r = r.distance(target, exclude, keepNA, unit, remove_zero, method, values, ops); + r = r.distance(target, exclude, keepNA, unit, remove_zero, method, values, threshold, ops); out.source[i] = r.source[0]; } if (!opt.get_filename().empty()) { @@ -640,7 +645,7 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std: return out; } if (!(values || is_lonlat())) { // && std::isnan(target) && std::isnan(exclude)) { - return proximity(target, exclude, keepNA, unit, false, 0, remove_zero, opt); + return proximity(target, exclude, keepNA, unit, false, threshold, remove_zero, opt); } bool setNA = false; @@ -656,9 +661,9 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std: } if (values) { std::vector> vv = extractXY(p[0], p[1], "", false); - return distance_crds_vals(p[0], p[1], vv[0], method, true, setNA, unit, opt); + return distance_crds_vals(p[0], p[1], vv[0], method, true, setNA, unit, threshold, opt); } else { - return distance_crds(p[0], p[1], method, true, setNA, unit, opt); + return distance_crds(p[0], p[1], method, true, setNA, unit, threshold, opt); } } else { x = replaceValues({exclude, target}, {NAN, NAN}, 1, false, NAN, false, ops); @@ -681,9 +686,9 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std: } if (values) { std::vector> vv = extractXY(p[0], p[1], "", false); - out = out.distance_crds_vals(p[0], p[1], vv[0], method, true, setNA, unit, opt); + out = out.distance_crds_vals(p[0], p[1], vv[0], method, true, setNA, unit, threshold, opt); } else { - out = out.distance_crds(p[0], p[1], method, true, setNA, unit, opt); + out = out.distance_crds(p[0], p[1], method, true, setNA, unit, threshold, opt); } return out; } diff --git a/src/distValueRaster.cpp b/src/distValueRaster.cpp index f0e60f1aa..f3c31ebff 100644 --- a/src/distValueRaster.cpp +++ b/src/distValueRaster.cpp @@ -223,7 +223,7 @@ void dist_only_vals(std::vector &d, std::vector &dv, const std:: } -SpatRaster SpatRaster::distance_crds_vals(std::vector& x, std::vector& y, const std::vector& v, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt) { +SpatRaster SpatRaster::distance_crds_vals(std::vector& x, std::vector& y, const std::vector& v, const std::string& method, bool skip, bool setNA, std::string unit, double maxdist, SpatOptions &opt) { SpatRaster out = geometry(); if (x.empty()) { @@ -236,6 +236,12 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector& x, std::vector 1000) { @@ -286,6 +292,17 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector& x, std::vector d, dv; dist_only_vals(d, dv, x, y, v, rxy[0], rxy[1], oldfirst, last, lonlat, dlast, dvlast, true, rv, method, setNA); oldfirst = first; + if (maxdist > 0) { + if (m != 1) { + for (size_t j=0; j maxdist) dv[j] = NAN; + } + } else { + for (size_t j=0; j maxdist) dv[j] = NAN; + } + } + } if (!out.writeBlock(dv, i)) return out; } readStop(); @@ -304,6 +321,17 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector& x, std::vector d, dv; dist_only_vals(d, dv, x, y, v, rxy[0], rxy[1], oldfirst, last, lonlat, dlast, dvlast, false, rv, method, setNA); oldfirst = first; + if (maxdist > 0) { + if (m != 1) { + for (size_t j=0; j maxdist) dv[j] = NAN; + } + } else { + for (size_t j=0; j maxdist) dv[j] = NAN; + } + } + } if (!out.writeBlock(dv, i)) return out; } } diff --git a/src/extract.cpp b/src/extract.cpp index b469f1f94..894e36de2 100644 --- a/src/extract.cpp +++ b/src/extract.cpp @@ -176,7 +176,7 @@ std::vector circ_dist(double xres, double yres, double d, size_t nrows, v[v.size()/2] = 1; SpatOptions opt; x.setValues(v, opt); - x = x.distance(NAN, NAN, false, "m", false, "cosine", false, opt); + x = x.distance(NAN, NAN, false, "m", false, "cosine", false, -1, opt); std::vector out; x.getValuesSource(0, out); diff --git a/src/gdal_algs.cpp b/src/gdal_algs.cpp index 07a9d315a..a7b29c0a6 100644 --- a/src/gdal_algs.cpp +++ b/src/gdal_algs.cpp @@ -1816,6 +1816,10 @@ SpatRaster SpatRaster::proximity(double target, double exclude, bool keepNA, std SpatRaster x; bool mask = false; std::vector mvals; + if ((!buffer) && (maxdist > 0)) { + papszOptions = CSLSetNameValue(papszOptions, "MAXDIST", doubleToAlmostChar(maxdist).c_str()); + } + if (buffer) { if (remove_zero) { x = isnotnan(true, ops); @@ -1833,6 +1837,7 @@ SpatRaster SpatRaster::proximity(double target, double exclude, bool keepNA, std } else { //option for keepNA does not work, perhaps because of int conversion //papszOptions = CSLSetNameValue(papszOptions, "USE_INPUT_NODATA", "YES"); + if (std::isnan(exclude)) { if (keepNA) { x = replaceValues({target, NAN}, {0, 0}, 1, true, 1, false, ops); diff --git a/src/spatRaster.h b/src/spatRaster.h index 3672855e7..31d1fb887 100644 --- a/src/spatRaster.h +++ b/src/spatRaster.h @@ -655,7 +655,7 @@ class SpatRaster { SpatRaster proximity(double target, double exclude, bool keepNA, std::string unit, bool buffer, double maxdist, bool remove_zero, SpatOptions &opt); SpatRaster fillNA(double missing, double maxdist, int niter, SpatOptions &opt); - SpatRaster distance(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, std::string method, bool values, SpatOptions &opt); + SpatRaster distance(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, std::string method, bool values, double threshold, SpatOptions &opt); SpatRaster nearest(double target, double exclude, bool keepNA, std::string unit, bool remove_zero, std::string method, SpatOptions &opt); // SpatRaster distance_spatvector(SpatVector p, std::string unit, const std::string& method, SpatOptions &opt); @@ -665,8 +665,8 @@ class SpatRaster { SpatRaster direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, const std::string& method, SpatOptions &opt); - SpatRaster distance_crds(std::vector& x, std::vector& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt); - SpatRaster distance_crds_vals(std::vector& x, std::vector& y, const std::vector& v, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt); + SpatRaster distance_crds(std::vector& x, std::vector& y, const std::string& method, bool skip, bool setNA, std::string unit,double threshold, SpatOptions &opt); + SpatRaster distance_crds_vals(std::vector& x, std::vector& y, const std::vector& v, const std::string& method, bool skip, bool setNA, std::string unit, double threshold, SpatOptions &opt); SpatRaster dn_crds(std::vector& x, std::vector& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt);