From e0aeea9e023dfccc4d6f9667192117dce15351d4 Mon Sep 17 00:00:00 2001 From: rhijmans Date: Thu, 26 Dec 2024 10:44:25 -0800 Subject: [PATCH] m --- R/RcppExports.R | 52 +---------- R/distance.R | 14 +-- R/ncdf.R | 2 +- R/panel.R | 6 +- man/direction.Rd | 7 +- man/distance.Rd | 8 +- man/panel.Rd | 3 +- src/RcppExports.cpp | 213 ++------------------------------------------ src/RcppModule.cpp | 2 +- src/distRaster.cpp | 112 +++++++---------------- src/distance.cpp | 109 ++++++++++++++++------- src/distance.h | 2 +- src/geosphere.cpp | 141 ++++++++++++----------------- src/geosphere.h | 25 ++++++ src/spatRaster.h | 6 +- src/spatVector.h | 2 +- 16 files changed, 225 insertions(+), 479 deletions(-) create mode 100644 src/geosphere.h diff --git a/R/RcppExports.R b/R/RcppExports.R index 4192394de..1cc594430 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -117,56 +117,8 @@ rgb2hex <- function(x) { .Call(`_terra_uniqueSymmetricRows`, x, y) } -dist_lonlat <- function(lon1, lat1, lon2, lat2) { - .Call(`_terra_dist_lonlat`, lon1, lat1, lon2, lat2) -} - -dist_cosine <- function(lon1, lat1, lon2, lat2, r = 6378137) { - .Call(`_terra_dist_cosine`, lon1, lat1, lon2, lat2, r) -} - -dist_cosine_rad <- function(lon1, lat1, lon2, lat2, r = 6378137) { - .Call(`_terra_dist_cosine_rad`, lon1, lat1, lon2, lat2, r) -} - -dest_lonlat <- function(slon, slat, sazi, dist, dlon, dlat, dazi) { - invisible(.Call(`_terra_dest_lonlat`, slon, slat, sazi, dist, dlon, dlat, dazi)) -} - -dir_lonlat <- function(lon1, lat1, lon2, lat2) { - .Call(`_terra_dir_lonlat`, lon1, lat1, lon2, lat2) -} - -dir_rad <- function(lon1, lat1, lon2, lat2) { - .Call(`_terra_dir_rad`, lon1, lat1, lon2, lat2) -} - -dist2track <- function(lon1, lat1, lon2, lat2, plon, plat, sign, r = 6378137) { - .Call(`_terra_dist2track`, lon1, lat1, lon2, lat2, plon, plat, sign, r) -} - -dist2track_cosine_rad <- function(lon1, lat1, lon2, lat2, plon, plat, sign, r = 6378137) { - .Call(`_terra_dist2track_cosine_rad`, lon1, lat1, lon2, lat2, plon, plat, sign, r) -} - -alongTrackDistance <- function(lon1, lat1, lon2, lat2, plon, plat, r = 6378137) { - .Call(`_terra_alongTrackDistance`, lon1, lat1, lon2, lat2, plon, plat, r) -} - -alongTrackDistance_rad <- function(lon1, lat1, lon2, lat2, plon, plat, r = 6378137) { - .Call(`_terra_alongTrackDistance_rad`, lon1, lat1, lon2, lat2, plon, plat, r) -} - -dist2segment <- function(plon, plat, lon1, lat1, lon2, lat2) { - .Call(`_terra_dist2segment`, plon, plat, lon1, lat1, lon2, lat2) -} - -dist2segment_cosine_rad <- function(plon, plat, lon1, lat1, lon2, lat2, r = 6378137) { - .Call(`_terra_dist2segment_cosine_rad`, plon, plat, lon1, lat1, lon2, lat2, r) -} - -dist2segmentPoint <- function(plon, plat, lon1, lat1, lon2, lat2, ilon, ilat) { - .Call(`_terra_dist2segmentPoint`, plon, plat, lon1, lat1, lon2, lat2, ilon, ilat) +dist2segmentPoint_geo <- function(plon, plat, lon1, lat1, lon2, lat2, ilon, ilat) { + .Call(`_terra_dist2segmentPoint_geo`, plon, plat, lon1, lat1, lon2, lat2, ilon, ilat) } intermediate <- function(lon1, lat1, lon2, lat2, n, distance) { diff --git a/R/distance.R b/R/distance.R index 88418bf53..3566f09dd 100644 --- a/R/distance.R +++ b/R/distance.R @@ -66,7 +66,7 @@ setMethod("gridDist", signature(x="SpatRaster"), setMethod("distance", signature(x="SpatRaster", y="SpatVector"), - function(x, y, unit="m", rasterize=FALSE, method="haversine", filename="", ...) { + function(x, y, unit="m", rasterize=FALSE, method="cosine", filename="", ...) { opt <- spatOptions(filename, ...) unit <- as.character(unit[1]) @@ -85,7 +85,7 @@ setMethod("distance", signature(x="SpatRaster", y="SpatVector"), setMethod("distance", signature(x="SpatRaster", y="sf"), - function(x, y, unit="m", rasterize=FALSE, method="haversine", filename="", ...) { + function(x, y, unit="m", rasterize=FALSE, method="cosine", filename="", ...) { distance(x, vect(y), unit=unit, rasterize=rasterize, method=method, filename=filename, ...) } ) @@ -113,7 +113,7 @@ mat2wide <- function(m, sym=TRUE, keep=NULL) { } setMethod("distance", signature(x="SpatVector", y="ANY"), - function(x, y, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m", method="geo") { + function(x, y, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m", method="cosine") { if (!missing(y)) { error("distance", "If 'x' is a SpatVector, 'y' should be a SpatVector or missing") } @@ -144,7 +144,7 @@ setMethod("distance", signature(x="SpatVector", y="ANY"), setMethod("distance", signature(x="SpatVector", y="SpatVector"), - function(x, y, pairwise=FALSE, unit="m", method = "geo") { + function(x, y, pairwise=FALSE, unit="m", method = "cosine") { unit <- as.character(unit[1]) d <- x@pntr$distance_other(y@pntr, pairwise, unit, method) messages(x, "distance") @@ -205,15 +205,15 @@ setMethod("distance", signature(x="matrix", y="missing"), setMethod("distance", signature(x="data.frame", y="missing"), function(x, y, lonlat=NULL, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m", method="geo") { - distance(as.matrix(x), lonlat=lonlat, sequential=sequential, pairs=pairs, symmetrical=symmetrical, unit=unit, method="geo") + distance(as.matrix(x), lonlat=lonlat, sequential=sequential, pairs=pairs, symmetrical=symmetrical, unit=unit, method=method) } ) setMethod("direction", signature(x="SpatRaster"), - function(x, from=FALSE, degrees=FALSE, filename="", ...) { + function(x, from=FALSE, degrees=FALSE, method="cosine", filename="", ...) { opt <- spatOptions(filename, ...) - x@pntr <- x@pntr$rastDirection(from[1], degrees[1], NA, NA, opt) + x@pntr <- x@pntr$rastDirection(from[1], degrees[1], NA, NA, method, opt) messages(x, "direction") } ) diff --git a/R/ncdf.R b/R/ncdf.R index 00d3c66f6..34131a8d2 100644 --- a/R/ncdf.R +++ b/R/ncdf.R @@ -220,7 +220,7 @@ write_tags <- function(tags, nc, varid, prefix="TAG_") { } write_tags(metags(y), ncobj, ncvars[[i]], "") } - if (progress) utils::close(pb) + if (progress) close(pb) ncdf4::ncatt_put(ncobj, 0, "Conventions", "CF-1.4", prec="text") pkgversion <- drop(read.dcf(file=system.file("DESCRIPTION", package="terra"), fields=c("Version"))) diff --git a/R/panel.R b/R/panel.R index eee3484b1..8d3b1f140 100644 --- a/R/panel.R +++ b/R/panel.R @@ -1,7 +1,7 @@ setMethod("panel", signature(x="SpatRaster"), function(x, main, loc.main="topleft", nc, nr, maxnl=16, maxcell=500000, - box=FALSE, pax=list(), plg=list(), range=NULL, ...) { + box=FALSE, pax=list(), plg=list(), range=NULL, halo=TRUE, ...) { dots <- list(...) if (!is.null(dots$type)) { @@ -97,10 +97,10 @@ setMethod("panel", signature(x="SpatRaster"), y <- x[[i]] levels(y) <- lv plot(y, 1, main=main[i], mar=mar, legend=legend[i], pax=pax, box=box, - loc.main=loc.main, halo=TRUE, plg=plg, type="classes", ...) + loc.main=loc.main, halo=halo, plg=plg, type="classes", ...) } else { plot(x, i, main=main[i], mar=mar, legend=legend[i], range=range, pax=pax, box=box, - loc.main=loc.main, halo=TRUE, plg=plg, type=ptype, ...) + loc.main=loc.main, halo=halo, plg=plg, type=ptype, ...) } } } diff --git a/man/direction.Rd b/man/direction.Rd index 75056413a..2a15a4eea 100644 --- a/man/direction.Rd +++ b/man/direction.Rd @@ -10,14 +10,15 @@ The direction (azimuth) to or from the nearest cell that is not \code{NA}. The d } \usage{ -\S4method{direction}{SpatRaster}(x, from=FALSE, degrees=FALSE, filename="", ...) +\S4method{direction}{SpatRaster}(x, from=FALSE, degrees=FALSE, method="cosine", filename="", ...) } \arguments{ \item{x}{SpatRaster} -\item{filename}{Character. Output filename (optional)} -\item{degrees}{Logical. If \code{FALSE} (the default) the unit of direction is radians.} \item{from}{Logical. Default is \code{FALSE}. If \code{TRUE}, the direction from (instead of to) the nearest cell that is not \code{NA} is returned} +\item{degrees}{Logical. If \code{FALSE} (the default) the unit of direction is radians.} +\item{method}{character. Should be "geo", or "cosine". With "geo" the most precise but slower geodesic method of Karney (2003) is used. The "cosine" method is faster but less precise} +\item{filename}{Character. Output filename (optional)} \item{...}{Additional arguments as for \code{\link{writeRaster}}} } diff --git a/man/distance.Rd b/man/distance.Rd index 3239b0d40..9badd4a5e 100644 --- a/man/distance.Rd +++ b/man/distance.Rd @@ -48,11 +48,11 @@ 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", filename="", ...) -\S4method{distance}{SpatRaster,SpatVector}(x, y, unit="m", rasterize=FALSE, method="haversine", filename="", ...) +\S4method{distance}{SpatRaster,SpatVector}(x, y, unit="m", rasterize=FALSE, method="cosine", filename="", ...) -\S4method{distance}{SpatVector,ANY}(x, y, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m", method="geo") +\S4method{distance}{SpatVector,ANY}(x, y, sequential=FALSE, pairs=FALSE, symmetrical=TRUE, unit="m", method="cosine") -\S4method{distance}{SpatVector,SpatVector}(x, y, pairwise=FALSE, unit="m", method="geo") +\S4method{distance}{SpatVector,SpatVector}(x, y, pairwise=FALSE, unit="m", method="cosine") \S4method{distance}{matrix,matrix}(x, y, lonlat, pairwise=FALSE, unit="m", method="geo") @@ -66,7 +66,7 @@ If \code{y} is missing, the distance between each points in \code{x} with all ot \item{target}{numeric. The value of the cells for which distances to cells that are not \code{NA} should be computed} \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", "haversine", "cosine". With "geo" the most precise but slower method of Karney (2003) is used. The other two methods are faster but less precise} + \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{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/man/panel.Rd b/man/panel.Rd index 334d85ee5..3762f9bd8 100644 --- a/man/panel.Rd +++ b/man/panel.Rd @@ -13,7 +13,7 @@ Show multiple maps that share a single legend. \usage{ \S4method{panel}{SpatRaster}(x, main, loc.main="topleft", nc, nr, maxnl=16, - maxcell=500000, box=FALSE, pax=list(), plg=list(), range=NULL, ...) + maxcell=500000, box=FALSE, pax=list(), plg=list(), range=NULL, halo=TRUE, ...) } \arguments{ @@ -28,6 +28,7 @@ Show multiple maps that share a single legend. \item{plg}{see \code{\link{plot}}} \item{pax}{see \code{\link{plot}}} \item{range}{numeric. minimum and maximum values to be used for the continuous legend } + \item{halo}{logical. Use a halo around main (the title)?} \item{...}{arguments passed to \code{plot("SpatRaster", "numeric")} and additional graphical arguments} } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1455ae19f..585fc7977 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -336,200 +336,9 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// dist_lonlat -double dist_lonlat(const double& lon1, const double& lat1, const double& lon2, const double& lat2); -RcppExport SEXP _terra_dist_lonlat(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const double& >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< const double& >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< const double& >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< const double& >::type lat2(lat2SEXP); - rcpp_result_gen = Rcpp::wrap(dist_lonlat(lon1, lat1, lon2, lat2)); - return rcpp_result_gen; -END_RCPP -} -// dist_cosine -double dist_cosine(double lon1, double lat1, double lon2, double lat2, const double& r); -RcppExport SEXP _terra_dist_cosine(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< const double& >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(dist_cosine(lon1, lat1, lon2, lat2, r)); - return rcpp_result_gen; -END_RCPP -} -// dist_cosine_rad -double dist_cosine_rad(const double& lon1, const double& lat1, const double& lon2, const double& lat2, const double& r); -RcppExport SEXP _terra_dist_cosine_rad(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const double& >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< const double& >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< const double& >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< const double& >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< const double& >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(dist_cosine_rad(lon1, lat1, lon2, lat2, r)); - return rcpp_result_gen; -END_RCPP -} -// dest_lonlat -void dest_lonlat(double slon, double slat, double sazi, double dist, double& dlon, double& dlat, double& dazi); -RcppExport SEXP _terra_dest_lonlat(SEXP slonSEXP, SEXP slatSEXP, SEXP saziSEXP, SEXP distSEXP, SEXP dlonSEXP, SEXP dlatSEXP, SEXP daziSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type slon(slonSEXP); - Rcpp::traits::input_parameter< double >::type slat(slatSEXP); - Rcpp::traits::input_parameter< double >::type sazi(saziSEXP); - Rcpp::traits::input_parameter< double >::type dist(distSEXP); - Rcpp::traits::input_parameter< double& >::type dlon(dlonSEXP); - Rcpp::traits::input_parameter< double& >::type dlat(dlatSEXP); - Rcpp::traits::input_parameter< double& >::type dazi(daziSEXP); - dest_lonlat(slon, slat, sazi, dist, dlon, dlat, dazi); - return R_NilValue; -END_RCPP -} -// dir_lonlat -double dir_lonlat(double lon1, double lat1, double lon2, double lat2); -RcppExport SEXP _terra_dir_lonlat(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - rcpp_result_gen = Rcpp::wrap(dir_lonlat(lon1, lat1, lon2, lat2)); - return rcpp_result_gen; -END_RCPP -} -// dir_rad -double dir_rad(double lon1, double lat1, double lon2, double lat2); -RcppExport SEXP _terra_dir_rad(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - rcpp_result_gen = Rcpp::wrap(dir_rad(lon1, lat1, lon2, lat2)); - return rcpp_result_gen; -END_RCPP -} -// dist2track -double dist2track(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r); -RcppExport SEXP _terra_dist2track(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP plonSEXP, SEXP platSEXP, SEXP signSEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< bool >::type sign(signSEXP); - Rcpp::traits::input_parameter< double >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(dist2track(lon1, lat1, lon2, lat2, plon, plat, sign, r)); - return rcpp_result_gen; -END_RCPP -} -// dist2track_cosine_rad -double dist2track_cosine_rad(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r); -RcppExport SEXP _terra_dist2track_cosine_rad(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP plonSEXP, SEXP platSEXP, SEXP signSEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< bool >::type sign(signSEXP); - Rcpp::traits::input_parameter< double >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(dist2track_cosine_rad(lon1, lat1, lon2, lat2, plon, plat, sign, r)); - return rcpp_result_gen; -END_RCPP -} -// alongTrackDistance -double alongTrackDistance(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r); -RcppExport SEXP _terra_alongTrackDistance(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP plonSEXP, SEXP platSEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< double >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(alongTrackDistance(lon1, lat1, lon2, lat2, plon, plat, r)); - return rcpp_result_gen; -END_RCPP -} -// alongTrackDistance_rad -double alongTrackDistance_rad(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r); -RcppExport SEXP _terra_alongTrackDistance_rad(SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP plonSEXP, SEXP platSEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< double >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(alongTrackDistance_rad(lon1, lat1, lon2, lat2, plon, plat, r)); - return rcpp_result_gen; -END_RCPP -} -// dist2segment -double dist2segment(double plon, double plat, double lon1, double lat1, double lon2, double lat2); -RcppExport SEXP _terra_dist2segment(SEXP plonSEXP, SEXP platSEXP, SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - rcpp_result_gen = Rcpp::wrap(dist2segment(plon, plat, lon1, lat1, lon2, lat2)); - return rcpp_result_gen; -END_RCPP -} -// dist2segment_cosine_rad -double dist2segment_cosine_rad(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double r); -RcppExport SEXP _terra_dist2segment_cosine_rad(SEXP plonSEXP, SEXP platSEXP, SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP rSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type plon(plonSEXP); - Rcpp::traits::input_parameter< double >::type plat(platSEXP); - Rcpp::traits::input_parameter< double >::type lon1(lon1SEXP); - Rcpp::traits::input_parameter< double >::type lat1(lat1SEXP); - Rcpp::traits::input_parameter< double >::type lon2(lon2SEXP); - Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); - Rcpp::traits::input_parameter< double >::type r(rSEXP); - rcpp_result_gen = Rcpp::wrap(dist2segment_cosine_rad(plon, plat, lon1, lat1, lon2, lat2, r)); - return rcpp_result_gen; -END_RCPP -} -// dist2segmentPoint -double dist2segmentPoint(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double& ilon, double& ilat); -RcppExport SEXP _terra_dist2segmentPoint(SEXP plonSEXP, SEXP platSEXP, SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP ilonSEXP, SEXP ilatSEXP) { +// dist2segmentPoint_geo +double dist2segmentPoint_geo(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double& ilon, double& ilat); +RcppExport SEXP _terra_dist2segmentPoint_geo(SEXP plonSEXP, SEXP platSEXP, SEXP lon1SEXP, SEXP lat1SEXP, SEXP lon2SEXP, SEXP lat2SEXP, SEXP ilonSEXP, SEXP ilatSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -541,7 +350,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type lat2(lat2SEXP); Rcpp::traits::input_parameter< double& >::type ilon(ilonSEXP); Rcpp::traits::input_parameter< double& >::type ilat(ilatSEXP); - rcpp_result_gen = Rcpp::wrap(dist2segmentPoint(plon, plat, lon1, lat1, lon2, lat2, ilon, ilat)); + rcpp_result_gen = Rcpp::wrap(dist2segmentPoint_geo(plon, plat, lon1, lat1, lon2, lat2, ilon, ilat)); return rcpp_result_gen; END_RCPP } @@ -594,19 +403,7 @@ static const R_CallMethodDef CallEntries[] = { {"_terra_pearson_cor", (DL_FUNC) &_terra_pearson_cor, 3}, {"_terra_weighted_pearson_cor", (DL_FUNC) &_terra_weighted_pearson_cor, 4}, {"_terra_uniqueSymmetricRows", (DL_FUNC) &_terra_uniqueSymmetricRows, 2}, - {"_terra_dist_lonlat", (DL_FUNC) &_terra_dist_lonlat, 4}, - {"_terra_dist_cosine", (DL_FUNC) &_terra_dist_cosine, 5}, - {"_terra_dist_cosine_rad", (DL_FUNC) &_terra_dist_cosine_rad, 5}, - {"_terra_dest_lonlat", (DL_FUNC) &_terra_dest_lonlat, 7}, - {"_terra_dir_lonlat", (DL_FUNC) &_terra_dir_lonlat, 4}, - {"_terra_dir_rad", (DL_FUNC) &_terra_dir_rad, 4}, - {"_terra_dist2track", (DL_FUNC) &_terra_dist2track, 8}, - {"_terra_dist2track_cosine_rad", (DL_FUNC) &_terra_dist2track_cosine_rad, 8}, - {"_terra_alongTrackDistance", (DL_FUNC) &_terra_alongTrackDistance, 7}, - {"_terra_alongTrackDistance_rad", (DL_FUNC) &_terra_alongTrackDistance_rad, 7}, - {"_terra_dist2segment", (DL_FUNC) &_terra_dist2segment, 6}, - {"_terra_dist2segment_cosine_rad", (DL_FUNC) &_terra_dist2segment_cosine_rad, 7}, - {"_terra_dist2segmentPoint", (DL_FUNC) &_terra_dist2segmentPoint, 8}, + {"_terra_dist2segmentPoint_geo", (DL_FUNC) &_terra_dist2segmentPoint_geo, 8}, {"_terra_intermediate", (DL_FUNC) &_terra_intermediate, 6}, {"_rcpp_module_boot_spat", (DL_FUNC) &_rcpp_module_boot_spat, 0}, {NULL, NULL, 0} diff --git a/src/RcppModule.cpp b/src/RcppModule.cpp index ba8716da9..758ca9980 100644 --- a/src/RcppModule.cpp +++ b/src/RcppModule.cpp @@ -460,7 +460,7 @@ RCPP_MODULE(spat){ .method("linesList", &SpatVector::linesList) .method("polygonsList", &SpatVector::polygonsList) .method("linesNA", &SpatVector::linesNA) - .method("linedistlonlat", &SpatVector::linedistLonLat) + .method("linedistlonlat", &SpatVector::point2lineDistLonLat) .method("add_column_empty", (void (SpatVector::*)(unsigned dtype, std::string name))( &SpatVector::add_column)) .method("add_column_double", (bool (SpatVector::*)(std::vector, std::string name))( &SpatVector::add_column)) diff --git a/src/distRaster.cpp b/src/distRaster.cpp index d5603123f..54cd9fdf7 100644 --- a/src/distRaster.cpp +++ b/src/distRaster.cpp @@ -37,9 +37,9 @@ inline void shortDistPoints(std::vector &d, const std::vector &x } } -inline void shortDirectPoints(std::vector &d, const std::vector &x, const std::vector &y, const std::vector &px, const std::vector &py, const bool& lonlat, bool &from, bool °rees) { +inline void shortDirectPoints(std::vector &d, std::vector &x, std::vector &y, std::vector &px, std::vector &py, const bool& lonlat, bool &from, bool °rees, const std::string &method) { if (lonlat) { - directionToNearest_lonlat(d, x, y, px, py, degrees, from); + directionToNearest_lonlat(d, x, y, px, py, degrees, from, method); } else { directionToNearest_plane(d, x, y, px, py, degrees, from); } @@ -78,6 +78,7 @@ std::vector dist_bounds(const std::vector& vx, const std::vector first = vx.size(); last = 0; + if (lonlat) { std::function dfun; if (method == "haversine") { @@ -462,7 +463,7 @@ SpatRaster SpatRaster::distance_rasterize(SpatVector p, double target, double ex -SpatRaster SpatRaster::direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, SpatOptions &opt) { +SpatRaster SpatRaster::direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, const std::string &method, SpatOptions &opt) { SpatRaster out = geometry(); if (source[0].srs.wkt.empty()) { @@ -498,9 +499,6 @@ SpatRaster SpatRaster::direction_rasterize(SpatVector p, bool from, bool degrees out.setError("no locations to compute from"); return(out); } - - - unsigned nc = ncol(); if (!readStart()) { out.setError(getError()); @@ -520,69 +518,25 @@ SpatRaster SpatRaster::direction_rasterize(SpatVector p, bool from, bool degrees std::iota(cells.begin(), cells.end(), out.bs.row[i] * nc); -/* - if (gtype == "points") { - readBlock(v, out.bs, i); - if (std::isnan(target)) { - if (std::isnan(exclude)) { - for (size_t j=0; j> xy = xyFromCell(cells); - shortDirectPoints(vals, xy[0], xy[1], pxy[0], pxy[1], lonlat, from, degrees); + shortDirectPoints(vals, xy[0], xy[1], pxy[0], pxy[1], lonlat, from, degrees, method); if (!out.writeBlock(vals, i)) return out; } @@ -721,7 +675,7 @@ SpatRaster SpatRaster::distance(double target, double exclude, bool keepNA, std: -SpatRaster SpatRaster::direction(bool from, bool degrees, double target, double exclude, SpatOptions &opt) { +SpatRaster SpatRaster::direction(bool from, bool degrees, double target, double exclude, const std::string &method, SpatOptions &opt) { SpatRaster out = geometry(1); if (!hasValues()) { out.setError("SpatRaster has no values"); @@ -740,7 +694,7 @@ SpatRaster SpatRaster::direction(bool from, bool degrees, double target, double std::vector lyr = {i}; SpatRaster r = subset(lyr, ops); ops.names = {nms[i]}; - r = r.direction(from, degrees, target, exclude, ops); + r = r.direction(from, degrees, target, exclude, method, ops); out.source[i] = r.source[0]; } if (!opt.get_filename().empty()) { @@ -761,7 +715,7 @@ SpatRaster SpatRaster::direction(bool from, bool degrees, double target, double out.setError("no cells to compute direction from or to"); return(out); } - return direction_rasterize(p, from, degrees, target, exclude, opt); + return direction_rasterize(p, from, degrees, target, exclude, method, opt); } @@ -995,7 +949,7 @@ std::vector SpatVector::distance(SpatVector x, bool pairwise, std::strin } if (pairwise && (s != sx ) && (s > 1) && (sx > 1)) { - setError("Can only do pairwise distance if geometries match, or if one is a single geometry"); + setError("For pairwise distance, the number of geometries must match, or one should have a single geometry"); return(d); } @@ -1012,30 +966,24 @@ std::vector SpatVector::distance(SpatVector x, bool pairwise, std::strin if ((gtype != "points") || (xtype != "points")) { if (lonlat) { - if (xtype == "points") { - return linedistLonLat(x, unit); - } else if (gtype == "points") { - return x.linedistLonLat(*this, unit); + if (gtype == "points") { + return x.point2lineDistLonLat(*this, unit, method); } else { - // not good enough - // need fixing - SpatVector tmp = x.as_points(false, true); - return x.linedistLonLat(tmp, unit); + return point2lineDistLonLat(x, unit, method); } + } else { + std::string distfun=""; + d = geos_distance(x, pairwise, distfun); + if (m != 1) { + for (double &i : d) i *= m; + } + return d; } - - std::string distfun=""; - d = geos_distance(x, pairwise, distfun); - if (m != 1) { - for (double &i : d) i *= m; - } - return d; } std::vector> p = coordinates(); std::vector> px = x.coordinates(); - return pointdistance(p[0], p[1], px[0], px[1], pairwise, m, lonlat, method); } diff --git a/src/distance.cpp b/src/distance.cpp index 838de4e0f..0bee63e53 100644 --- a/src/distance.cpp +++ b/src/distance.cpp @@ -15,7 +15,6 @@ // You should have received a copy of the GNU General Public License // along with spat. If not, see . - #ifndef M_PI #define M_PI (3.14159265358979323846) #endif @@ -24,6 +23,7 @@ //#include #include #include "geodesic.h" +#include "geosphere.h" #include "recycle.h" #include #include @@ -33,7 +33,6 @@ double toRad(double °) { return( deg * 0.0174532925199433 ); } - double toDeg(double &rad) { return( rad * 57.2957795130823 ); } @@ -218,45 +217,91 @@ std::vector direction_lonlat(std::vector lon1, std::vector &azi, const std::vector &lon1, const std::vector &lat1, const std::vector &lon2, const std::vector &lat2, bool& degrees, bool& from) { +inline void deg2rad(std::vector &x) { + const double f = 0.0174532925199433; + for (double& d : x) d *= f; +} - double a = 6378137.0; - double f = 1/298.257223563; - double azi1, azi2, s12, dist; - size_t n = lon1.size(); - size_t m = lon2.size(); - azi.resize(n, NAN); +void directionToNearest_lonlat(std::vector &azi, std::vector &lon1, std::vector &lat1, std::vector &lon2, std::vector &lat2, bool& degrees, bool& from, const std::string &method) { - struct geod_geodesic g; - geod_init(&g, a, f); + if (method == "geo") { + double a = 6378137.0; + double f = 1/298.257223563; - for (size_t i=0; i < n; i++) { - if (std::isnan(lat1[i])) { - azi[i] = NAN; - continue; - } - geod_inverse(&g, lat1[i], lon1[i], lat2[0], lon2[0], &dist, &azi1, &azi2); - size_t minj=0; - azi[i] = azi1; - for (size_t j=1; j direction_lonlat(std::vector lon1, std::vector lat1, std::vector lon2, std::vector lat2, bool degrees, const std::string& method); -void directionToNearest_lonlat(std::vector &azi, const std::vector &lon1, const std::vector &lat1, const std::vector &lon2, const std::vector &lat2, bool °rees, bool &from); +void directionToNearest_lonlat(std::vector &azi, std::vector &lon1, std::vector &lat1, std::vector &lon2, std::vector &lat2, bool °rees, bool &from, const std::string &method); void distanceCosineToNearest_lonlat(std::vector &d, const std::vector &lon1, const std::vector &lat1, const std::vector &lon2, const std::vector &lat2); diff --git a/src/geosphere.cpp b/src/geosphere.cpp index 579187c4a..2eccfca5c 100644 --- a/src/geosphere.cpp +++ b/src/geosphere.cpp @@ -18,6 +18,7 @@ #include "spatVector.h" #include "geodesic.h" #include "recycle.h" +#include "geosphere.h" #ifndef M_PI #define M_PI (3.1415926535897932384626433) @@ -67,10 +68,7 @@ inline double get_sign(const double &x) { } -// [[Rcpp::export]] -double dist_lonlat(const double &lon1, const double &lat1, const double &lon2, const double &lat2) { - //double a = 6378137.0; - //double f = 1/298.257223563; +double distance_geo(const double &lon1, const double &lat1, const double &lon2, const double &lat2) { double s12, azi1, azi2; struct geod_geodesic g; geod_init(&g, WGS84_a, WGS84_f); @@ -79,8 +77,7 @@ double dist_lonlat(const double &lon1, const double &lat1, const double &lon2, c } -// [[Rcpp::export]] -double dist_cosine(double lon1, double lat1, double lon2, double lat2, const double &r = 6378137) { +double distance_cosdeg(double lon1, double lat1, double lon2, double lat2, const double &r = 6378137) { deg2rad(lon1); deg2rad(lon2); deg2rad(lat1); @@ -88,27 +85,16 @@ double dist_cosine(double lon1, double lat1, double lon2, double lat2, const dou return r * acos((sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1-lon2))); } -// [[Rcpp::export]] -double dist_cosine_rad(const double &lon1, const double &lat1, const double &lon2, const double &lat2, const double &r = 6378137) { - return r * acos((sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1-lon2))); -} -// [[Rcpp::export]] -void dest_lonlat(double slon, double slat, double sazi, double dist, double &dlon, double &dlat, double &dazi) { - //double a = 6378137.0; - //double f = 1/298.257223563; +void dest_geo(double slon, double slat, double sazi, double dist, double &dlon, double &dlat, double &dazi) { struct geod_geodesic g; geod_init(&g, WGS84_a, WGS84_f); geod_direct(&g, slat, slon, sazi, dist, &dlat, &dlon, &dazi); } -// [[Rcpp::export]] -double dir_lonlat(double lon1, double lat1, double lon2, double lat2) { - //double a = 6378137.0; - //double f = 1/298.257223563; - +double direction_geo(double lon1, double lat1, double lon2, double lat2) { double s12, azi1, azi2; struct geod_geodesic g; geod_init(&g, WGS84_a, WGS84_f); @@ -117,15 +103,12 @@ double dir_lonlat(double lon1, double lat1, double lon2, double lat2) { } -// [[Rcpp::export]] -double dir_rad(double lon1, double lat1, double lon2, double lat2) { - - if ((lon1 == lon2) && (lat1 == lat2)) return 0; // NAN? +double direction_cos(double lon1, double lat1, double lon2, double lat2) { + if ((lon1 == lon2) && (lat1 == lat2)) return 0; // NAN? double dLon = lon2 - lon1; double y = sin(dLon) * cos(lat2); double x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon); - double azm = atan2(y, x); azm = fmod(azm+M_PI, M_PI); return azm > M_PI ? -(M_PI - azm) : azm; @@ -133,14 +116,11 @@ double dir_rad(double lon1, double lat1, double lon2, double lat2) { -// [[Rcpp::export]] -double dist2track(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r=6378137) { +double dist2track_geo(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r=6378137) { double a = 1; double f = 0; struct geod_geodesic geod; geod_init(&geod, a, f); - - double d, b2, b3, azi; geod_inverse(&geod, lat1, lon1, lat2, lon2, &d, &b2, &azi); geod_inverse(&geod, lat1, lon1, plat, plon, &d, &b3, &azi); @@ -152,21 +132,16 @@ double dist2track(double lon1, double lat1, double lon2, double lat2, double plo } -// [[Rcpp::export]] -double dist2track_cosine_rad(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r=6378137) { - // same results as above but probably much faster - double b2 = dir_rad(lon1, lat1, lon2, lat2); - double b3 = dir_rad(lon1, lat1, plon, plat); - double d = dist_cosine_rad(lon1, lat1, plon, plat, 1); +inline double dist2track_cos(double lon1, double lat1, double lon2, double lat2, double plon, double plat, bool sign, double r=6378137) { + double b2 = direction_cos(lon1, lat1, lon2, lat2); + double b3 = direction_cos(lon1, lat1, plon, plat); + double d = distance_cos(lon1, lat1, plon, plat, 1); double xtr = asin(sin(b3-b2) * sin(d)) * r; return sign ? xtr : fabs(xtr); } - - -// [[Rcpp::export]] -double alongTrackDistance(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r=6378137) { +double alongTrackDistance_geo(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r=6378137) { double a = 1; double f = 0; struct geod_geodesic geod; @@ -184,13 +159,11 @@ double alongTrackDistance(double lon1, double lat1, double lon2, double lat2, do } +double alongTrackDistance_cos(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r=6378137) { -// [[Rcpp::export]] -double alongTrackDistance_rad(double lon1, double lat1, double lon2, double lat2, double plon, double plat, double r=6378137) { - - double tc = dir_rad(lon1, lat1, lon2, lat2); // * toRad - double tcp = dir_rad(lon1, lat1, plon, plat); // * toRad - double dp = dist_cosine_rad(lon1, lat1, plon, plat, 1); + double tc = direction_cos(lon1, lat1, lon2, lat2); // * toRad + double tcp = direction_cos(lon1, lat1, plon, plat); // * toRad + double dp = distance_cos(lon1, lat1, plon, plat, 1); double xtr = asin(sin(tcp-tc) * sin(dp)); // +1/-1 for ahead/behind [lat1,lon1] @@ -205,49 +178,45 @@ double alongTrackDistance_rad(double lon1, double lat1, double lon2, double lat2 -// [[Rcpp::export]] -double dist2segment(double plon, double plat, double lon1, double lat1, double lon2, double lat2) { - // the alongTrackDistance is the length of the path along the great circle to the point of intersection // there are two, depending on which node you start // we want to use the min, but the max needs to be < segment length - double seglength = dist_lonlat(lon1, lat1, lon2, lat2); - double trackdist1 = alongTrackDistance(lon1, lat1, lon2, lat2, plon, plat); - double trackdist2 = alongTrackDistance(lon2, lat2, lon1, lat1, plon, plat); + +double dist2segment_geo(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double notused=0) { + + double seglength = distance_geo(lon1, lat1, lon2, lat2); + double trackdist1 = alongTrackDistance_geo(lon1, lat1, lon2, lat2, plon, plat); + double trackdist2 = alongTrackDistance_geo(lon2, lat2, lon1, lat1, plon, plat); if ((trackdist1 >= seglength) || (trackdist2 >= seglength)) { - double d1 = dist_lonlat(lon1, lat1, plon, plat); - double d2 = dist_lonlat(lon2, lat2, plon, plat); + double d1 = distance_geo(lon1, lat1, plon, plat); + double d2 = distance_geo(lon2, lat2, plon, plat); return d1 < d2 ? d1 : d2; } - return dist2track(lon1, lat1, lon2, lat2, plon, plat, false); + return dist2track_geo(lon1, lat1, lon2, lat2, plon, plat, false); } -// [[Rcpp::export]] -double dist2segment_cosine_rad(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double r=6378137) { -// the alongTrackDistance is the length of the path along the great circle to the point of intersection -// there are two, depending on which node you start -// we want to use the min, but the max needs to be < segment length - double seglength = dist_cosine_rad(lon1, lat1, lon2, lat2, r); - double trackdist1 = alongTrackDistance_rad(lon1, lat1, lon2, lat2, plon, plat, r); - double trackdist2 = alongTrackDistance_rad(lon2, lat2, lon1, lat1, plon, plat, r); +double dist2segment_cos(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double r=6378137) { + double seglength = distance_cos(lon1, lat1, lon2, lat2, r); + double trackdist1 = alongTrackDistance_cos(lon1, lat1, lon2, lat2, plon, plat, r); + double trackdist2 = alongTrackDistance_cos(lon2, lat2, lon1, lat1, plon, plat, r); if ((trackdist1 >= seglength) || (trackdist2 >= seglength)) { - double d1 = dist_cosine_rad(lon1, lat1, plon, plat, r); - double d2 = dist_cosine_rad(lon2, lat2, plon, plat, r); + double d1 = distance_cos(lon1, lat1, plon, plat, r); + double d2 = distance_cos(lon2, lat2, plon, plat, r); return d1 < d2 ? d1 : d2; } - return dist2track_cosine_rad(lon1, lat1, lon2, lat2, plon, plat, false, r); + return dist2track_cos(lon1, lat1, lon2, lat2, plon, plat, false, r); } // [[Rcpp::export]] -double dist2segmentPoint(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double &ilon, double &ilat) { +double dist2segmentPoint_geo(double plon, double plat, double lon1, double lat1, double lon2, double lat2, double &ilon, double &ilat) { - double seglength = dist_lonlat(lon1, lat1, lon2, lat2); - double trackdist1 = alongTrackDistance(lon1, lat1, lon2, lat2, plon, plat); - double trackdist2 = alongTrackDistance(lon2, lat2, lon1, lat1, plon, plat); + double seglength = distance_geo(lon1, lat1, lon2, lat2); + double trackdist1 = alongTrackDistance_geo(lon1, lat1, lon2, lat2, plon, plat); + double trackdist2 = alongTrackDistance_geo(lon2, lat2, lon1, lat1, plon, plat); if ((trackdist1 >= seglength) || (trackdist2 >= seglength)) { - double d1 = dist_lonlat(lon1, lat1, plon, plat); - double d2 = dist_lonlat(lat2, lat2, plon, plat); + double d1 = distance_geo(lon1, lat1, plon, plat); + double d2 = distance_geo(lat2, lat2, plon, plat); if (d1 < d2) { ilon = lon1; ilat = lat1; @@ -259,19 +228,19 @@ double dist2segmentPoint(double plon, double plat, double lon1, double lat1, dou } } double azi; - double crossd = dist2track(lon1, lat1, lon2, lat2, plon, plat, false); + double crossd = dist2track_geo(lon1, lat1, lon2, lat2, plon, plat, false); if (trackdist1 < trackdist2) { - double bear = dir_lonlat(lon1, lat1, lon2, lat2); - dest_lonlat(lon1, lat1, bear, trackdist1, ilon, ilat, azi); + double bear = direction_geo(lon1, lat1, lon2, lat2); + dest_geo(lon1, lat1, bear, trackdist1, ilon, ilat, azi); } else { - double bear = dir_lonlat(lon2, lat2, lon1, lat1); - dest_lonlat(lon2, lat2, bear, trackdist2, ilon, ilat, azi); + double bear = direction_geo(lon2, lat2, lon1, lat1); + dest_geo(lon2, lat2, bear, trackdist2, ilon, ilat, azi); } - return(crossd); + return crossd; } -std::vector SpatVector::linedistLonLat(SpatVector x, std::string unit) { +std::vector SpatVector::point2lineDistLonLat(SpatVector x, std::string unit, std::string method) { std::vector d; @@ -289,8 +258,16 @@ std::vector SpatVector::linedistLonLat(SpatVector x, std::string unit) { } std::vector> pxy = x.coordinates(); - deg2rad(pxy[0]); - deg2rad(pxy[1]); + + std::function d2seg; + + if (method == "cosine") { + deg2rad(pxy[0]); + deg2rad(pxy[1]); + d2seg = dist2segment_cos; + } else { + d2seg = dist2segment_geo; + } size_t np = pxy[0].size(); size_t ng = size(); @@ -317,7 +294,7 @@ std::vector SpatVector::linedistLonLat(SpatVector x, std::string unit) { } for (size_t j=0; j SpatVector::linedistLonLat(SpatVector x, std::string unit) { if ((d[i] != 0) && (insect[i] == 0)) { for (size_t j=0; j SpatVector::linedistLonLat(SpatVector x, std::string unit) { for (size_t i=0; i. + + +//double distance_cos(const double &lon1, const double &lat1, const double &lon2, const double &lat2, const double &r = 6378137); + +inline double distance_cos(const double &lon1, const double &lat1, const double &lon2, const double &lat2, const double &r = 6378137) { + return r * acos((sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1-lon2))); +} + +double direction_cos(double lon1, double lat1, double lon2, double lat2); diff --git a/src/spatRaster.h b/src/spatRaster.h index 803e47230..1a66e8cf7 100644 --- a/src/spatRaster.h +++ b/src/spatRaster.h @@ -647,13 +647,13 @@ class SpatRaster { SpatRaster fillNA(double missing, double maxdist, int niter, SpatOptions &opt); SpatRaster distance_rasterize(SpatVector p, double target, double exclude, std::string unit, const std::string& method, SpatOptions &opt); - SpatRaster direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, SpatOptions &opt); + SpatRaster direction_rasterize(SpatVector p, bool from, bool degrees, double target, double exclude, const std::string& method, SpatOptions &opt); SpatRaster distance_spatvector(SpatVector p, std::string unit, 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 direction(bool from, bool degrees, double target, double exclude, SpatOptions &opt); - SpatRaster direction_vector(SpatVector p, bool from, bool degrees, SpatOptions &opt); + SpatRaster direction(bool from, bool degrees, double target, double exclude, const std::string& method, SpatOptions &opt); + SpatRaster direction_vector(SpatVector p, bool from, bool degrees, const std::string& method, SpatOptions &opt); SpatRaster clumps(int directions, bool zeroAsNA, SpatOptions &opt); SpatRaster patches(size_t directions, SpatOptions &opt); diff --git a/src/spatVector.h b/src/spatVector.h index 9a77a607f..0278db76f 100644 --- a/src/spatVector.h +++ b/src/spatVector.h @@ -203,7 +203,7 @@ class SpatVector { // std::vector pointdistance_seq(const std::vector& px, const std::vector& py, double m, bool lonlat); - std::vector linedistLonLat(SpatVector x, std::string unit); + std::vector point2lineDistLonLat(SpatVector x, std::string unit, std::string method); std::vector> knearest(size_t k);