Skip to content

Implements min_number for segmentation #406

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
runarberg wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Implements min_number for segmentation #406

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 109 additions & 14 deletions src/transformations/segmentize.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,12 @@ export segmentize
export LinearSegments, GeodesicSegments

#=
This function "segmentizes" or "densifies" 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.
This function "segmentizes" or "densifies" a geometry by adding extra
vertices to the geometry so that either a) no segment is longer than a
given distance, or b) total number of vertices is no less then a given
number. 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.

!!! info
We plan to add interpolated segmentization from DataInterpolations.jl in the future,
Expand All @@ -26,6 +27,16 @@ collect(GI.getpoint(linear))
```
You can see that this geometry was segmentized correctly, and now has 8 vertices where it previously had only 4.

If you want to _increase_ the number of vertices above to a specific number you can use the `max_number` keyword
argument instead of `min_distance`.

```@example segmentize
rectangle = GI.Wrappers.Polygon([[(0.0, 50.0), (7.071, 57.07), (0, 64.14), (-7.07, 57.07), (0.0, 50.0)]])
linear = GO.segmentize(rectangle; min_number = 9)
collect(GI.getpoint(linear))
```
This gives you the exact same segmentation as above.

Now, we'll also segmentize this using the geodesic method, which is more accurate for lat/lon coordinates.

```@example segmentize
Expand Down Expand Up @@ -163,32 +174,57 @@ end
# ## Implementation

"""
segmentize([method = Planar()], geom; max_distance::Real, threaded)
segmentize([method = Planar()], geom; max_distance::Union{Real, Nothing}, min_number::Union{Int, Nothing}, threaded)

Segmentize a geometry by adding extra vertices to the geometry so that no segment is longer than a given distance.
Segmentize a geometry by adding extra vertices to the geometry so that no segment is longer than a given distance or contains at least a given number of vertices.
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`](@ref) (assumes a flat plane) and [`Geodesic`](@ref) (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.**
- `max_distance::Union{Real, Nothing}`: 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.**
- `min_number::Union{Int, Nothing}`: The minimum number of vertices in the geometry.

Returns a geometry of similar type to the input geometry, but resampled.
"""
function segmentize(geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = False())
return segmentize(Planar(), geom; max_distance, threaded = booltype(threaded))
function segmentize(
geom;
max_distance::Union{Real, Nothing} = nothing,
min_number::Union{Int, Nothing} = nothing,
threaded::Union{Bool, BoolsAsTypes} = False()
)
return segmentize(Planar(), geom; max_distance, min_number, threaded = booltype(threaded))
end

# allow three-arg method as well, just in case
segmentize(geom, max_distance::Real; threaded = False()) = segmentize(Planar(), geom, max_distance; threaded)
segmentize(method::Manifold, geom, max_distance::Real; threaded = False()) = segmentize(Planar(), geom; max_distance, threaded)

# generic implementation
function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = False())
if max_distance <= 0
throw(ArgumentError("`max_distance` should be positive and nonzero! Found $(max_distance)."))
function segmentize(
method::Manifold,
geom;
max_distance::Union{Real,Nothing} = nothing,
min_number::Union{Int,Nothing} = nothing,
threaded::Union{Bool, BoolsAsTypes} = False(),
)
if count(isnothing, (max_distance, min_number)) ≠ 1
throw(ArgumentError("Must provide one of `max_distance`, or `min_number` keywords"))
end
if !isnothing(max_distance)
if max_distance ≤ 0
throw(ArgumentError("`max_distance` should be positive and nonzero! Found $(max_distance)."))
end
elseif !isnothing(min_number)
if min_number ≤ 0
throw(ArgumentError("`min_number` should be positive and nonzero! Found $(min_number)."))
end
end
_segmentize_function(geom) = if !isnothing(max_distance)
_segmentize(method, geom, GI.trait(geom); max_distance)
else
_segmentize(geom, GI.trait(geom); min_number)
end
_segmentize_function(geom) = _segmentize(method, geom, GI.trait(geom); max_distance)
return apply(_segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded)
end

Expand Down Expand Up @@ -217,6 +253,65 @@ function _segmentize(method::Union{Planar, Spherical}, geom, T::Union{GI.LineStr
return rebuild(geom, new_coords)
end

function _segmentize(geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; min_number)
new_coords = GI.getgeom(geom)

if min_number <= length(new_coords)
return rebuild(geom, new_coords)
end

distances = Vector()
first_coord = GI.getpoint(geom, 1)
x1, y1 = GI.x(first_coord), GI.y(first_coord)
for coord in Iterators.drop(GI.getpoint(geom), 1)
x2, y2 = GI.x(coord), GI.y(coord)
dx, dy = x2 - x1, y2 - y1
distance = hypot(dx, dy)
push!(distances, distance)
x1, y1 = x2, y2
end

if new_coords isa Vector{Tuple{Int, Int}}
new_coords = convert(Vector{Tuple{Float64, Float64}}, new_coords)
elseif new_coords isa Vector{Vector{Int}}
new_coords = convert(Vector{Vector{Float64}}, new_coords)
end

point = if first_coord isa Tuple
(a, b) -> (a, b)
else
(a, b) -> [a, b]
end

while length(new_coords) < min_number
max_2 = 0
i, max_1 = 0, 0
for (j, dist) in Iterators.enumerate(distances)
if dist > max_1
max_2 = max_1
i, max_1 = j, dist
elseif dist > max_2
max_2 = dist
end
end
x1, y1 = new_coords[i]
x2, y2 = new_coords[i + 1]
dx, dy = x2 - x1, y2 - y1

n_segments = 1 + min(ceil(Int, max_1 / (max_2 + 1)), min_number - length(new_coords))
new_dist = max_1 / n_segments
distances[i] = new_dist
for j in 1:(n_segments - 1)
t = j / n_segments
new_coord = point(x1 + t * dx, y1 + t * dy)
insert!(new_coords, i + j, new_coord)
insert!(distances, i + j, new_dist)
end
end

return rebuild(geom, new_coords)
end

function _fill_linear_kernel!(::Planar, new_coords::Vector, x1, y1, x2, y2; max_distance)
dx, dy = x2 - x1, y2 - y1
distance = hypot(dx, dy) # this is a more stable way to compute the Euclidean distance
Expand Down
81 changes: 80 additions & 1 deletion test/transformations/segmentize.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ using ..TestHelpers
mp = GI.MultiPolygon([p, p, p])
mls = GI.MultiLineString([ls, ls, ls])

@testset_implementations "LinearSegments" begin
@testset_implementations "LinearSegments max_distance" begin
@test GO.segmentize($ls; max_distance = 0.5) isa GI.LineString
if GI.trait($lr) isa GI.LinearRingTrait
@test GO.segmentize($lr; max_distance = 0.5) isa GI.LinearRing
Expand All @@ -30,6 +30,25 @@ using ..TestHelpers
@test GI.ngeom(GO.segmentize($mls; max_distance = 0.5)) == 3
end

@testset_implementations "LinearSegments min_number" begin
@test GO.segmentize($ls; min_number = 7) isa GI.LineString
if GI.trait($lr) isa GI.LinearRingTrait
@test GO.segmentize($lr; min_number = 7) isa GI.LinearRing
end
# Test that linear rings are closed after segmentization
segmentized_ring = GO.segmentize($lr; min_number = 7)
@test GI.getpoint(segmentized_ring, 1) ==
GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring))

@test GO.segmentize($p; min_number = 7) isa GI.Polygon
@test GO.segmentize($mp; min_number = 7) isa GI.MultiPolygon
@test GI.ngeom(GO.segmentize($mp; min_number = 7)) == 3

# Now test multilinestrings
@test GO.segmentize($mls; min_number = 7) isa GI.MultiLineString
@test GI.ngeom(GO.segmentize($mls; min_number = 7)) == 3
end

@testset_implementations "GeodesicSegments" begin
@test GO.segmentize(GO.Geodesic(), $ls; max_distance = 0.5*900) isa GI.LineString
if GI.trait($lr) isa GI.LinearRingTrait
Expand Down Expand Up @@ -60,9 +79,69 @@ lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
end
end

@testset_implementations "Planar min_number" begin
ct = GO.centroid($lr)
ar = GO.area($lr)
for min_number in 5:15
segmentized = GO.segmentize(GO.Planar(), $lr; min_number)
@test all(GO.centroid(segmentized) .≈ ct)
@test GO.area(segmentized) ≈ ar
end
end

lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
@testset_implementations "Geodesic" begin
for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) .* 900
@test_nowarn segmentized = GO.segmentize(GO.Geodesic(), $lr; max_distance)
end
end

@testset "Segmentation with min_number" begin
@testset_implementations "Simple line" begin
# If there is no line to be the second longest we risk going
# into an infinate loop.
ls = GI.LineString([(0, 0), (10, 2)])
segmentized = GO.segmentize(ls; min_number = 3)
@test GI.ngeom(segmentized) == 3
@test GI.getgeom(segmentized, 2) == (5.0, 1.0)
end

@testset_implementations "Geom is larger then min_number" begin
ls = GI.LineString([(0, 0), (5, 5), (10, 10)])
segmentized = GO.segmentize(ls; min_number = 2)
@test GI.ngeom(segmentized) == 3
@test GI.getgeom(segmentized, 2) == (5, 5)
end

@testset_implementations "Cuts long segments in n-parts when second longest segment allows" begin
ls = GI.LineString([(0.0, 0.0), (7.5, 7.5), (10.0, 10.0)])
segmentized = GO.segmentize(ls; min_number = 5)
@test GI.ngeom(segmentized) == 5
@test GI.getgeom(segmentized, 2) == (2.5, 2.5)
@test GI.getgeom(segmentized, 3) == (5.0, 5.0)
end

@testset_implementations "Points as vectors" begin
# If there is no line to be the second longest we risk going
# into an infinate loop.
ls = GI.LineString([[0, 0], [10, 10]])
segmentized = GO.segmentize(ls; min_number = 3)
@test GI.ngeom(segmentized) == 3
@test GI.getgeom(segmentized, 2) == [5.0, 5.0]
end

@testset_implementations "Points as Float32" begin
# If there is no line to be the second longest we risk going
# into an infinate loop.
ls = GI.LineString(Vector{Tuple{Float32, Float32}}([(0.0, 0.0), (10.0, 10.0)]))
segmentized = GO.segmentize(ls; min_number = 3)
@test GI.ngeom(segmentized) == 3
@test GI.getgeom(segmentized, 2) isa Tuple{Float32, Float32}
end
end

@testset_implementations "Wrong number of keyword arguments" begin
ls = GI.LineString([(0, 0), (1, 1), (2, 2), (3, 3)])
@test_throws ArgumentError GO.segmentize(ls)
@test_throws ArgumentError GO.segmentize(ls; max_distance = 0.5, min_number = 5)
end
Loading