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, only Planar (assumes a flat plane) 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. Beware: for Planar, this is in the units of the geometry, but for Geodesic and Spherical it's in units of the radius of the sphere.

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

source

Implementation

The implementation uses PROJ's geodesic calculations to:

1. Compute the geodesic distance between points

2. Calculate intermediate points along the great circle path

3. 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.