From e70d20ea40adc67df1d9cf132d721f70b1e6c23d Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 18 Mar 2025 16:51:22 -0400 Subject: [PATCH 1/4] implement an antimeridian cutting fix ported from the python package --- docs/src/experiments/regridding/regridding.jl | 118 ++++++ .../experiments/regridding/sphericalpoints.jl | 116 ++++++ src/GeometryOps.jl | 1 + src/transformations/correction/closed_ring.jl | 2 +- .../correction/cut_at_antimeridian.jl | 353 ++++++++++++++++++ .../correction/geometry_correction.jl | 37 +- .../correction/intersecting_polygons.jl | 4 +- 7 files changed, 620 insertions(+), 11 deletions(-) create mode 100644 docs/src/experiments/regridding/regridding.jl create mode 100644 docs/src/experiments/regridding/sphericalpoints.jl create mode 100644 src/transformations/correction/cut_at_antimeridian.jl diff --git a/docs/src/experiments/regridding/regridding.jl b/docs/src/experiments/regridding/regridding.jl new file mode 100644 index 000000000..6833d387e --- /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.fix +polys2 = GI.Polygon.(GI.LinearRing.(eachcol(faces2))) .|> 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..43bd7be9b --- /dev/null +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -0,0 +1,353 @@ +#= +# 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::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 + +# 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 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 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(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 + + 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(::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 diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index 5cccb5ea2..eb1ce5737 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[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) - 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) From fa14077962f120fe21d4ba6327d6f5bf4266f262 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 4 Apr 2025 08:54:06 -0400 Subject: [PATCH 2/4] antimeridian cutting on any domain --- .../correction/cut_at_antimeridian.jl | 97 ++++++++++--------- 1 file changed, 50 insertions(+), 47 deletions(-) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index 43bd7be9b..2cbc4239f 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -26,11 +26,11 @@ Base.@kwdef struct CutAtAntimeridianAndPoles <: GeometryCorrection "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!!! + left::Float64 = -180.0 "The value at the right edge of the antimeridian, in your angular units" - right::Float64 = 180.0 # TODO not used!!! + 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 # TODO not used!!! + 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" @@ -77,7 +77,7 @@ function spherical_degrees_to_cartesian(point::Tuple{Float64,Float64})::Tuple{Fl end # Calculate crossing latitude using great circle method -function crossing_latitude_great_circle(start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64})::Float64 +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(start) p2 = spherical_degrees_to_cartesian(endpoint) @@ -85,8 +85,8 @@ function crossing_latitude_great_circle(start::Tuple{Float64,Float64}, 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) + # Unit vector that defines the meridian plane + n2 = spherical_degrees_to_cartesian(c.left, 0.0) # Intersection of planes defined by cross product intersection = _cross(n1, n2) @@ -98,31 +98,31 @@ function crossing_latitude_great_circle(start::Tuple{Float64,Float64}, endpoint: 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 +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] + (180.0 - start[1]) * lat_delta / (endpoint[1] + 360.0 - start[1]), + start[2] + (c.right - start[1]) * lat_delta / (stop[1] + c.period - start[1]), digits=7 ) else round( - start[2] + (start[1] + 180.0) * lat_delta / (start[1] + 360.0 - endpoint[1]), + 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(start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64}, great_circle::Bool)::Float64 +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(start, endpoint) : crossing_latitude_flat(start, endpoint) + 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(coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{Float64,Float64}} +function normalize_coords(c::CutAtAntimeridianAndPoles, coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{Float64,Float64}} normalized = deepcopy(coords) all_on_antimeridian = true @@ -130,20 +130,20 @@ function normalize_coords(coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{ 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]) + 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] = (180.0, point[2]) + normalized[i] = (c.right, point[2]) end - elseif isapprox(point[1], -180) - if abs(point[2]) != 90 && isapprox(prev_point[1], 180) - normalized[i] = (180.0, point[2]) + 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] = (-180.0, point[2]) + normalized[i] = (c.left, point[2]) end else - normalized[i] = (((point[1] + 180) % 360) - 180, point[2]) + normalized[i] = (((point[1] - c.left) % c.period) + c.left, point[2]) all_on_antimeridian = false end end @@ -152,7 +152,7 @@ function normalize_coords(coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{ 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}}} +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}[] @@ -163,12 +163,12 @@ function segment_ring(coords::Vector{Tuple{Float64,Float64}}, great_circle::Bool # Check for antimeridian crossing if (endpoint[1] - start[1] > 180) && (endpoint[1] - start[1] != 360) # left crossing - lat = crossing_latitude(start, endpoint, great_circle) + 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(endpoint, start, great_circle) + lat = crossing_latitude(c, endpoint, start, great_circle) push!(current_segment, (180.0, lat)) push!(segments, current_segment) current_segment = [(-180.0, lat)] @@ -190,8 +190,8 @@ function segment_ring(coords::Vector{Tuple{Float64,Float64}}, great_circle::Bool end # Check if a segment is self-closing -function is_self_closing(segment::Vector{Tuple{Float64,Float64}})::Bool - is_right = segment[end][1] == 180 +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]) @@ -199,23 +199,23 @@ function is_self_closing(segment::Vector{Tuple{Float64,Float64}})::Bool end # Build polygons from segments -function build_polygons(segments::Vector{Vector{Tuple{Float64,Float64}}})::Vector{GI.Polygon} +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] == 180 + is_right = segment[end][1] == c.right candidates = Tuple{Union{Nothing,Int},Float64}[] - if is_self_closing(segment) + 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(s) || s[end][2] < segment[1][2])) || + (!is_self_closing(c, 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])) + (!is_self_closing(c, s) || s[end][2] > segment[1][2])) push!(candidates, (i, s[1][2])) end end @@ -230,10 +230,10 @@ function build_polygons(segments::Vector{Vector{Tuple{Float64,Float64}}})::Vecto # Join segments and recurse segment = vcat(segment, segments[index]) segments[index] = segment - return build_polygons(segments) + return build_polygons(c, segments) else # Handle self-joining segment - polygons = build_polygons(segments) + polygons = build_polygons(c, segments) if !all(p == segment[1] for p in segment) push!(polygons, GI.Polygon([segment])) end @@ -245,9 +245,12 @@ function build_polygons(segments::Vector{Vector{Tuple{Float64,Float64}}})::Vecto end # Main function to cut a polygon at the antimeridian -cut_at_antimeridian(x) = cut_at_antimeridian(GI.trait(x), x) +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, @@ -257,16 +260,16 @@ function cut_at_antimeridian( ) where {T} # Get exterior ring exterior = GO.tuples(GI.getexterior(polygon)).geom - exterior = normalize_coords(exterior) + exterior = normalize_coords(c, exterior) # Segment the exterior ring - segments = segment_ring(exterior, great_circle) + 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.(coord_vecs)) + return GI.Polygon(normalize_coords.((c,), coord_vecs)) end return polygon end @@ -275,13 +278,13 @@ function cut_at_antimeridian( 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) + 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(reverse(hole), great_circle) + hole_segments = segment_ring(c, reverse(hole), great_circle) end end append!(segments, hole_segments) @@ -310,9 +313,9 @@ function cut_at_antimeridian( 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} +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(coords, great_circle) + segments = segment_ring(c, coords, great_circle) if isempty(segments) return line @@ -322,10 +325,10 @@ function cut_at_antimeridian(::GI.AbstractCurveTrait, line::T; great_circle::Boo end -function cut_at_antimeridian(::GI.MultiPolygonTrait, x; kwargs...) +function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.MultiPolygonTrait, x; kwargs...) results = GI.Polygon[] for poly in GI.getgeom(x) - result = cut_at_antimeridian(GI.PolygonTrait(), poly; kwargs...) + result = cut_at_antimeridian(c, GI.PolygonTrait(), poly; kwargs...) if result isa GI.Polygon push!(results, result) elseif result isa GI.MultiPolygon @@ -335,11 +338,11 @@ function cut_at_antimeridian(::GI.MultiPolygonTrait, x; kwargs...) return GI.MultiPolygon(results) end -function cut_at_antimeridian(::GI.MultiLineStringTrait, multiline::T; great_circle::Bool=true) where {T} +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(GI.LineStringTrait(), line; great_circle) + fixed = cut_at_antimeridian(c, GI.LineStringTrait(), line; great_circle) if fixed isa GI.LineString push!(linestrings, GO.tuples(fixed).geom) else From 7c43e2fbcb3199414c9172f4bf0e9839ddee6f96 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 8 Apr 2025 07:37:04 -0400 Subject: [PATCH 3/4] Fix the code --- docs/src/experiments/regridding/regridding.jl | 4 ++-- .../correction/cut_at_antimeridian.jl | 12 +++++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/docs/src/experiments/regridding/regridding.jl b/docs/src/experiments/regridding/regridding.jl index 6833d387e..2a4f0c8fc 100644 --- a/docs/src/experiments/regridding/regridding.jl +++ b/docs/src/experiments/regridding/regridding.jl @@ -93,8 +93,8 @@ grid2 = rand(FullGaussianGrid, 4 + 100) faces1 = SpeedyWeatherGeoMakieExt.get_faces(grid1) faces2 = SpeedyWeatherGeoMakieExt.get_faces(grid2) -polys1 = GI.Polygon.(GI.LinearRing.(eachcol(faces1))) .|> GO.fix -polys2 = GI.Polygon.(GI.LinearRing.(eachcol(faces2))) .|> GO.fix +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) diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl index 2cbc4239f..6a88040ac 100644 --- a/src/transformations/correction/cut_at_antimeridian.jl +++ b/src/transformations/correction/cut_at_antimeridian.jl @@ -54,6 +54,7 @@ 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}) @@ -65,7 +66,8 @@ function _cross(a::Tuple{Float64,Float64,Float64}, b::Tuple{Float64,Float64,Floa end # Convert spherical degrees to cartesian coordinates -function spherical_degrees_to_cartesian(point::Tuple{Float64,Float64})::Tuple{Float64,Float64,Float64} +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) @@ -79,14 +81,14 @@ 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(start) - p2 = spherical_degrees_to_cartesian(endpoint) + 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.left, 0.0) + n2 = spherical_degrees_to_cartesian(c, (c.left, 0.0)) # Intersection of planes defined by cross product intersection = _cross(n1, n2) @@ -294,7 +296,7 @@ function cut_at_antimeridian( end # Build final polygons - result_polygons = build_polygons(segments) + result_polygons = build_polygons(c, segments) # Add holes to appropriate polygons for hole in holes From 3d35b666e6c4071502493bba67f59bb1dc10629c Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 8 Apr 2025 07:37:35 -0400 Subject: [PATCH 4/4] Don't use antimeridian cutting by default it applies itself to planar geoms where it should not --- src/transformations/correction/geometry_correction.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl index eb1ce5737..a9e227f13 100644 --- a/src/transformations/correction/geometry_correction.jl +++ b/src/transformations/correction/geometry_correction.jl @@ -52,7 +52,7 @@ 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[CutAtAntimeridianAndPoles(), ClosedRing()], kwargs...) +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,