st_geod_segmentize

[edzer]

In order to be able to drop dependency on lwgeom entirely, sf needs a replacement of lwgeom::st_geod_segmentize, to use in sf::st_segmentize for geographic data. S2Polylines have an Interpolate method that could accomplish this. A tricky thing is that segmentize adds points, at regular intervals bound by a maximum distance, to the existing points of a linestring.

[paleolimbot]

I've seen simplifiers, but no segmentizers yet! I'll keep an eye out.

[paleolimbot]

There's a reference to s2builderutil::EdgeSplitter() in the docs, but this doesn't seem to exist. Maybe it's the EdgeTesselator?

    // Self-intersections can also arise when importing data from a 2D
    // projection.  You can minimize this problem by subdividing the input
    // edges so that the S2 edges (which are geodesics) stay close to the
    // original projected edges (which are curves on the sphere).  This can
    // be done using s2builderutil::EdgeSplitter(), for example.

[edzer]

That seems to be dealing with crossing edges, and I also cannot find any EdgeSplitter in the code.

Although s2_segmentize would be nice to have, it is not that urgent since lwgeom::st_geod_segmentize now does the same thing; also, this is mainly useful for plotting after projecting, when great circles have to become curves, and that is both plotting and projecting are done outside s2. I won't drop lwgeom straight away but give the option to use s2 in sf, and switch back to lwgeom, in any case for some period of time. See https://github.com/r-spatial/sf/blob/s2/vignettes/sf7.Rmd (WIP) for a write-up (I added you as co-author). Next 0.9-x release of sf will not use s2 by default but can be switched on, sf 1.0-x will use it as default.

For segmentizing in S2 closest I have seen is S2Polyline::Interpolate(), but this wouldn't work for polygons, where it would also be needed.

[paleolimbot]

It's a bit more "manual", but one could also loop through edges and subdivide them...the S2Builder interface provides an AddEdge() method, and one can loop through shapes in a shape index (and edges in a shape) to make this happen. Battle for another day, but definitely possible.

[paleolimbot]

Also, the blog post looks great!

[paleolimbot]

I'll have a closer look at how this is done in GEOS, but because a S2ShapeIndex is just a bunch of edges, I think that one can "just" loop through edges in the index and replace each edge with n new edges if needed and then pass to the builder to take care of the geometry generation.

Alternatively, one could just work with the vertices (something like export to WKB, split edges, then re-import), although my sense is that working with the edges will result in cleaner code (even though it will take a while to figure out).

[paleolimbot]

As far as I can tell there's no segmentizer built-in...the code is mostly focused on removing as many vertices as possible whenever possible. The algorithm isn't complex...given a 0...1 vector t_norm01 it's just the weighted average of the xyz (S2Point) representation of the point.

  interp_xyz <- matrix(
    p2_xyz * rep(t_norm01, each = 3) +
      p1_xyz * rep(1 - t_norm01, each = 3),
    ncol = 3,
    byrow = TRUE
  )

[paleolimbot]

Going the other direction (projected -> spherical), this is where the EdgeTesselator comes in. It's very similar to D3's adaptive resample algorithm...adding points opportunistically based on some tolereance. Both are similar in that they work with edges not geographies.

[paleolimbot]

This is now possible after #115. It doesn't make much sense to have an s2_tessellate() function because an s2_geography() vector already has implicit geodesic edges...it's the import and export that need this option. With the current implementation, one could implement an sf operator (see below). Is an explicit segmentize needed, or is tessellation enough?

library(s2)
library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
#> Linking to GEOS 3.8.1, GDAL 3.2.0, PROJ 7.2.0

tessellate_geod <- function(x, distance, radius = s2_earth_radius_meters()) {
  coords_tes <- wk::wk_handle(
    x,
    s2_unprojection_filter(
      s2_projection_filter(
        wk::sfc_writer(),
        tessellate_tol = distance / radius
      )
    )
  )
  
  coords_tes %>% st_set_crs(st_crs(x))
}

plot(tessellate_geod(sf::st_as_sfc("LINESTRING (0 0, 0 45, -60 45)"), 100))

^{Created on 2021-06-05 by the reprex package (v0.3.0)}

[edzer]

Trying to get this into sf; I'm not getting what I do wrong here:

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
sf = st_sf(a=1, geom=st_sfc(st_linestring(rbind(c(0,0),c(1,1)))), crs = 4326)
library(s2)

tessellate_geod <- function(x, distance, radius = s2::s2_earth_radius_meters()) {
    if (!requireNamespace("wk", quietly = TRUE))
        stop("package wk not available: install it first?")
    tol = if (units::ud_are_convertible(units(distance), "rad"))
            units::set_units(distance, "rad", mode = "standard")
        else
            units::set_units(distance, "m", mode = "standard") / radius
    coords_tes <- wk::wk_handle(
        x,
        s2::s2_unprojection_filter(
            s2::s2_projection_filter(
                wk::sfc_writer(),
                tessellate_tol = units::drop_units(tol)
            )
        )
    )
    st_set_crs(coords_tes, st_crs(x))
}

(seg = tessellate_geod(sf, units::set_units(100, km)))
# Geometry set for 1 feature 
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
# Geodetic CRS:  WGS 84
# LINESTRING (0 0, 1 1)

[edzer]

sorry for the noise, fixed.

[edzer]

I now understand what this does, and it seems useful, for instance when plotting, but not a replacement / drop-in for st_segmentize. I'll leave the lwgeom stuff in place, there.