Converting line to set of H3 hexagons?

Question

What is a proper way to convert a GeoJSON linestring to a set of H3 hexagons?

I came up with an algorithm where I

1. extract points from a line,
2. calculate an H3 hexagon of these points geo_to_h3
3. connect the pairs of H3 hexagons h3-line
4. concatenate the pair paths removing duplicating hexagons

I wonder if there is a better way to do the conversion?

Answer 1

I think this is fine, with two caveats:

• You should preserve order by using an ordered data structure (and not an unordered set).
• You should allow for some repetitions by using a data structure that allows for some repeated DGGS cells: this allows for the representation of self-intersecting lines. You can do this by eliminating only sequential duplicates (which are redundant, unless perhaps you wanted to record the original vertex density, which is possible).

I have an implementation of this in Python:

from typing import Iterator, Union

import h3
from shapely.geometry.linestring import LineString
from shapely.geometry.multilinestring import MultiLineString

def sequential_deduplication(func: Iterator[str]) -> Iterator[str]:
    '''
    Decorator that doesn't permit two consecutive items to be the same
    '''
    def inner(*args):
        iterable = func(*args)
        last = None
        while (cell := next(iterable, None)) is not None:
            if cell != last:
                yield cell
            last = cell
    return inner

@sequential_deduplication
def h3polyline(line: Union[LineString, MultiLineString], resolution: int) -> Iterator[str]:
    '''
    Iterator yeilding H3 cells representing a (multi)line,
    retaining order and self-intersections
    '''
    if line.geom_type == 'MultiLineString':
        # Recurse after getting component linestrings from the multiline
        for l in map(lambda geom: h3polyline(geom, resolution), line.geoms):
            yield from l
    else:
        coords = zip(line.coords, line.coords[1:])
        while (vertex_pair := next(coords, None)) is not None:
            i, j = vertex_pair
            a = h3.geo_to_h3(*i[::-1], resolution)
            b = h3.geo_to_h3(*j[::-1], resolution)
            yield from h3.h3_line(a, b) # inclusive of a and b

Although it accepts multilinestring input, it doesn't really retain the distinction of multiple parts (I just didn't need that distinction).

That said, you might not care about order or self-intersections. This really depends on the intended application of your H3-discretised representation of a line. For lines with really long distances between vertices, perhaps first projecting and densifying a line would be sensible.

Answer 2

Same problem, similar solution. Only one additional step that seemed to provide a smoother / more accurate result (i.e., H3s that more closely follow the shape of the line).

For each set of points along the line, create a "sub-line", take the centerpoint of that sub-line, and run h3-line from start --> centerpoint and centerpoint --> end.

More detail for clarity:

1. Extract points from linestring
2. Calculate H3 for starting point
3. Calculate H3 for "next" point (point+1) -- note: make sure not to end up going back to the first point when you point+1 on the last point in the linestring
4. Create "line" between starting point and next point
5. Get centroid of line
6. Calculate H3 for centroid
7. Use h3-line for start-->mid & mid-->next
8. Concatenate / remove duplicate h3s

Answer 3

I want to share a discussion from a Github issue. I think it is valuable regarding the question.

Nick Rabinowitz said

Your approach is what I would have recommended. h3_line has some limitations here - it's describing a path between two cells, so it may not overlap perfectly with the line it's representing. The only guarantees it offers are 1) the path it returns is the shortest path (though there may be many shortest paths of equal length) and 2) each successive cell touches the previous cells but none of the cells in the line before that. The shortest path guarantee usually means that h3_line follows great arcs rather than cartesian lines, which may or may not be what you need.

If you want to make sure that the line is covered by H3 cells, you might try e.g. sampling indexes along the length of the line, buffering them by 1, and then removing indexes that the line does not intersect. This is a much more expensive operation than h3_line, but it would ensure that the line is fully covered.

Kevin Sahr added

Just to clarify, within each icosahedron face the grid shortest path will follow a great circle arc, which maps to a cartesian line in the H3 projection. If you want to precisely capture a great circle arc/line that crosses icosahedron edges, you would need to break-up the original arc/line wherever it intersects an edge, and then perform your algorithm on each of the resulting line segments.

Answer 4

Solution above proposed by @alphabetasoup works very well with H3 version 3.x

With H3 version 4.x you need to follow the migration documentation and make the following replacements:

• replace geo_to_h3 with latlng_to_cell.
• replace h3_line with grid_path_cells.