diff --git a/docs/src/experiments/regridding/regridding.jl b/docs/src/experiments/regridding/regridding.jl new file mode 100644 index 000000000..2a4f0c8fc --- /dev/null +++ b/docs/src/experiments/regridding/regridding.jl @@ -0,0 +1,118 @@ +import GeoInterface as GI, GeometryOps as GO +using SortTileRecursiveTree: STRtree +using SparseArrays: spzeros +using Extents + +using CairoMakie, GeoInterfaceMakie + +include("sphericalpoints.jl") + +function area_of_intersection_operator(grid1, grid2; nodecapacity1 = 10, nodecapacity2 = 10) + area_of_intersection_operator(GO.Planar(), grid1, grid2; nodecapacity1 = nodecapacity1, nodecapacity2 = nodecapacity2) +end + +function area_of_intersection_operator(m::GO.Manifold, grid1, grid2; nodecapacity1 = 10, nodecapacity2 = 10) # grid1 and grid2 are both vectors of polygons + A = spzeros(Float64, length(grid1), length(grid2)) + # Prepare STRtrees for the two grids, to speed up intersection queries + # we may want to separately tune nodecapacity if one is much larger than the other. + # specifically we may want to tune leaf node capacity via Hilbert packing while still + # constraining inner node capacity. But that can come later. + tree1 = STRtree(grid1; nodecapacity = nodecapacity1) + tree2 = STRtree(grid2; nodecapacity = nodecapacity2) + # Do the dual query, which is the most efficient way to do this, + # by iterating down both trees simultaneously, rejecting pairs of nodes that do not intersect. + # when we find an intersection, we calculate the area of the intersection and add it to the result matrix. + GO.SpatialTreeInterface.do_dual_query(Extents.intersects, tree1, tree2) do i1, i2 + p1, p2 = grid1[i1], grid2[i2] + # may want to check if the polygons intersect first, + # to avoid antimeridian-crossing multipolygons viewing a scanline. + intersection_polys = try # can remove this now, got all the errors cleared up in the fix. + # At some future point, we may want to add the manifold here + # but for right now, GeometryOps only supports planar polygons anyway. + GO.intersection(p1, p2; target = GI.PolygonTrait()) + catch e + @error "Intersection failed!" i1 i2 + rethrow(e) + end + + area_of_intersection = GO.area(m, intersection_polys) + if area_of_intersection > 0 + A[i1, i2] += area_of_intersection + end + end + + return A +end + +grid1 = begin + gridpoints = [(i, j) for i in 0:2, j in 0:2] + [GI.Polygon([GI.LinearRing([gridpoints[i, j], gridpoints[i, j+1], gridpoints[i+1, j+1], gridpoints[i+1, j], gridpoints[i, j]])]) for i in 1:size(gridpoints, 1)-1, j in 1:size(gridpoints, 2)-1] |> vec +end + +grid2 = begin + diamondpoly = GI.Polygon([GI.LinearRing([(0, 1), (1, 2), (2, 1), (1, 0), (0, 1)])]) + trianglepolys = GI.Polygon.([ + [GI.LinearRing([(0, 0), (1, 0), (0, 1), (0, 0)])], + [GI.LinearRing([(0, 1), (0, 2), (1, 2), (0, 1)])], + [GI.LinearRing([(1, 2), (2, 1), (2, 2), (1, 2)])], + [GI.LinearRing([(2, 1), (2, 0), (1, 0), (2, 1)])], + ]) + [diamondpoly, trianglepolys...] +end + +A = area_of_intersection_operator(grid1, grid2) + +# Now, let's perform some interpolation! +area1 = vec(sum(A, dims=2)) +# test: @assert area1 == GO.area.(grid1) +area2 = vec(sum(A, dims=1)) +# test: @assert area2 == GO.area.(grid2) + +values_on_grid2 = [0, 0, 5, 0, 0] +poly(grid2; color = values_on_grid2, strokewidth = 2, strokecolor = :red) + +values_on_grid1 = A * values_on_grid2 ./ area1 +@assert sum(values_on_grid1 .* area1) == sum(values_on_grid2 .* area2) +poly(grid1; color = values_on_grid1, strokewidth = 2, strokecolor = :blue) + +values_back_on_grid2 = A' * values_on_grid1 ./ area2 +@assert sum(values_back_on_grid2 .* area2) == sum(values_on_grid2 .* area2) +poly(grid2; color = values_back_on_grid2, strokewidth = 2, strokecolor = :green) +# We can see here that some data has diffused into the central diamond cell of grid2, +# since it was overlapped by the top left cell of grid1. + + +using SpeedyWeather +using GeoMakie + +SpeedyWeatherGeoMakieExt = Base.get_extension(SpeedyWeather, :SpeedyWeatherGeoMakieExt) + +grid1 = rand(OctaHEALPixGrid, 5 + 100) +grid2 = rand(FullGaussianGrid, 4 + 100) + +faces1 = SpeedyWeatherGeoMakieExt.get_faces(grid1) +faces2 = SpeedyWeatherGeoMakieExt.get_faces(grid2) + +polys1 = GI.Polygon.(GI.LinearRing.(eachcol(faces1))) .|> GO.CutAtAntimeridianAndPoles() .|> GO.fix +polys2 = GI.Polygon.(GI.LinearRing.(eachcol(faces2))) .|> GO.CutAtAntimeridianAndPoles() .|> GO.fix + +A = @time area_of_intersection_operator(polys1, polys2) + +p1 = polys1[93] +p2 = polys2[105] + +f, a, p = poly(p1) +poly!(a, p2) +f + +# bug found in Foster Hormann tracing +# but geos also does the same thing +boxpoly = GI.Polygon([GI.LinearRing([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)])]) +diamondpoly = GI.Polygon([GI.LinearRing([(0, 1), (1, 2), (2, 1), (1, 0), (0, 1)])]) + +diffpoly = GO.difference(boxpoly, diamondpoly; target = GI.PolygonTrait()) |> only +cutpolys = GO.cut(diffpoly, GI.Line([(0, 0), (4, 0)])) # even cut misbehaves! + + + + diff --git a/docs/src/experiments/regridding/sphericalpoints.jl b/docs/src/experiments/regridding/sphericalpoints.jl new file mode 100644 index 000000000..4444a1ca9 --- /dev/null +++ b/docs/src/experiments/regridding/sphericalpoints.jl @@ -0,0 +1,116 @@ +import GeoInterface as GI +import GeometryOps as GO +import LinearAlgebra +import LinearAlgebra: dot, cross + +## Get the area of a LinearRing with coordinates in radians +struct SphericalPoint{T <: Real} + data::NTuple{3, T} +end +SphericalPoint(x, y, z) = SphericalPoint((x, y, z)) + +# define the 4 basic mathematical operators elementwise on the data tuple +Base.:+(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .+ q.data) +Base.:-(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .- q.data) +Base.:*(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .* q.data) +Base.:/(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data ./ q.data) +# Define sum on a SphericalPoint to sum across its data +Base.sum(p::SphericalPoint) = sum(p.data) + +# define dot and cross products +LinearAlgebra.dot(p::SphericalPoint, q::SphericalPoint) = sum(p * q) +function LinearAlgebra.cross(a::SphericalPoint, b::SphericalPoint) + a1, a2, a3 = a.data + b1, b2, b3 = b.data + SphericalPoint((a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1)) +end + +# Using Eriksson's formula for the area of spherical triangles: https://www.jstor.org/stable/2691141 +# This melts down only when two points are antipodal, which in our case will not happen. +function _unit_spherical_triangle_area(a, b, c) + #t = abs(dot(a, cross(b, c))) + #t /= 1 + dot(b,c) + dot(c, a) + dot(a, b) + t = abs(dot(a, (cross(b - a, c - a))) / dot(b + a, c + a)) + 2*atan(t) +end + +_lonlat_to_sphericalpoint(p) = _lonlat_to_sphericalpoint(GI.x(p), GI.y(p)) +function _lonlat_to_sphericalpoint(lon, lat) + lonsin, loncos = sincosd(lon) + latsin, latcos = sincosd(lat) + x = latcos * loncos + y = latcos * lonsin + z = latsin + return SphericalPoint(x,y,z) +end + + + +# Extend area to spherical + +# TODO: make this the other way around, but that can wait. +GO.area(m::GO.Planar, geoms) = GO.area(geoms) + +function GO.area(m::GO.Spherical, geoms) + return GO.applyreduce(+, GI.PolygonTrait(), geoms) do poly + GO.area(m, GI.PolygonTrait(), poly) + end * ((-m.radius^2)/ 2) # do this after the sum, to increase accuracy and minimize calculations. +end + +function GO.area(m::GO.Spherical, ::GI.PolygonTrait, poly) + area = abs(_ring_area(m, GI.getexterior(poly))) + for interior in GI.gethole(poly) + area -= abs(_ring_area(m, interior)) + end + return area +end + +function _ring_area(m::GO.Spherical, ring) + # convert ring to a sequence of SphericalPoints + points = _lonlat_to_sphericalpoint.(GI.getpoint(ring))[1:end-1] # deliberately drop the closing point + p1, p2, p3 = points[1], points[2], points[3] + # For a spherical polygon, we can compute the area by splitting it into triangles + # and summing their areas. We use the first point as a common vertex for all triangles. + area = 0.0 + # Sum areas of triangles formed by first point and consecutive pairs of points + np = length(points) + for i in 1:np + p1, p2, p3 = p2, p3, points[i] + area += _unit_spherical_triangle_area(p1, p2, p3) + end + return area +end + +function _ring_area(m::GO.Spherical, ring) + # Convert ring points to spherical coordinates + points = GO.tuples(ring).geom + + # Remove last point if it's the same as first (closed ring) + if points[end] == points[1] + points = points[1:end-1] + end + + n = length(points) + if n < 3 + return 0.0 + end + + area = 0.0 + + # Use L'Huilier's formula to sum the areas of spherical triangles + # formed by first point and consecutive pairs of points + for i in 1:n + p1, p2, p3 = points[mod1(i-1, n)], points[mod1(i, n)], points[mod1(i+1, n)] + area += sind(GI.y(p2)) * (GI.x(p3) - GI.x(p1)) + end + + return area +end + + + + +# Test the area calculation +p1 = GI.Polygon([GI.LinearRing(Point2f[(0, 0), (1, 0), (0, 1), (0, 0)] .- (Point2f(0.5, 0.5),))]) + +GO.area(GO.Spherical(), p1) \ No newline at end of file diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 23972c504..b2faaf35e 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -93,6 +93,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/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 new file mode 100644 index 000000000..6a88040ac --- /dev/null +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -0,0 +1,358 @@ +#= +# 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 + "The value at the right edge of the antimeridian, in your angular units" + right::Float64 = 180.0 + "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 + "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::CutAtAntimeridianAndPoles) + print(io, "CutAtAntimeridianAndPoles(left=$(cut.left), right=$(cut.right))") +end +Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridianAndPoles) = Base.show(io, cut) + +application_level(::CutAtAntimeridianAndPoles) = TraitTarget(GI.PolygonTrait(), GI.LineStringTrait(), GI.MultiLineStringTrait(), GI.MultiPolygonTrait()) + +function (c::CutAtAntimeridianAndPoles)(trait::Union{GI.PolygonTrait, GI.MultiPolygonTrait, GI.LineStringTrait, GI.MultiLineStringTrait}, geom) + return _AntimeridianHelpers.cut_at_antimeridian(trait, geom) +end + +module _AntimeridianHelpers + +import GeoInterface as GI +import ..GeometryOps as GO +using ..GeometryOps: CutAtAntimeridianAndPoles + +# 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(c::CutAtAntimeridianAndPoles, point::Tuple{Float64,Float64})::Tuple{Float64,Float64,Float64} + # TODO: handle non-degree domains somehow + 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(c::CutAtAntimeridianAndPoles, start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64})::Float64 + # Convert points to 3D vectors + p1 = spherical_degrees_to_cartesian(c, start) + p2 = spherical_degrees_to_cartesian(c, endpoint) + + # Cross product defines plane through both points + n1 = _cross(p1, p2) + + # Unit vector that defines the meridian plane + n2 = spherical_degrees_to_cartesian(c, (c.left, 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(c::CutAtAntimeridianAndPoles, start::Tuple{Float64,Float64}, stop::Tuple{Float64,Float64})::Float64 + lat_delta = stop[2] - start[2] + if stop[1] > 0 + round( + start[2] + (c.right - start[1]) * lat_delta / (stop[1] + c.period - start[1]), + digits=7 + ) + else + round( + start[2] + (start[1] - c.left) * lat_delta / (start[1] + c.period - stop[1]), + digits=7 + ) + end +end + +# Main crossing latitude calculation function +function crossing_latitude(c::CutAtAntimeridianAndPoles, 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(c, start, endpoint) : crossing_latitude_flat(c, start, endpoint) +end + +# Normalize coordinates to ensure longitudes are between -180 and 180 +function normalize_coords(c::CutAtAntimeridianAndPoles, 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], c.right) + if abs(point[2]) != c.northpole && isapprox(prev_point[1], c.left) + normalized[i] = (c.left, point[2]) + else + normalized[i] = (c.right, point[2]) + end + elseif isapprox(point[1], c.left) + if abs(point[2]) != c.northpole && isapprox(prev_point[1], c.right) + normalized[i] = (c.right, point[2]) + else + normalized[i] = (c.left, point[2]) + end + else + normalized[i] = (((point[1] - c.left) % c.period) + c.left, 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(c::CutAtAntimeridianAndPoles, 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(c, 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(c, 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 segments +end + +# Check if a segment is self-closing +function is_self_closing(c::CutAtAntimeridianAndPoles, segment::Vector{Tuple{Float64,Float64}})::Bool + is_right = segment[end][1] == c.right + 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(c::CutAtAntimeridianAndPoles, segments::Vector{Vector{Tuple{Float64,Float64}}})::Vector{GI.Polygon} + isempty(segments) && return GI.Polygon[] + + segment = pop!(segments) + is_right = segment[end][1] == c.right + candidates = Tuple{Union{Nothing,Int},Float64}[] + + if is_self_closing(c, 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(c, s) || s[end][2] < segment[1][2])) || + (!is_right && s[1][2] < segment[end][2] && + (!is_self_closing(c, 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(c, segments) + else + # Handle self-joining segment + polygons = build_polygons(c, segments) + if !all(p == segment[1] for p in segment) + push!(polygons, GI.Polygon([segment])) + end + return filter(polygons) do p + x1 = GI.x(first(GI.getpoint(p))) + !all(p -> GI.x(p) == x1, GI.getpoint(p)) + end + end +end + +# Main function to cut a polygon at the antimeridian +cut_at_antimeridian(x) = cut_at_antimeridian(CutAtAntimeridianAndPoles(), GI.trait(x), x) +cut_at_antimeridian(c::CutAtAntimeridianAndPoles, x) = cut_at_antimeridian(c, GI.trait(x), x) +cut_at_antimeridian(t::GI.AbstractTrait, x) = cut_at_antimeridian(CutAtAntimeridianAndPoles(), t, x) + +function cut_at_antimeridian( + c::CutAtAntimeridianAndPoles, + ::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(c, exterior) + + # Segment the exterior ring + segments = segment_ring(c, 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.((c,), 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(c, 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(c, reverse(hole), great_circle) + end + end + append!(segments, hole_segments) + else + push!(holes, hole) + end + end + + # Build final polygons + result_polygons = build_polygons(c, 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 + + filter!(poly -> !iszero(GO.area(poly)), result_polygons) + + return length(result_polygons) == 1 ? result_polygons[1] : GI.MultiPolygon(result_polygons) +end + +function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.AbstractCurveTrait, line::T; great_circle::Bool=true) where {T} + coords = GO.tuples(line).geom + segments = segment_ring(c, coords, great_circle) + + if isempty(segments) + return line + else + return GI.MultiLineString(segments) + end +end + + +function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.MultiPolygonTrait, x; kwargs...) + results = GI.Polygon[] + for poly in GI.getgeom(x) + result = cut_at_antimeridian(c, 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(c::CutAtAntimeridianAndPoles, t::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(c, 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 diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index 5cccb5ea2..a9e227f13 100644 --- a/src/transformations/correction/geometry_correction.jl +++ b/src/transformations/correction/geometry_correction.jl @@ -33,28 +33,49 @@ 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 +# 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) (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()], 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) - 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 + =# + return final_geoms end # ## Available corrections 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)