Geodesic segmentization via PROJ.

GeometryOps.segmentize Function

segmentize([method = Planar()], geom; max_distance::Real, threaded)

Segmentize a geometry by adding extra vertices to the geometry so that no segment is longer than a given distance. This is useful for plotting geometries with a limited number of vertices, or for ensuring that a geometry is not too "coarse" for a given application.

Arguments

• method::Manifold = Planar(): The method to use for segmentizing the geometry. At the moment, Planar (assumes a flat plane), Spherical (assumes geometry on a sphere and interpolates along great circles) and Geodesic (assumes geometry on the ellipsoidal Earth and uses Vincenty's formulae) are available.

• geom: The geometry to segmentize. Must be a LineString, LinearRing, Polygon, MultiPolygon, or GeometryCollection, or some vector or table of those.

• max_distance::Real: The maximum distance between vertices in the geometry. For Planar manifolds, this is in the units of the geometry. For Spherical and Geodesic, this is in the units of the manifold (the units of the radius of the sphere/ellipsoid). By default this is meters.

Spherical and Geodesic both assume that the input geometry is in lon/lat coordinates, in degrees.

Returns a geometry of similar type to the input geometry, but resampled.

source

Implementation

The implementation uses PROJ's geodesic calculations to: 2. Compute the geodesic distance between points

3. Calculate intermediate points along the great circle path

4. Add points at regular intervals based on the maximum distance parameter

Key features:

• Uses PROJ's geod_geodesic for accurate ellipsoidal calculations

• Configurable equatorial radius and flattening

• Thread-safe implementation

• Supports both LineString and LinearRing geometries

The function creates a geodesic line between each pair of points and interpolates positions along that line at the specified maximum distance intervals.

This holds the segmentize geodesic functionality.

import GeometryOps: GeodesicSegments, _segmentize, _fill_linear_kernel!
import Proj

function GeometryOps.GeodesicSegments(; max_distance, equatorial_radius::Real=6378137, flattening::Real=1/298.257223563, geodesic::Proj.geod_geodesic = Proj.geod_geodesic(equatorial_radius, flattening))
    return GeometryOps.GeodesicSegments{Proj.geod_geodesic}(geodesic, max_distance)
end

This is the same method as in transformations/segmentize.jl, but it constructs a Proj geodesic line every time. Maybe this should be better...

function _segmentize(method::Geodesic, geom, ::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance)
    proj_geodesic = Proj.geod_geodesic(method.semimajor_axis #= same thing as equatorial radius =#, 1/method.inv_flattening)
    first_coord = GI.getpoint(geom, 1)
    x1, y1 = GI.x(first_coord), GI.y(first_coord)
    new_coords = NTuple{2, Float64}[]
    sizehint!(new_coords, GI.npoint(geom))
    push!(new_coords, (x1, y1))
    for coord in Iterators.drop(GI.getpoint(geom), 1)
        x2, y2 = GI.x(coord), GI.y(coord)
        _fill_linear_kernel!(method, new_coords, x1, y1, x2, y2; max_distance, proj_geodesic)
        x1, y1 = x2, y2
    end
    return rebuild(geom, new_coords)
end

function GeometryOps._fill_linear_kernel!(method::Geodesic, new_coords::Vector, x1, y1, x2, y2; max_distance, proj_geodesic)
    geod_line = Proj.geod_inverseline(proj_geodesic, y1, x1, y2, x2)

This is the distance in meters computed between the two points. It's s13 because geod_inverseline sets point 3 to the second input point.

    distance = geod_line.s13
    if distance > max_distance
        n_segments = ceil(Int, distance / max_distance)
        for i in 1:(n_segments - 1)
            y, x, _ = Proj.geod_position(geod_line, i / n_segments * distance)
            push!(new_coords, (x, y))
        end
    end

End the line with the original coordinate, to avoid any multiplication errors.

    push!(new_coords, (x2, y2))
    return nothing
end

---

This page was generated using Literate.jl.