Skip to content

Commit

Permalink
m
Browse files Browse the repository at this point in the history
  • Loading branch information
@rhijmans
rhijmans committed Dec 27, 2024
1 parent b14c13e commit 4eaf505
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 59 deletions.
7 changes: 2 additions & 5 deletions man/distance.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ If \code{y} is \code{missing} this method computes the distance, for all cells t

If \code{y} is a numeric value, the cells with that value are ignored. That is, distance to or from these cells is not computed (only if \code{grid=FALSE}).

If \code{y} is a SpatVector, the distance to that SpatVector is computed for all cells. For lines and polygons this is done after rasterization; and only the overlapping areas of the vector and raster are considered (for now).
If \code{y} is a SpatVector, the distance to that SpatVector is computed for all cells, optionally after rasterization.

The distance is always expressed in meter if the coordinate reference system is longitude/latitude, and in map units otherwise. Map units are typically meter, but inspect \code{crs(x)} if in doubt.

Expand All @@ -33,9 +33,6 @@ If \code{y} is \code{missing}, a distance matrix between all objects in \code{x}

If \code{y} is a SpatVector the geographic distance between all objects is computed (and a matrix is returned). If both sets have the same number of points, and \code{pairwise=TRUE}, the distance between each pair of objects is computed, and a vector is returned.

The distance is always expressed in meter, except when the coordinate reference system is longitude/latitude AND one of the SpatVector(s) consists of lines or polygons. In that case the distance is in degrees, and thus not very useful (this will be fixed soon). Otherwise, results are more precise, sometimes much more precise, when using longitude/latitude rather than a planar coordinate reference system, as these distort distance.


If \code{x} is a \bold{matrix}:

\code{x} should consist of two columns, the first with "x" (or longitude) and the second with "y" coordinates (or latitude). If \code{y} is a also a matrix, the distance between each points in \code{x} and all points in \code{y} is computed, unless \code{pairwise=TRUE}
Expand Down Expand Up @@ -78,7 +75,7 @@ If \code{y} is missing, the distance between each points in \code{x} with all ot
}

\value{
SpatRaster or numeric or matrix or distance matrix (object of class "dist")
SpatRaster, numeric, matrix, or a distance matrix (object of class "dist")
}

\note{
Expand Down
145 changes: 91 additions & 54 deletions src/distVector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -237,74 +237,111 @@ std::vector<double> SpatVector::distance(bool sequential, std::string unit, cons
return(d);
}
std::string gtype = type();
if (gtype != "points") {
return geos_distance(sequential, "", m);
} else {
if (sequential) {
std::vector<std::vector<double>> p = coordinates();
size_t n = p[0].size();
d.reserve(n);
d.push_back(0);
n -= 1;
if (lonlat) {
std::function<double(double, double, double, double)> dfun;
if (method == "haversine") {
dfun = distHaversine;
} else if (method == "cosine") {
dfun = distCosine;
} else if (method == "geo") {
dfun = distLonlat;
} else {
setError("invalid lonlat distance method. Should be 'geo', 'cosine', or 'haversine'");
return(d);
}
for (size_t i=0; i<n; i++) {
d.push_back(
dfun(p[0][i], p[1][i], p[0][i+1], p[1][i+1]) * m
);
}
if (gtype == "points") {
if (lonlat) {
std::function<double(double, double, double, double)> dfun;
if (method == "haversine") {
dfun = distHaversine;
} else if (method == "cosine") {
dfun = distCosine;
} else if (method == "geo") {
dfun = distLonlat;
} else {
for (size_t i=0; i<n; i++) {
d.push_back(
distance_plane(p[0][i], p[1][i], p[0][i+1], p[1][i+1]) * m
);
}
setError("invalid lonlat distance method. Should be 'geo', 'cosine', or 'haversine'");
return(d);
}

} else {
size_t s = size();
size_t n = ((s-1) * s)/2;
d.reserve(n);
std::vector<std::vector<double>> p = coordinates();
if (lonlat) {
std::function<double(double, double, double, double)> dfun;
if (method == "haversine") {
dfun = distHaversine;
} else if (method == "cosine") {
dfun = distCosine;
} else if (method == "geo") {
dfun = distLonlat;
if (sequential) {
std::vector<std::vector<double>> p = coordinates();
size_t n = p[0].size();
d.reserve(n);
d.push_back(0);
n -= 1;
if (lonlat) {
for (size_t i=0; i<n; i++) {
d.push_back(
dfun(p[0][i], p[1][i], p[0][i+1], p[1][i+1]) * m
);
}
} else {
setError("invalid lonlat distance method. Should be 'geo', 'cosine', or 'haversine'");
return(d);
}

for (size_t i=0; i<(s-1); i++) {
for (size_t j=(i+1); j<s; j++) {
for (size_t i=0; i<n; i++) {
d.push_back(
dfun(p[0][i], p[1][i], p[0][j], p[1][j]) * m
distance_plane(p[0][i], p[1][i], p[0][i+1], p[1][i+1]) * m
);
}
}

} else {
size_t s = size();
size_t n = ((s-1) * s)/2;
d.reserve(n);
std::vector<std::vector<double>> p = coordinates();
if (lonlat) {
for (size_t i=0; i<(s-1); i++) {
for (size_t j=(i+1); j<s; j++) {
d.push_back(
dfun(p[0][i], p[1][i], p[0][j], p[1][j]) * m
);
}
}
} else {
for (size_t i=0; i<(s-1); i++) {
for (size_t j=(i+1); j<s; j++) {
d.push_back(
distance_plane(p[0][i], p[1][i], p[0][j], p[1][j]) * m
);
}
}
}
}
}
} else {
if (lonlat) {
std::function<double(double, double, double, double)> dfun;
if (method == "haversine") {
dfun = distHaversine;
} else if (method == "cosine") {
dfun = distCosine;
} else if (method == "geo") {
dfun = distLonlat;
} else {
setError("invalid lonlat distance method. Should be 'geo', 'cosine', or 'haversine'");
return(d);
}

size_t n = size();
d.reserve(n);
if (sequential) {
n -= 1;
SpatVector tmp1 = subset_rows({0});
std::vector<std::vector<double>> xy1 = tmp1.coordinates();
for (size_t i=0; i<n; i++) {
SpatVector tmp2 = subset_rows( {(int)i+1} );
std::vector<double> d1 = tmp2.nearestDistLonLat(xy1[0], xy1[1], unit, method);
std::vector<std::vector<double>> xy2 = tmp2.coordinates();
std::vector<double> d2 = tmp1.nearestDistLonLat(xy2[0], xy2[1], unit, method);
d.push_back(std::min(vmin(d1, false), vmin(d2, false)));
tmp1 = tmp2;
xy1 = xy2;
}
} else {
size_t s = size();
size_t n = ((s-1) * s)/2;
d.reserve(n);
for (size_t i=0; i<(s-1); i++) {
SpatVector tmp1 = subset_rows({(int)i});
std::vector<std::vector<double>> xy1 = tmp1.coordinates();
for (size_t j=(i+1); j<s; j++) {
d.push_back(
distance_plane(p[0][i], p[1][i], p[0][j], p[1][j]) * m
);
SpatVector tmp2 = subset_rows( {(int)j} );
std::vector<double> d1 = tmp2.nearestDistLonLat(xy1[0], xy1[1], unit, method);
std::vector<std::vector<double>> xy2 = tmp2.coordinates();
std::vector<double> d2 = tmp1.nearestDistLonLat(xy2[0], xy2[1], unit, method);
d.push_back(std::min(vmin(d1, false), vmin(d2, false)));
}
}
}
} else {
return geos_distance(sequential, "", m);
}
}

Expand Down

0 comments on commit 4eaf505

Please sign in to comment.