From f7170d0224c7bf106c07538e2d31f49c1a6c06c9 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 16:51:22 -0400 Subject: [PATCH 1/7] implement an antimeridian cutting fix ported from the python package --- src/GeometryOps.jl | 1 + .../correction/cut_at_antimeridian.jl | 348 ++++++++++++++++++ 2 files changed, 349 insertions(+) create mode 100644 src/transformations/correction/cut_at_antimeridian.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 152899db9..53f4deb70 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -75,6 +75,7 @@ include("transformations/forcedims.jl") include("transformations/correction/geometry_correction.jl") include("transformations/correction/closed_ring.jl") include("transformations/correction/intersecting_polygons.jl") +include("transformations/correction/cut_at_antimeridian.jl") # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`. for name in names(GeoInterface) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl new file mode 100644 index 000000000..ca60a3d7f --- /dev/null +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -0,0 +1,348 @@ +#= +# Antimeridian Cutting +=# +export CutAtAntimeridianAndPoles # TODO: too wordy? + +#= +This correction cuts the geometry at the antimeridian and the poles, and returns a fixed geometry. + +The implementation is translated from the implementation in https://github.com/gadomski/antimeridian, +which is a Python package. Several ports of that algorithm exist in e.g. Go, Rust, etc. + +At some point we will have to go in and clean up the implementation, remove all the hardcoded code, +and make it more efficient by using raw geointerface, and not allocating so much (perhaps, by passing allocs around). + +But for right now it works, that's all I need. +=# + +""" + CutAtAntimeridianAndPoles(; kwargs...) <: GeometryCorrection + +This correction cuts the geometry at the antimeridian and the poles, and returns a fixed geometry. +""" +Base.@kwdef struct CutAtAntimeridianAndPoles <: GeometryCorrection + "The value at the north pole, in your angular units" + northpole::Float64 = 90.0 # TODO not used!!! + "The value at the south pole, in your angular units" + southpole::Float64 = -90.0 # TODO not used!!! + "The value at the left edge of the antimeridian, in your angular units" + left::Float64 = -180.0 # TODO not used!!! + "The value at the right edge of the antimeridian, in your angular units" + right::Float64 = 180.0 # TODO not used!!! + "The period of the cyclic / cylindrical coordinate system's x value, usually computed automatically so you don't have to provide it." + period::Float64 = right - left # TODO not used!!! + "If the polygon is known to enclose the north pole, set this to true" + force_north_pole::Bool=false # TODO not used!!! + "If the polygon is known to enclose the south pole, set this to true" + force_south_pole::Bool=false # TODO not used!!! + "If true, use the great circle method to find the antimeridian crossing, otherwise use the flat projection method. There is no reason to set this to be off." + great_circle = true +end + +function Base.show(io::IO, cut::CutAtAntimeridian) + print(io, "CutAtAntimeridian(left=$(cut.left), right=$(cut.right))") +end +Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridian) = Base.show(io, cut) + +application_level(::CutAtAntimeridian) = Union{GI.PolygonTrait, GI.LineStringTrait, GI.MultiLineStringTrait, GI.MultiPolygonTrait} + +function (c::CutAtAntimeridianAndPoles)(trait::GI.AbstractTrait, geom) + return _AntimeridianHelpers.cut_at_antimeridian(trait, geom) +end + +module _AntimeridianHelpers + +import GeoInterface as GI +import ..GeometryOps as GO + +# Custom cross product for 3D tuples +function _cross(a::Tuple{Float64,Float64,Float64}, b::Tuple{Float64,Float64,Float64}) + return ( + a[2]*b[3] - a[3]*b[2], + a[3]*b[1] - a[1]*b[3], + a[1]*b[2] - a[2]*b[1] + ) +end + +# Convert spherical degrees to cartesian coordinates +function spherical_degrees_to_cartesian(point::Tuple{Float64,Float64})::Tuple{Float64,Float64,Float64} + lon, lat = point + slon, clon = sincosd(lon) + slat, clat = sincosd(lat) + return ( + clon * clat, + slon * clat, + slat + ) +end + +# Calculate crossing latitude using great circle method +function crossing_latitude_great_circle(start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64})::Float64 + # Convert points to 3D vectors + p1 = spherical_degrees_to_cartesian(start) + p2 = spherical_degrees_to_cartesian(endpoint) + + # Cross product defines plane through both points + n1 = _cross(p1, p2) + + # Unit vector -Y defines meridian plane + n2 = (0.0, -1.0, 0.0) + + # Intersection of planes defined by cross product + intersection = _cross(n1, n2) + norm = sqrt(sum(intersection .^ 2)) + intersection = intersection ./ norm + + # Convert back to spherical coordinates (just need latitude) + round(asind(intersection[3]), digits=7) +end + +# Calculate crossing latitude using flat projection method +function crossing_latitude_flat(start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64})::Float64 + lat_delta = endpoint[2] - start[2] + if endpoint[1] > 0 + round( + start[2] + (180.0 - start[1]) * lat_delta / (endpoint[1] + 360.0 - start[1]), + digits=7 + ) + else + round( + start[2] + (start[1] + 180.0) * lat_delta / (start[1] + 360.0 - endpoint[1]), + digits=7 + ) + end +end + +# Main crossing latitude calculation function +function crossing_latitude(start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64}, great_circle::Bool)::Float64 + abs(start[1]) == 180 && return start[2] + abs(endpoint[1]) == 180 && return endpoint[2] + + return great_circle ? crossing_latitude_great_circle(start, endpoint) : crossing_latitude_flat(start, endpoint) +end + +# Normalize coordinates to ensure longitudes are between -180 and 180 +function normalize_coords(coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{Float64,Float64}} + normalized = deepcopy(coords) + all_on_antimeridian = true + + for i in eachindex(normalized) + point = normalized[i] + prev_point = normalized[mod1(i-1, length(normalized))] + + if isapprox(point[1], 180) + if abs(point[2]) != 90 && isapprox(prev_point[1], -180) + normalized[i] = (-180.0, point[2]) + else + normalized[i] = (180.0, point[2]) + end + elseif isapprox(point[1], -180) + if abs(point[2]) != 90 && isapprox(prev_point[1], 180) + normalized[i] = (180.0, point[2]) + else + normalized[i] = (-180.0, point[2]) + end + else + normalized[i] = (((point[1] + 180) % 360) - 180, point[2]) + all_on_antimeridian = false + end + end + + return all_on_antimeridian ? coords : normalized +end + +# Segment a ring of coordinates at antimeridian crossings +function segment_ring(coords::Vector{Tuple{Float64,Float64}}, great_circle::Bool)::Vector{Vector{Tuple{Float64,Float64}}} + segments = Vector{Vector{Tuple{Float64,Float64}}}() + current_segment = Tuple{Float64,Float64}[] + + for i in 1:length(coords)-1 + start = coords[i] + endpoint = coords[i+1] + push!(current_segment, start) + + # Check for antimeridian crossing + if (endpoint[1] - start[1] > 180) && (endpoint[1] - start[1] != 360) # left crossing + lat = crossing_latitude(start, endpoint, great_circle) + push!(current_segment, (-180.0, lat)) + push!(segments, current_segment) + current_segment = [(180.0, lat)] + elseif (start[1] - endpoint[1] > 180) && (start[1] - endpoint[1] != 360) # right crossing + lat = crossing_latitude(endpoint, start, great_circle) + push!(current_segment, (180.0, lat)) + push!(segments, current_segment) + current_segment = [(-180.0, lat)] + end + end + + # Handle last point and segment + if isempty(segments) + return Vector{Vector{Tuple{Float64,Float64}}}() # No crossings + elseif coords[end] == segments[1][1] + # Join polygons + segments[1] = vcat(current_segment, segments[1]) + else + push!(current_segment, coords[end]) + push!(segments, current_segment) + end + + return normalize_coords.(segments) +end + +# Check if a segment is self-closing +function is_self_closing(segment::Vector{Tuple{Float64,Float64}})::Bool + is_right = segment[end][1] == 180 + return segment[1][1] == segment[end][1] && ( + (is_right && segment[1][2] > segment[end][2]) || + (!is_right && segment[1][2] < segment[end][2]) + ) +end + +# Build polygons from segments +function build_polygons(segments::Vector{Vector{Tuple{Float64,Float64}}})::Vector{GI.Polygon} + isempty(segments) && return GI.Polygon[] + + segment = pop!(segments) + is_right = segment[end][1] == 180 + candidates = Tuple{Union{Nothing,Int},Float64}[] + + if is_self_closing(segment) + push!(candidates, (nothing, segment[1][2])) + end + + for (i, s) in enumerate(segments) + if s[1][1] == segment[end][1] + if (is_right && s[1][2] > segment[end][2] && + (!is_self_closing(s) || s[end][2] < segment[1][2])) || + (!is_right && s[1][2] < segment[end][2] && + (!is_self_closing(s) || s[end][2] > segment[1][2])) + push!(candidates, (i, s[1][2])) + end + end + end + + # Sort candidates + sort!(candidates, by=x->x[2], rev=!is_right) + + index = isempty(candidates) ? nothing : candidates[1][1] + + if !isnothing(index) + # Join segments and recurse + segment = vcat(segment, segments[index]) + segments[index] = segment + return build_polygons(segments) + else + # Handle self-joining segment + polygons = build_polygons(segments) + if !all(p == segment[1] for p in segment) + push!(polygons, GI.Polygon([segment])) + end + return polygons + end +end + +# Main function to cut a polygon at the antimeridian +cut_at_antimeridian(x) = cut_at_antimeridian(GI.trait(x), x) + +function cut_at_antimeridian( + ::GI.PolygonTrait, + polygon::T; + force_north_pole::Bool=false, + force_south_pole::Bool=false, + fix_winding::Bool=true, + great_circle::Bool=true +) where {T} + # Get exterior ring + exterior = GO.tuples(GI.getexterior(polygon)).geom + exterior = normalize_coords(exterior) + + # Segment the exterior ring + segments = segment_ring(exterior, great_circle) + + if isempty(segments) + # No antimeridian crossing + if fix_winding && GO.isclockwise(GI.LinearRing(exterior)) + coord_vecs = reverse.(getproperty.(GO.tuples.(GI.getring(polygon)), :geom)) + return GI.Polygon(normalize_coords.(coord_vecs)) + end + return polygon + end + + # Handle holes + holes = Vector{Vector{Tuple{Float64,Float64}}}() + for hole_idx in 1:GI.nhole(polygon) + hole = GO.tuples(GI.gethole(polygon, hole_idx)).geom + hole_segments = segment_ring(hole, great_circle) + + if !isempty(hole_segments) + if fix_winding + unwrapped = [(x % 360, y) for (x, y) in hole] + if !GO.isclockwise(GI.LineString(unwrapped)) + hole_segments = segment_ring(reverse(hole), great_circle) + end + end + append!(segments, hole_segments) + else + push!(holes, hole) + end + end + + # Build final polygons + result_polygons = build_polygons(segments) + + # Add holes to appropriate polygons + for hole in holes + for (i, poly) in enumerate(result_polygons) + if GO.contains(poly, GI.Point(hole[1])) + rings = GI.getring(poly) + push!(rings, hole) + result_polygons[i] = GI.Polygon(rings) + break + end + end + end + + return length(result_polygons) == 1 ? result_polygons[1] : GI.MultiPolygon(result_polygons) +end + +function cut_at_antimeridian(::GI.AbstractCurveTrait, line::T; great_circle::Bool=true) where {T} + coords = GO.tuples(line).geom + segments = segment_ring(coords, great_circle) + + if isempty(segments) + return line + else + return GI.MultiLineString(segments) + end +end + + +function cut_at_antimeridian(::GI.MultiPolygonTrait, x; kwargs...) + results = GI.Polygon[] + for poly in GI.getgeom(x) + result = cut_at_antimeridian(GI.PolygonTrait(), poly; kwargs...) + if result isa GI.Polygon + push!(results, result) + elseif result isa GI.MultiPolygon + append!(results, result.geom) + end + end + return GI.MultiPolygon(results) +end + +function cut_at_antimeridian(::GI.MultiLineStringTrait, multiline::T; great_circle::Bool=true) where {T} + linestrings = Vector{Vector{Tuple{Float64,Float64}}}() + + for line in GI.getgeom(multiline) + fixed = cut_at_antimeridian(GI.LineStringTrait(), line; great_circle) + if fixed isa GI.LineString + push!(linestrings, GO.tuples(fixed).geom) + else + append!(linestrings, GO.tuples.(GI.getgeom(fixed)).geom) + end + end + + return GI.MultiLineString(linestrings) +end + +end \ No newline at end of file From 129564343db2a1ccf3d1f42db563c3f4a72d67fc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 16:51:45 -0400 Subject: [PATCH 2/7] use traittarget in application_level for fix --- src/transformations/correction/closed_ring.jl | 2 +- .../correction/cut_at_antimeridian.jl | 2 +- .../correction/geometry_correction.jl | 23 ++++++++++++------- .../correction/intersecting_polygons.jl | 4 ++-- 4 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/transformations/correction/closed_ring.jl b/src/transformations/correction/closed_ring.jl index b8076b821..7f4c7e905 100644 --- a/src/transformations/correction/closed_ring.jl +++ b/src/transformations/correction/closed_ring.jl @@ -46,7 +46,7 @@ See also [`GeometryCorrection`](@ref). """ struct ClosedRing <: GeometryCorrection end -application_level(::ClosedRing) = GI.PolygonTrait +application_level(::ClosedRing) = TraitTarget(GI.PolygonTrait()) function (::ClosedRing)(::GI.PolygonTrait, polygon) exterior = _close_linear_ring(GI.getexterior(polygon)) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index ca60a3d7f..7ea1e1e28 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -44,7 +44,7 @@ function Base.show(io::IO, cut::CutAtAntimeridian) end Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridian) = Base.show(io, cut) -application_level(::CutAtAntimeridian) = Union{GI.PolygonTrait, GI.LineStringTrait, GI.MultiLineStringTrait, GI.MultiPolygonTrait} +application_level(::CutAtAntimeridian) = TraitTarget(GI.PolygonTrait(), GI.LineStringTrait(), GI.MultiLineStringTrait(), GI.MultiPolygonTrait()) function (c::CutAtAntimeridianAndPoles)(trait::GI.AbstractTrait, geom) return _AntimeridianHelpers.cut_at_antimeridian(trait, geom) diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index 5cccb5ea2..c252eb822 100644 --- a/src/transformations/correction/geometry_correction.jl +++ b/src/transformations/correction/geometry_correction.jl @@ -33,8 +33,13 @@ This abstract type represents a geometry correction. ## Interface Any `GeometryCorrection` must implement two functions: - * `application_level(::GeometryCorrection)::AbstractGeometryTrait`: This function should return the `GeoInterface` trait that the correction is intended to be applied to, like `PointTrait` or `LineStringTrait` or `PolygonTrait`. - * `(::GeometryCorrection)(::AbstractGeometryTrait, geometry)::(some_geometry)`: This function should apply the correction to the given geometry, and return a new geometry. + - `application_level(::GeometryCorrection)::TraitTarget`: This function should + return the `GeoInterface` trait that the correction is intended to be applied to, + like `PointTrait` or `LineStringTrait` or `PolygonTrait`. It can also return a + union of traits via `TraitTarget`, but that behaviour is a bit tricky... + - `(::GeometryCorrection)(::AbstractGeometryTrait, geometry)::(some_geometry)`: + This function should apply the correction to the given geometry, and return a new + geometry. """ abstract type GeometryCorrection end @@ -44,15 +49,17 @@ application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc) (gc::GeometryCorrection)(trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).") -function fix(geometry; corrections = GeometryCorrection[ClosedRing(),], kwargs...) - traits = application_level.(corrections) +function fix(geometry; corrections = GeometryCorrection[ClosedRing(), CutAtAntimeridianAndPoles()], kwargs...) + application_levels = application_level.(corrections) final_geometry = geometry - for Trait in (GI.PointTrait, GI.MultiPointTrait, GI.LineStringTrait, GI.LinearRingTrait, GI.MultiLineStringTrait, GI.PolygonTrait, GI.MultiPolygonTrait) - available_corrections = findall(x -> x == Trait, traits) + for trait in (GI.PointTrait(), GI.MultiPointTrait(), GI.LineStringTrait(), GI.LinearRingTrait(), GI.MultiLineStringTrait(), GI.PolygonTrait(), GI.MultiPolygonTrait()) + available_corrections = findall(x -> trait in x, application_levels) isempty(available_corrections) && continue - @debug "Correcting for $(Trait)" + @debug "Correcting for $(trait), with corrections: " available_corrections net_function = reduce(∘, corrections[available_corrections]) - final_geometry = apply(net_function, Trait, final_geometry; kwargs...) + # TODO: this allocates too much, because it keeps reconstructing higher level geoms. + # We might want some way to embed the fixes in reconstruct/rebuild, which would imply a modified apply pipeline... + final_geometry = apply(net_function, trait, final_geometry; kwargs...) end return final_geometry end diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl index 2223f6b89..73d857525 100644 --- a/src/transformations/correction/intersecting_polygons.jl +++ b/src/transformations/correction/intersecting_polygons.jl @@ -56,7 +56,7 @@ See also [`GeometryCorrection`](@ref). """ struct UnionIntersectingPolygons <: GeometryCorrection end -application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait +application_level(::UnionIntersectingPolygons) = TraitTarget(GI.MultiPolygonTrait()) function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly) union_multipoly = tuples(multipoly) @@ -99,7 +99,7 @@ See also [`GeometryCorrection`](@ref), [`UnionIntersectingPolygons`](@ref). """ struct DiffIntersectingPolygons <: GeometryCorrection end -application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait +application_level(::DiffIntersectingPolygons) = TraitTarget(GI.MultiPolygonTrait()) function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly) diff_multipoly = tuples(multipoly) From e8f8f9f196c19c2e103df207d49faa1b61e82ba1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 17:06:33 -0400 Subject: [PATCH 3/7] fix typos --- src/transformations/correction/cut_at_antimeridian.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index 7ea1e1e28..1e8a92a52 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -39,12 +39,12 @@ Base.@kwdef struct CutAtAntimeridianAndPoles <: GeometryCorrection great_circle = true end -function Base.show(io::IO, cut::CutAtAntimeridian) - print(io, "CutAtAntimeridian(left=$(cut.left), right=$(cut.right))") +function Base.show(io::IO, cut::CutAtAntimeridianAndPoles) + print(io, "CutAtAntimeridianAndPoles(left=$(cut.left), right=$(cut.right))") end -Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridian) = Base.show(io, cut) +Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridianAndPoles) = Base.show(io, cut) -application_level(::CutAtAntimeridian) = TraitTarget(GI.PolygonTrait(), GI.LineStringTrait(), GI.MultiLineStringTrait(), GI.MultiPolygonTrait()) +application_level(::CutAtAntimeridianAndPoles) = TraitTarget(GI.PolygonTrait(), GI.LineStringTrait(), GI.MultiLineStringTrait(), GI.MultiPolygonTrait()) function (c::CutAtAntimeridianAndPoles)(trait::GI.AbstractTrait, geom) return _AntimeridianHelpers.cut_at_antimeridian(trait, geom) From 183fc30220be47b1887fe39b3acc2518a4d1dff5 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 18:09:31 -0400 Subject: [PATCH 4/7] make geomcorrections broadcast as scalars enables geoms .|> fix since fix doesnt have apply integrated(!) we should fix --- src/transformations/correction/geometry_correction.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index c252eb822..abffbb862 100644 --- a/src/transformations/correction/geometry_correction.jl +++ b/src/transformations/correction/geometry_correction.jl @@ -43,6 +43,9 @@ Any `GeometryCorrection` must implement two functions: """ abstract type GeometryCorrection end +# Make sure that geometry corrections are treated as scalars when broadcasting. +Base.Broadcast.broadcastable(c::GeometryCorrection) = (c,) + application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc)") (gc::GeometryCorrection)(geometry) = gc(GI.trait(geometry), geometry) From 4e87297684122e2f4c8d8f71fa3049942ce8fcac Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 18:10:44 -0400 Subject: [PATCH 5/7] apply fixes iteratively, not by trait this is worse for performance, better for predictability and correctness (the most complex trait will be chosen for the application rather than the least complex) --- .../correction/geometry_correction.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index abffbb862..eb1ce5737 100644 --- a/src/transformations/correction/geometry_correction.jl +++ b/src/transformations/correction/geometry_correction.jl @@ -52,7 +52,16 @@ application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc) (gc::GeometryCorrection)(trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).") -function fix(geometry; corrections = GeometryCorrection[ClosedRing(), CutAtAntimeridianAndPoles()], kwargs...) +function fix(geometry; corrections = GeometryCorrection[CutAtAntimeridianAndPoles(), ClosedRing()], kwargs...) + final_geoms = geometry + # Iterate through the corrections and apply them to the input. + # This allocates a _lot_, especially when reconstructing tables, + # but it's the only fully general way to do this that I can think of. + for correction in corrections + final_geoms = apply(correction, application_level(correction), final_geoms; kwargs...) + end + #= + # This was the old implementation application_levels = application_level.(corrections) final_geometry = geometry for trait in (GI.PointTrait(), GI.MultiPointTrait(), GI.LineStringTrait(), GI.LinearRingTrait(), GI.MultiLineStringTrait(), GI.PolygonTrait(), GI.MultiPolygonTrait()) @@ -65,6 +74,8 @@ function fix(geometry; corrections = GeometryCorrection[ClosedRing(), CutAtAntim final_geometry = apply(net_function, trait, final_geometry; kwargs...) end return final_geometry + =# + return final_geoms end # ## Available corrections From 917194a81c4bdf7051508101c9cdb528f10dc0d9 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 18:11:02 -0400 Subject: [PATCH 6/7] Filter out zero area polygonal shards --- src/transformations/correction/cut_at_antimeridian.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index 1e8a92a52..bcac8a0d3 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -301,6 +301,8 @@ function cut_at_antimeridian( end end end + + filter!(poly -> !iszero(GO.area(poly)), result_polygons) return length(result_polygons) == 1 ? result_polygons[1] : GI.MultiPolygon(result_polygons) end From 4d5bfc65f7c0d4284f43fc57fc7a0a63c90b42e9 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 18:11:16 -0400 Subject: [PATCH 7/7] don't unnecessarily normalize - that wasn't the place that needed it --- src/transformations/correction/cut_at_antimeridian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index bcac8a0d3..a5ad34d9a 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -186,7 +186,7 @@ function segment_ring(coords::Vector{Tuple{Float64,Float64}}, great_circle::Bool push!(segments, current_segment) end - return normalize_coords.(segments) + return segments end # Check if a segment is self-closing