Meridian cutting issue

[simone-silvestri]

I think we are still hitting a meridian-cutting issue. Here is a minimum working example (unfortunately I need to write here the methods to get the faces)

using Oceananigans
using Oceananigans.Grids: λnode, φnode
using SpeedyWeather
using ConservativeRegridding

function get_faces(grid::Oceananigans.AbstractGrid)
    Nx, Ny, _ = size(grid)
    cell_matrix = Array{Tuple{Float64, Float64}}(undef, 5, Nx*Ny)

    @inbounds for i in 1:Nx, j in 1:Ny
        linear_idx = i + (j - 1) * Nx
    
        λ⁻⁻ = λnode(i,   j,   1, grid, Face(), Face(), nothing)
        λ⁺⁻ = λnode(i,   j+1, 1, grid, Face(), Face(), nothing)
        λ⁻⁺ = λnode(i+1, j,   1, grid, Face(), Face(), nothing)
        λ⁺⁺ = λnode(i+1, j+1, 1, grid, Face(), Face(), nothing)
        φ⁻⁻ = φnode(i,   j,   1, grid, Face(), Face(), nothing)
        φ⁺⁻ = φnode(i,   j+1, 1, grid, Face(), Face(), nothing)
        φ⁻⁺ = φnode(i+1, j,   1, grid, Face(), Face(), nothing)
        φ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), nothing)

        cell_matrix[1, linear_idx] = (λ⁻⁻, φ⁻⁻)
        cell_matrix[2, linear_idx] = (λ⁻⁺, φ⁻⁺)
        cell_matrix[3, linear_idx] = (λ⁺⁺, φ⁺⁺)
        cell_matrix[4, linear_idx] = (λ⁺⁻, φ⁺⁻)
        cell_matrix[5, linear_idx] = (λ⁻⁻, φ⁻⁻)
    end
    
    return cell_matrix
end

function get_faces(Grid, nlat_half::Integer)
    npoints = RingGrids.get_npoints2D(Grid, nlat_half)
    E, S, W, N = RingGrids.get_vertices(Grid, nlat_half)
    faces = Matrix{NTuple{2, Float64}}(undef, 5, npoints)

    @inbounds for ij in 1:npoints
        faces[1, ij] = (E[1, ij], E[2, ij])  # clockwise
        faces[2, ij] = (S[1, ij], S[2, ij])
        faces[3, ij] = (W[1, ij], W[2, ij])
        faces[4, ij] = (N[1, ij], N[2, ij])
        faces[5, ij] = (E[1, ij], E[2, ij])  # back to east to close the polygon        
    end

    return faces
end

get_faces(grid::SpeedyWeather.AbstractGridArray) = get_faces(typeof(grid), grid.nlat_half)

grid = LatitudeLongitudeGrid(Oceananigans.CPU(); size=(50, 50, 1), latitude=(-80, 80), longitude=(0, 360), z = (0, 1))

src = rand(OctaHEALPixGrid, 5 + 100) .+ 100
dst = CenterField(grid)

src_cells = get_faces(src)
dst_cells = get_faces(grid)

regridder = ConservativeRegridding.Regridder(dst_cells, src_cells)
ConservativeRegridding.regrid!(vec(interior(dst)), regridder, src)

heatmap(dst)
heatmap(src)

The output of this regridding is:

[simone-silvestri]

With @milankl we are looking at the areas that are computed by compute_intersection_areas for
rand(OctahedralGaussianGrid, 24) and we get the following:

So it looks like there is a meridian cutting problem at 0-360

EDIT: Note this is the OctahedralGaussianGrid which has cell centres on the prime meridian, so there's polygons with vertices that are on either side of 0˚/360˚E, these are the ones you can see in yellow

[simone-silvestri]

This is what happens for a OctaminimalGaussianGrid, EDIT by milan: This is a grid which has cells with vertices on the meridian but not the centres, not sure why the area "blows" up sometimes east/west of the meridian

[glwagner]

@simone-silvestri @asinghvi17 what is the status of this issue? Is there anything I can do to help resolve it?

[simone-silvestri]

I am not sure, I think the status is the same as my last post, but I think it will depend on the last commits in this PR JuliaGeo/GeometryOps.jl#283.

[glwagner]

There haven't been any commits in a few months on that PR. @asinghvi17 is there any way I can help?

[glwagner]

@milankl may have information too. This is blocking some work on coupling which we hope to make progress on by the end of this summer.

[simone-silvestri]

I tried again the MWE at the top of the issue, updating the packages, and now I encounter the following

julia> regridder = ConservativeRegridding.Regridder(dst_cells, src_cells)
ERROR: UndefVarError: `do_dual_query` not defined
Stacktrace:
 [1] getproperty
   @ ./Base.jl:31 [inlined]
 [2] compute_intersection_areas!(areas::SparseArrays.SparseMatrixCSC{…}, grid1::Vector{…}, grid2::Vector{…}; manifold::GeometryOpsCore.Planar, nodecapacity1::Int64, nodecapacity2::Int64)
   @ ConservativeRegridding ~/development/ConservativeRegridding.jl/src/regridder.jl:73
 [3] compute_intersection_areas!
   @ ~/development/ConservativeRegridding.jl/src/regridder.jl:56 [inlined]
 [4] #intersection_areas#1
   @ ~/development/ConservativeRegridding.jl/src/regridder.jl:53 [inlined]
 [5] intersection_areas(src_field::Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{…}, Nothing, Nothing}}, dst_field::Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{…}, Nothing, Nothing}})
   @ ConservativeRegridding ~/development/ConservativeRegridding.jl/src/regridder.jl:44
 [6] ConservativeRegridding.Regridder(dst_vertices::Matrix{Tuple{Float64, Float64}}, src_vertices::Matrix{Tuple{Float64, Float64}}; normalize::Bool, kwargs::@Kwargs{})
   @ ConservativeRegridding ~/development/ConservativeRegridding.jl/src/regridder.jl:116
 [7] ConservativeRegridding.Regridder(dst_vertices::Matrix{Tuple{Float64, Float64}}, src_vertices::Matrix{Tuple{Float64, Float64}})
   @ ConservativeRegridding ~/development/ConservativeRegridding.jl/src/regridder.jl:106
 [8] top-level scope
   @ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types.

[asinghvi17]

Ah we need to update that. I haven't had much time to look into this in the last few months but should have some next week! Will keep you posted

[asinghvi17]

The function was changed to dual_depth_first_search iirc

[simone-silvestri]

Ah we need to update that. I haven't had much time to look into this in the last few months but should have some next week! Will keep you posted

Sounds good, thanks @asinghvi17!