Trying to plot heat map of field in halo region using Makie gives `BoundsError`

[matt-graham]

I am not sure if this is intended behaviour or not, but when trying to plot a heatmap of a field with indices corresponding to a slice in the halo region, I get a BoundsError due to the line

https://github.com/CliMA/Oceananigans.jl/blob/df5150aebcdeb9bdbeddbbf5dd91cef0cc1d9a0e/ext/OceananigansMakieExt.jl#L149

which gets nodes for field, as with_halos=true is not set.

Minimal reproducing example

using Oceananigans
using CairoMakie

grid = RectilinearGrid(size=(5, 5, 5), halo=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1))
field = Field{Center, Center, Center}(grid, indices=(:, :, grid.Nz + 1))
heatmap(field)

gives error / stack trace

ERROR: BoundsError: attempt to access 5-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64} at index [6:6]
Stacktrace:
  [1] throw_boundserror(A::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, I::Tuple{UnitRange{Int64}})
    @ Base ./essentials.jl:15
  [2] checkbounds
    @ ./abstractarray.jl:699 [inlined]
  [3] getindex(r::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, s::UnitRange{Int64})
    @ Base ./twiceprecision.jl:496
  [4] view
    @ ./subarray.jl:234 [inlined]
  [5] rnodes
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Grids/vertical_discretization.jl:201 [inlined]
  [6] rnodes
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Grids/vertical_discretization.jl:202 [inlined]
  [7] znodes
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Grids/vertical_discretization.jl:213 [inlined]
  [8] #nodes#71
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Grids/rectilinear_grid.jl:489 [inlined]
  [9] nodes
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Grids/rectilinear_grid.jl:486 [inlined]
 [10] nodes
    @ ~/.julia/packages/Oceananigans/VsQ7B/src/Fields/field.jl:863 [inlined]
 [11] convert_field_argument(f::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Nothing})
    @ OceananigansMakieExt ~/.julia/packages/Oceananigans/VsQ7B/ext/OceananigansMakieExt.jl:149
 [12] _create_plot(F::Function, attributes::Dict{Symbol, Any}, f::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Nothing})
    @ OceananigansMakieExt ~/.julia/packages/Oceananigans/VsQ7B/ext/OceananigansMakieExt.jl:55
 [13] #heatmap#49
    @ ~/.julia/packages/Makie/TOy8O/src/recipes.jl:530 [inlined]
 [14] heatmap(args::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Nothing})
    @ Makie ~/.julia/packages/Makie/TOy8O/src/recipes.jl:528
 [15] top-level scope
    @ REPL[41]:1
Some type information was truncated. Use `show(err)` to see complete types.

[glwagner]

Is this maybe a problem with nodes? I think @tomchor was working on making nodes respect Field.indices.

[tomchor]

Is this maybe a problem with nodes? I think @tomchor was working on making nodes respect Field.indices.

Yes. This was already merged in https://github.com/CliMA/Oceananigans.jl/pull/4814. @matt-graham have you tried with the latest version?

[glwagner]

@tomchor did you specifically treat / test the case that the indices involve halo regions? I guess we interpret (:, :, :) as "interior", but then if we have (:, :, Nz+1) this might require a z slice inside the halo regions.

[matt-graham]

Yes. This was already merged in #4814. @matt-graham have you tried with the latest version?

Can confirm the above minimal example still errors for me with v0.102.1 with same stack trace as above.

[glwagner]

Yes. This was already merged in #4814. @matt-graham have you tried with the latest version?

Can confirm the above minimal example still errors for me with v0.102.1 with same stack trace as above.

Can you make a more minimal example that isolates the nodes issue without requiring CairoMakie? Will help for debugging and solving! Also, this must not be tested.

[matt-graham]

Yes. This was already merged in #4814. @matt-graham have you tried with the latest version?

Can confirm the above minimal example still errors for me with v0.102.1 with same stack trace as above.

Can you make a more minimal example that isolates the nodes issue without requiring CairoMakie? Will help for debugging and solving! Also, this must not be tested.

So the following also reproduces the error without requiring CairoMakie

using Oceananigans

grid = RectilinearGrid(size=(5, 5, 5), halo=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1))
field = Field{Center, Center, Center}(grid, indices=(:, :, grid.Nz + 1))
nodes(field)

and is directly analogous to the code path from which the error originates in the CairoMakie example. However in this case, I would say that its unclear if the behaviour is a bug, as there is a with_halos argument to nodes which defaults to false and so it would seem the 'correct' example in that case would be

using Oceananigans

grid = RectilinearGrid(size=(5, 5, 5), halo=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1))
field = Field{Center, Center, Center}(grid, indices=(:, :, grid.Nz + 1))
nodes(field, with_halos=true)

which returns a 3-tuple ([-0.1, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1], [-0.1, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1], [1.1]) of node values as expected. I guess it would be nice if the error reported in the original case specifically pointed out that using with_halos=true would resolve the error and the effect of with_halos was documented, but otherwise it seems the current behaviour is probably as expected?

For my use case the specific issue is that the nodes call in convert_field_argument in OceananigansMakieExt does not set with_halos=true and there isn't anyway to passthrough this keyword argument. A more faithful minimal example might therefore be something like

using Oceananigans

grid = RectilinearGrid(size=(5, 5, 5), halo=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1))
field = Field{Center, Center, Center}(grid, indices=(:, :, grid.Nz + 1))
convert_field_argument(field)

but I don't know how to load and use functionality in an extension module (to allow calling convert_field_argument) without explicitly importing the relevant dependency modules.

[tomchor]

@tomchor did you specifically treat / test the case that the indices involve halo regions? I guess we interpret (:, :, :) as "interior", but then if we have (:, :, Nz+1) this might require a z slice inside the halo regions.

I don't remember explicitly adding tests for this. Although if this was tested before, it should still be tested now since I didn't remove any tests.

[tomchor]

Thanks for the MWE @matt-graham. I'll take a look into it once I get back in office. As @glwagner said, we should test this (and any other potential cases the tests are missing).

As a side note, the Makie extension has very bare bones testing. We should probably bulk that sometime soon.

[glwagner]

However in this case, I would say that its unclear if the behaviour is a bug, as there is a with_halos argument to nodes

My opinion is that it is a bug. In fact, I think we should deprecate the with_halos argument. With a plotting extension and the ability to use nodes in abstract operations, there should be no need for users to call nodes themselves.