Sf: Feature Request: Mapbox Cheap Ruler

[mem48]

The mapbox cheap ruler is a way to approximate distances between lng/lat points much faster than the normal geodesic formulas. For short distances (less than 100 miles) it is accurate to 0.1% or better but is x30 faster.

For many applications, this is good enough and it would be useful to have this option in st_distance() and related functions.

See some related content:
https://github.com/ATFutures/dodgr/issues/57
https://cran.r-project.org/web/packages/geodist/vignettes/geodist.html

[edzer]

As nobody seems to pick this up, I'm closing here.

[mem48]

😪

[mem48]

@mpadge is there any way to implement this using geodist? Perhaps as a reccomended package for sf?

[mpadge]

cheapdist is the default method in geodist, but the standard reference still remains Karney's geodesic algorithms used in sf and bundled with geodist1. sf is intended to serve as a "reference" package for such calculations, and is right to stick with geodesic alone, I'd suggest. The main aim of geodist is to offer as many alternatives as possible, primarily for usages where efficiency is judged more important than absolute accuracy. sf is and should remain the benchmark for accuracy, so i personally would see it as scope-creep for it to start trading that off in key areas for the sake of efficiency gains. Also not doing it leaves scope for other packages and developers to implement and offer more efficient alternatives for those for whom efficient is important.

[mem48]

Valid points, perhaps then the solution is to mirror sf functions in geodists or another package. For example how could I get a cheap version of st_area()?

[mpadge]

Best answer i can give for the moment is to watch hypertidy developments and wait. We've got really efficient implementations of the clipper algorithm, adapted to efficiently calculate multiple areas, but currently for internal use and not yet exported. The main algorithm has been wrapped in the R package polyclip, but that version only exposes a limited set of the functionality, and not the area algorithm. And importantly, this is again scope creep for sf, because clipper always assumes and is applicable to planar geometry only. It can be a lot faster than sf, but is only applicable to this relatively very restricted special case.

The unfortunate thing here is that the most desirable way to proceed would be for the full functionality of clipper to be exposed by the polyclip package, which has rightly claimed this domain, but the github issue requesting such extension has lain dormant since Sep 2017. I'll chat to @mdsumner about this and see what he thinks, and whether we might be able to pick up the slack and offer something along these lines.

[Robinlovelace]

Another comment on this, seems that someone's come up with a fast 'cheap ruler' in Rust: https://github.com/Turbo87/flat-projection-rs

Rust can be called from R with this approach: https://github.com/r-rust/hellorust

[edzer]

sf spawns into lwgeom::lwgeom_distance, which uses the liblwgeom code, and that uses a "cheap" spherical distance if the datum of x and y is spherical; that should be pretty much faster than the default spheroid/ellipsoid computations; a simple benchmark with a made-up sphere:

wkt = 'GEOGCRS["spherical",
    DATUM["MySphere",
        ELLIPSOID["WGS 84Sphere", 6378137, 0,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]]
]'

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
sph = st_crs(wkt)

pt1 = st_point(c(0,0))
pt2 = st_point(c(10,10))

p1 = st_sfc(pt1, pt2, crs = sph)
p2 = st_sfc(pt1, pt2, crs = 4326)
s = sample(1:2, 5000, replace = TRUE)

system.time(st_distance(p1[s]))
#    user  system elapsed 
#   4.719   0.064   4.783 
system.time(st_distance(p2[s]))
#    user  system elapsed 
#  24.395   0.043  24.445 

[Robinlovelace]

Neat solution and doesn't involve writing another R package. Don't know how it compares with the other options on this thread but a 5 times speedup looks good to me!

[mpadge]

I didn't know about that, thanks @edzer. This is an addendum just for the sake of completeness:

# <all preceding code, plus>
library(geodist)
xy <- st_coordinates(p1)

rbenchmark::benchmark (
                       st_distance(p1[s]),
                       st_distance(p2[s]),
                       geodist (xy[s, ]),
                       replications = 1)
#>                 test replications elapsed relative user.self sys.self
#> 3   geodist(xy[s, ])            1   0.201    1.000     0.085    0.117
#> 1 st_distance(p1[s])            1   3.961   19.706     3.885    0.073
#> 2 st_distance(p2[s])            1  20.327  101.129    20.213    0.086

^{Created on 2020-04-09 by the reprex package (v0.3.0)}

[mem48]

Does this work for other types of function e.g st_area