Suggested feature: area_intersection(poly1, poly2)

[milankl]

Talking to @asinghvi17, for SpeedyWeather and Oceananigans (@simone-silvestri) we need (towards things like conservative interpolation between two grids on the sphere) we will need

GeometryOps.area(GeometryOps.intersection(GeometryOps.GEOS(), polygon1, polygon2))

to be as fast as possible as polygon1 and 2 will each come from large vectors of polygons (10k to 1mio) and all combinations need to be computed. See also SpeedyWeather/SpeedyWeather.jl#648. Our polygons are all simple and convex, having 4 vertices each, all much smaller than the radius of the sphere. However, we basically never need the actual intersection polygon materialized, we just need its area.

This is an issue to discuss how to get there, maybe an area_intersection(poly1, poly2) function? Or maybe other shortcuts that we could take.

[rafaqz]

It's probably much faster already not using GEOS

[milankl]

You mean just leave that argument out?

[milankl]

Throws an error

MethodError: no method matching TraitTarget(::Nothing)
The type `TraitTarget` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  TraitTarget(::T) where T<:GeoInterface.AbstractTrait
   @ GeometryOps ~/.julia/packages/GeometryOps/xlUi5/src/types.jl:41
  TraitTarget(::GeoInterface.AbstractTrait...)
   @ GeometryOps ~/.julia/packages/GeometryOps/xlUi5/src/types.jl:44
  TraitTarget(::Type{<:TraitTarget{T}}) where T
   @ GeometryOps ~/.julia/packages/GeometryOps/xlUi5/src/types.jl:43
  ...


Stacktrace:
 [1] intersection(geom_a::GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Point{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}, geom_b::GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Point{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}, ::Type{Float64}; target::Nothing, kwargs::@Kwargs{})
   @ GeometryOps ~/.julia/packages/GeometryOps/xlUi5/src/methods/clipping/intersection.jl:45
 [2] intersection(geom_a::GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Point{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}, geom_b::GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Point{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}, ::Type{Float64})
   @ GeometryOps ~/.julia/packages/GeometryOps/xlUi5/src/methods/clipping/intersection.jl:42
 [3] fraction(polygon::GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrapper

[milankl]

@simone-silvestri and I have put our grids on top of each other

to formulate our problem, we have two vectors of polygons on the sphere

grid1 = [poly11, poly12, poly13, ...]
grid2 = [poly21, poly22, poly23, ...]

now we need a matrix of size length(grid1) x length(grid2) where every element is the area of the intersect between a polygon from grid1 with a polygon from grid2. If we get this right, then summing the rows or columns would then yield the area of every polygon in grid1 or grid2. This matrix we can then use to interpolate conservatively from one grid to the other by simply "redistributing" values from one grid cell via the intersects into the overlapping grid cells from the other grid.

Let us consider the following example with a blue grid (4x 1x1 cells) and an orange grid (4 triangles and 1 diamond in the middle) numbered as denoted

the areas would be (calculating by hand)

areas_blue_grid = [1, 1, 1, 1]
areas_orange_grid = [1/2, 1/2, 2, 1/2, 1/2]

and the matrix of the intersect areas would be (calculating by hand!)

julia> using SparseArrays
julia> SparseMatrixCSC([0 0 1/2 1/2 0; 0 0 1/2 0 1/2; 1/2 0 1/2 0 0; 0 1/2 1/2 0 0])
4×5 SparseMatrixCSC{Float64, Int64} with 8 stored entries:
  ⋅    ⋅   0.5  0.5   ⋅
  ⋅    ⋅   0.5   ⋅   0.5
 0.5   ⋅   0.5   ⋅    ⋅
  ⋅   0.5  0.5   ⋅    ⋅

and summing rows/columns would yield the area vectors from above. So for any two grids this matrix is what we need for the conservative interpolation. If we can work together to calculate reasonably fast on the sphere that would be fantastic!!

[milankl]

To illlustrate how to do conservative interpolation between these two grids

1. define the areas

julia> area1 = vec(sum(A, dims=2))
4-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0

julia> area2 = vec(sum(A, dims=1))
5-element Vector{Float64}:
 0.5
 0.5
 2.0
 0.5
 0.5

2. start with some data on grid2 (the orange one)

julia> v2 = [5, 0, 0, 0, 0]
5-element Vector{Int64}:
 5
 0
 0
 0
 0

global integral is 2.5

julia> sum(v2 .* area2)
2.5

3. interpolate onto other grid via multiplication with A and normalize by area on that grid

julia> v1 = (A*v2) ./ area1
4-element Vector{Float64}:
 0.0
 0.0
 2.5
 0.0

global integral is 2.5 again

julia> sum(v1 .* area1)
2.5

4. interpolate back to old grid (this is diffusive mean variance not conserved, but mean should be!) via the transpose A' (no new interpolation operator needed!)

julia> v2back = (A'*v1) ./ area2
5-element Vector{Float64}:
 2.5
 0.0
 0.625
 0.0
 0.0

you can see how half of the 5 we started with diffused into another grid cell (as expected when they overlap). But the integral is still 2.5 🥳

julia> sum(v2back .* area2)
2.5

[simone-silvestri]

Probably interesting also for @gaelforget

[asinghvi17]

A working version of your example on the plane, using #274 for the spatial tree interface and dual tree query function:

import GeoInterface as GI, GeometryOps as GO
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

using SortTileRecursiveTree
using Extents, SparseArrays

function area_of_intersection_operator(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.  
    # 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]
        area_of_intersection = GO.area(GO.intersection(p1, p2; target = GI.PolygonTrait()))
        if area_of_intersection > 0
            A[i1, i2] += area_of_intersection
        end
    end

    return A
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.

[simone-silvestri]

Thanks @asinghvi17 this looks exciting, I'll give it a try with native spherical grids from Oceananigans and SpeedyWeather.

[asinghvi17]

Better implementation in #283