Skip to content

Commit

Permalink
maxdist
Browse files Browse the repository at this point in the history
  • Loading branch information
@rhijmans
rhijmans committed Jan 25, 2025
1 parent 1091750 commit d95d4be
Show file tree
Hide file tree
Showing 7 changed files with 64 additions and 23 deletions.
7 changes: 4 additions & 3 deletions R/distance.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'")
Expand All @@ -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")
}
)
Expand Down
8 changes: 5 additions & 3 deletions man/distance.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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="", ...)

Expand All @@ -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{
Expand All @@ -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}}}
Expand Down
29 changes: 17 additions & 12 deletions src/distRaster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,7 +244,7 @@ void dist_only(std::vector<double> &d, const std::vector<double>& vx, const std:
}


SpatRaster SpatRaster::distance_crds(std::vector<double>& x, std::vector<double>& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt) {
SpatRaster SpatRaster::distance_crds(std::vector<double>& x, std::vector<double>& y, const std::string& method, bool skip, bool setNA, std::string unit, double max_dist, SpatOptions &opt) {

SpatRaster out = geometry();
if (x.empty()) {
Expand All @@ -264,7 +264,6 @@ SpatRaster SpatRaster::distance_crds(std::vector<double>& x, std::vector<double>
return(out);
}


unsigned nc = ncol();
if (nrow() > 1000) {
opt.steps = std::max(opt.steps, (size_t) 4);
Expand Down Expand Up @@ -317,6 +316,9 @@ SpatRaster SpatRaster::distance_crds(std::vector<double>& x, std::vector<double>
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();
Expand All @@ -338,6 +340,9 @@ SpatRaster SpatRaster::distance_crds(std::vector<double>& x, std::vector<double>
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;
}
}
Expand Down Expand Up @@ -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) {
Expand All @@ -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 {

Expand Down Expand Up @@ -459,7 +464,7 @@ SpatRaster SpatRaster::distance_vector(SpatVector p, bool rasterize, std::string
std::vector<std::vector<double>> 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;
Expand Down Expand Up @@ -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()) {
Expand All @@ -631,7 +636,7 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std:
std::vector<size_t> 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()) {
Expand All @@ -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;
Expand All @@ -656,9 +661,9 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std:
}
if (values) {
std::vector<std::vector<double>> 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);
Expand All @@ -681,9 +686,9 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std:
}
if (values) {
std::vector<std::vector<double>> 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;
}
Expand Down
30 changes: 29 additions & 1 deletion src/distValueRaster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ void dist_only_vals(std::vector<double> &d, std::vector<double> &dv, const std::
}


SpatRaster SpatRaster::distance_crds_vals(std::vector<double>& x, std::vector<double>& y, const std::vector<double>& v, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt) {
SpatRaster SpatRaster::distance_crds_vals(std::vector<double>& x, std::vector<double>& y, const std::vector<double>& v, const std::string& method, bool skip, bool setNA, std::string unit, double maxdist, SpatOptions &opt) {

SpatRaster out = geometry();
if (x.empty()) {
Expand All @@ -236,6 +236,12 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector<double>& x, std::vector<do
permute(y, pm);

bool lonlat = is_lonlat();
double m=1;
if (!source[0].srs.m_dist(m, lonlat, unit)) {
out.setError("invalid unit");
return(out);
}


unsigned nc = ncol();
if (nrow() > 1000) {
Expand Down Expand Up @@ -286,6 +292,17 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector<double>& x, std::vector<do
std::vector<double> 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<dv.size(); j++) {
if ((d[j] / m) > maxdist) dv[j] = NAN;
}
} else {
for (size_t j=0; j<dv.size(); j++) {
if (d[j] > maxdist) dv[j] = NAN;
}
}
}
if (!out.writeBlock(dv, i)) return out;
}
readStop();
Expand All @@ -304,6 +321,17 @@ SpatRaster SpatRaster::distance_crds_vals(std::vector<double>& x, std::vector<do
std::vector<double> 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<dv.size(); j++) {
if ((d[j] / m) > maxdist) dv[j] = NAN;
}
} else {
for (size_t j=0; j<dv.size(); j++) {
if (d[j] > maxdist) dv[j] = NAN;
}
}
}
if (!out.writeBlock(dv, i)) return out;
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/extract.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ std::vector<double> 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<double> out;
x.getValuesSource(0, out);
Expand Down
5 changes: 5 additions & 0 deletions src/gdal_algs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1816,6 +1816,10 @@ SpatRaster SpatRaster::proximity(double target, double exclude, bool keepNA, std
SpatRaster x;
bool mask = false;
std::vector<double> mvals;
if ((!buffer) && (maxdist > 0)) {
papszOptions = CSLSetNameValue(papszOptions, "MAXDIST", doubleToAlmostChar(maxdist).c_str());
}

if (buffer) {
if (remove_zero) {
x = isnotnan(true, ops);
Expand All @@ -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);
Expand Down
6 changes: 3 additions & 3 deletions src/spatRaster.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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<double>& x, std::vector<double>& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt);
SpatRaster distance_crds_vals(std::vector<double>& x, std::vector<double>& y, const std::vector<double>& v, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt);
SpatRaster distance_crds(std::vector<double>& x, std::vector<double>& y, const std::string& method, bool skip, bool setNA, std::string unit,double threshold, SpatOptions &opt);
SpatRaster distance_crds_vals(std::vector<double>& x, std::vector<double>& y, const std::vector<double>& v, const std::string& method, bool skip, bool setNA, std::string unit, double threshold, SpatOptions &opt);

SpatRaster dn_crds(std::vector<double>& x, std::vector<double>& y, const std::string& method, bool skip, bool setNA, std::string unit, SpatOptions &opt);

Expand Down

0 comments on commit d95d4be

Please sign in to comment.