R Under development (unstable) (2023-12-13 r85679 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > suppressPackageStartupMessages(library(sf)) > suppressPackageStartupMessages(library(units)) > > x = st_sfc( + st_point(c(0,0)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + crs = 4326 + ) > > y = st_sfc( + st_point(c(0,10)), + st_point(c(1,0)), + st_point(c(2,0)), + st_point(c(3,0)), + st_point(c(4,0)), + crs = 4326 + ) > > if (suppressPackageStartupMessages(require(sp, quietly = TRUE))) { + d.sf = st_distance(x, y) + d.sp = spDists(as(x, "Spatial"), as(y, "Spatial")) + units(d.sp) = as_units("km") + print(round(d.sf - d.sp, 7)) + + #summary(unclass(d.sf) - d.sp) + + st_crs(x) = st_crs(y) = NA + d.sf = st_distance(x, y) + d.sp = spDists(as(x, "Spatial"), as(y, "Spatial")) + print(round(d.sf - d.sp, 7)) + } Units: [m] [,1] [,2] [,3] [,4] [,5] [1,] 6107.765 -124.3896 -248.7792 -373.1688 -497.5585 [2,] 6065.138 0.0000 -124.3896 -248.7792 -373.1688 [3,] 5940.569 -124.3896 0.0000 -124.3896 -248.7792 [4,] 5743.252 -248.7792 -124.3896 0.0000 -124.3896 [,1] [,2] [,3] [,4] [,5] [1,] 0 0 0 0 0 [2,] 0 0 0 0 0 [3,] 0 0 0 0 0 [4,] 0 0 0 0 0 > > # st_length: > st_crs(y) = 4326 > (z = st_sfc(st_linestring(rbind(c(0,10), c(1,0), c(2,0), c(3,0), c(4,0))), crs = 4326)) Geometry set for 1 feature Geometry type: LINESTRING Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 4 ymax: 10 Geodetic CRS: WGS 84 LINESTRING (0 10, 1 0, 2 0, 3 0, 4 0) > d = st_distance(y, y) > round(d, 7) Units: [m] [,1] [,2] [,3] [,4] [,5] [1,] 0 1117440.6 1133750.1 1160423.1 1196767.0 [2,] 1117441 0.0 111195.1 222390.2 333585.3 [3,] 1133750 111195.1 0.0 111195.1 222390.2 [4,] 1160423 222390.2 111195.1 0.0 111195.1 [5,] 1196767 333585.3 222390.2 111195.1 0.0 > st_length(z) 1451026 [m] > st_length(z) - sum(d[1,2], d[2,3], d[3,4], d[4,5]) 0 [m] > > # st_line_sample: > ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), + st_linestring(rbind(c(0,0),c(10,0)))) > # set.seed(135) > st_line_sample(ls, density = 1) Geometry set for 2 features Geometry type: MULTIPOINT Dimension: XY Bounding box: xmin: 0 ymin: 0 xmax: 9.5 ymax: 0.5 CRS: NA MULTIPOINT ((0 0.5)) MULTIPOINT ((0.5 0), (1.5 0), (2.5 0), (3.5 0),... > > ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))), + st_linestring(rbind(c(0,0),c(.1,0))), crs = 4326) > > st_length(ls) Units: [m] [1] 111195.10 11119.51 > try(st_line_sample(ls, density = 1/1000)) Error in st_line_sample(ls, density = 1/1000) : st_line_sample for longitude/latitude not supported; use st_segmentize? > x = st_line_sample(st_transform(ls, 3857), density = 1/1000) # one per km > > proc.time() user system elapsed 1.06 0.20 1.25