Skip to content

This holds the segmentize geodesic functionality.

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

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

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

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

This page was generated using Literate.jl.