Points closest to the calculated path

Question

As in the previous question: The closest point to the path

I would like to find all the centers closest to the path. However, some data are missing on the path, and I would like to make linear segments in order to interpolate between the points to “estimate” the possible path and still find probable centers that would be close to this “calculated path”.

set.seed(1) n <- 10000 x <- 100*cumprod(1 + rnorm(n, 0.0001, 0.002)) y <- 50*cumprod(1 + rnorm(n, 0.0001, 0.002)) # original path path <- data.frame(cbind(x=x, y=y)) # path with missing points/points every hundred path.w.missing <- path[seq(1,n,by=100),] centers <- expand.grid(x=seq(0, 500,by=0.5) + rnorm(1001), y=seq(0, 500, by=0.2) + rnorm(2501)) centers$id <- seq(nrow(centers)) 

Except for modeling millions of linear points between given waypoints ... I'm not sure how this could be done ...

It sounds a bit like finding the intersection of a line and a cell matrix ... sort ... but maybe I'm miles from ...

Any help would be greatly appreciated.

Answer 1

Did it help:

 ## create data set.seed(1) n <- 10000 x <- 100 * cumprod(1 + rnorm(n, 0.0001, 0.002)) y <- 50 * cumprod(1 + rnorm(n, 0.0001, 0.002)) # original path path <- data.frame(cbind(x=x, y=y)) # path with missing points/points every hundred path.w.missing <- path[seq(1,n,by=100),] ## proposed solution library(raster) l <- SpatialLines(list(Lines(list(Line(path.w.missing)), "1"))) r <- raster(extent(bbox(l)), res = c(0.5, 0.2)) r <- rasterize(l, r, field = 1) ## xy will be a matrix xy <- rasterToPoints(r) plot(r) points(xy) points(path.w.missing, col = "red") 

Answer 2

The marmap package does exactly what you need:

 library(marmap) bcenters <- as.bathy(centers) out <- path.profile(path.w.missing,bcenters) lon lat dist.km depth 2 100.0512 50.04129 0.00000 247451 3 101.2070 50.11430 82.87177 247450 4 101.3687 50.18730 95.33454 247449 5 101.5973 50.26030 112.81665 248453 6 102.6877 50.33330 190.48139 248454 7 103.2105 50.40630 228.37510 248458 

id is in the depth column. This function was encoded to extract depth in the bathymetric matrix along the path. It finds all the cells of the bathymetric matrix “below” the path (ie, the cells closest to the path) and displays the value of these cells. By converting your centers object to a bathy matrix, your id column is used as if it were depth. Thus, path.profile() displays the id cells closest to your path in the depth column.

Here you start with 100 points in your path.w.missing object. path.profile() gets the identifier of

 nrow(out) 795 

which are closest to this path.

Here is a chart of data path.w.missing from the nearest points, superimposed in red:

 plot(path.w.missing,type="o",cex=.3, lwd=5) points(out[,1],out[,2], col=rgb(1,0,0,.8), cex=.3) 

Answer 3

Use the smooth.spline function (set the smoothing spline)

 mx=path.w.missing$x my=path.w.missing$y s<-smooth.spline(mx,my) plot(mx,my) lines(s) #draw the estimated path 

 predict(s,1:10) #And now you can generate all the points on the estimated path > predict(s,1:10) $x [1] 1 2 3 4 5 6 7 8 9 10 $y [1] 2.418294 2.904019 3.389744 3.875469 4.361195 4.846920 5.332645 5.818370 6.304095 6.789820 

And you can use a linear model to match it.

 > summary(lm(my~mx)) Call: lm(formula = my ~ mx) Residuals: Min 1Q Median 3Q Max -12.772 -8.642 3.112 5.831 17.237 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -12.60098 3.46583 -3.636 0.000444 *** mx 0.56579 0.02138 26.469 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 7.896 on 98 degrees of freedom Multiple R-squared: 0.8773, Adjusted R-squared: 0.876 F-statistic: 700.6 on 1 and 98 DF, p-value: < 2.2e-16 

as you can see, p <0.05, so my and mx can fit as a string:

 my=-12.60098+0.56579*mx 

ggplot2 can easily execute this line:

  d=data.frame(mx,my) library(ggplot2) ggplot(d,aes(mx,my))+geom_point()+geom_smooth(method="lm") 

Suppose that the line Ax+By+C=0 , and the point (X0，Y0) , indicate the distance of the line:

 |AX0+BY0+C|/√(A²+B²) 

therefore, in this case 0.56579 * mx-my-12.60098 = 0, A = 0.56579, B = -1, C = -12.60098, it is easy to calculate the distance from point to line and find the nearest point to the line.

In addition, if you want to find the nearest point, remove the denominator √ (A² + B²), this will not affect the sorting, therefore optimizing the formula:

 |AX0+BY0+C| 

Result

 > for(i in 1:2503501){ + temp=abs(centers[[1]][i]*0.56579-centers[[2]][i]-12.60098) + if(m>temp){ + m=temp + pos=i + } + } > m [1] 2.523392e-05 > pos [1] 638133 

Use rcpp to speed up the test.cpp program

 #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] int closest(NumericVector a, NumericVector b) { int length = a.size(); double temp,m=10000; int pos=0; for(int i=0;i<length;i++){ temp=a[i]*0.56579-b[i]-12.60098; if(temp>=0){ if(m>temp){m=temp;pos=i;} }else{ if(m>-temp){m=-temp;pos=i;} } } return pos+1; } 

Result (Rcpp costs about 2 seconds):

 > sourceCpp("test.cpp") > closest(centers[[1]],centers[[2]]) [1] 974698 > abs(centers[[1]][974698]*0.56579-centers[[2]][974698]-12.60098) [1] 0.0002597022