retrieve_bathymetry fails on a 1 degree grid

[glwagner]

MWE:

julia> using Oceananigans, OrthogonalSphericalShellGrids, ClimaOcean
Precompiling ClimaOcean
  3 dependencies successfully precompiled in 9 seconds. 216 already precompiled.

julia> grid = TripolarGrid(size = (360, 170, 1),
                           z = (0, 1),
                           north_poles_latitude = 55,
                           first_pole_longitude = 75)
360×170×1 OrthogonalSphericalShellGrid{Float64, Periodic, RightConnected, Bounded} on CPU with 4×4×4 halo and with precomputed metrics
├── centered at (λ, φ) = (75.0, 1.8005)
├── longitude: Periodic  extent 360.156 degrees     variably spaced with min(Δλ)=0.00798438, max(Δλ)=1.05373
├── latitude:  RightConnected  extent 171.0 degrees variably spaced with min(Δφ)=0.0127345, max(Δφ)=1.00592
└── z:         Bounded  z ∈ [0.0, 1.0]              regularly spaced with Δz=1.0

julia> bathymetry = retrieve_bathymetry(grid)
[ Info: Regridding bathymetry from existing file /Users/gregorywagner/.julia/scratchspaces/0376089a-ecfe-4b0e-a64f-9c555d74d754/Bathymetry/ETOPO_2022_v1_60s_N90W180_surface.nc.
ERROR: ArgumentError: The northern latitude cannot be less than -90 degrees.
Stacktrace:
 [1] validate_lat_lon_grid_args
   @ ~/.julia/packages/Oceananigans/IsF1o/src/Grids/latitude_longitude_grid.jl:281 [inlined]
 [2] LatitudeLongitudeGrid(architecture::CPU, FT::DataType; size::Tuple{…}, longitude::Tuple{…}, latitude::Tuple{…}, z::Tuple{…}, radius::Float64, topology::Nothing, precompute_metrics::Bool, halo::Tuple{…})
   @ Oceananigans.Grids ~/.julia/packages/Oceananigans/IsF1o/src/Grids/latitude_longitude_grid.jl:189
 [3] LatitudeLongitudeGrid
   @ ~/.julia/packages/Oceananigans/IsF1o/src/Grids/latitude_longitude_grid.jl:174 [inlined]
 [4] regrid_bathymetry(target_grid::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{…}; height_above_water::Nothing, minimum_depth::Int64, dir::String, url::String, filename::String, interpolation_passes::Int64, major_basins::Float64)
   @ ClimaOcean.Bathymetry ~/Projects/ClimaOcean.jl/src/Bathymetry.jl:170
 [5] regrid_bathymetry
   @ ~/Projects/ClimaOcean.jl/src/Bathymetry.jl:80 [inlined]
 [6] retrieve_bathymetry(grid::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{…})
   @ ClimaOcean.Bathymetry ~/Projects/ClimaOcean.jl/src/Bathymetry.jl:366
 [7] top-level scope
   @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.

[glwagner]

I'm confused why the grid matters in this situation

[simone-silvestri]

It looks like a precision error. The latitude-longitude grid to load the bathymetry wants to be at latitudes
latitude = (φ₁_data, φ₂_data) = (-80.48333333333333, 90.00000000000001),
and that triggers the error in Oceananigans. It looks like it's not possible to retrieve the bathymetry on a grid that reaches 90 degrees north:

julia> grid = LatitudeLongitudeGrid(size = (100, 100, 1), z = (0, 1), latitude = (-85, 90), longitude = (0, 360))
100×100×1 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×1 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [0.0, 360.0)  regularly spaced with Δλ=3.6
├── latitude:  Bounded  φ ∈ [-85.0, 90.0] regularly spaced with Δφ=1.75
└── z:         Bounded  z ∈ [0.0, 1.0]    regularly spaced with Δz=1.0

julia> bathymetry = retrieve_bathymetry(grid, dir = "./")
[ Info: Regridding bathymetry from existing file ./ETOPO_2022_v1_60s_N90W180_surface.nc.
(φ_data[end], φ_data[end - 1], Δφ) = (89.99166666666667, 89.975, 0.01666666666667993)
(φ₁_data, φ₂_data) = (-85.00000000000001, 90.00000000000001)
ERROR: ArgumentError: The northern latitude cannot be less than -90 degrees.
Stacktrace:
 [1] validate_lat_lon_grid_args
   @ ~/.julia/packages/Oceananigans/s1DfC/src/Grids/latitude_longitude_grid.jl:281 [inlined]
 [2] LatitudeLongitudeGrid(architecture::CPU, FT::DataType; size::Tuple{…}, longitude::Tuple{…}, latitude::Tuple{…}, z::Tuple{…}, radius::Float64, topology::Nothing, precompute_metrics::Bool, halo::Tuple{…})
   @ Oceananigans.Grids ~/.julia/packages/Oceananigans/s1DfC/src/Grids/latitude_longitude_grid.jl:189
 [3] LatitudeLongitudeGrid
   @ ~/.julia/packages/Oceananigans/s1DfC/src/Grids/latitude_longitude_grid.jl:174 [inlined]
 [4] regrid_bathymetry(target_grid::LatitudeLongitudeGrid{…}; height_above_water::Nothing, minimum_depth::Int64, dir::String, url::String, filename::String, interpolation_passes::Int64, major_basins::Float64)
   @ ClimaOcean.Bathymetry ~/development/ClimaOcean.jl/src/Bathymetry.jl:173
 [5] regrid_bathymetry
   @ ~/development/ClimaOcean.jl/src/Bathymetry.jl:80 [inlined]
 [6] #retrieve_bathymetry#13
   @ ~/development/ClimaOcean.jl/src/Bathymetry.jl:366 [inlined]
 [7] top-level scope
   @ REPL[44]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia>

I wonder how we did not have hit this error before or if something did not change in the bathymetry.

[simone-silvestri]

This error happens because the last element in

ClimaOcean.jl/src/Bathymetry.jl

Line 114 in 69a111a

φ_data = dataset["lat"][:]

is

89.99166666666667

While it should be

89.99166666666666

to have the face of the last element at 90 instead of 90.00000000000001

[glwagner]

Yeah weird that it somehow depends on the grid?

But maybe the constructor for the lat lon grid should be allowed to fudge the bounds by very small amounts to ensure correctness?

[simone-silvestri]

It seems like changing

ClimaOcean.jl/src/Bathymetry.jl

Lines 114 to 115 in 69a111a

	φ_data = dataset["lat"][:]
	λ_data = dataset["lon"][:]

to

    φ_data = dataset["lat"][:] |> Array{BigFloat}
    λ_data = dataset["lon"][:] |> Array{BigFloat}

ClimaOcean.jl/src/Bathymetry.jl

Lines 160 to 165 in 69a111a

	λ₁_data = λ_data[1] - Δλ / 2
	λ₂_data = λ_data[end] + Δλ / 2
	φ₁_data = φ_data[1] - Δφ / 2
	φ₂_data = φ_data[end] + Δφ / 2

to

    λ₁_data = convert(Float64, λ_data[1]   - Δλ / 2)
    λ₂_data = convert(Float64, λ_data[end] + Δλ / 2)
    φ₁_data = convert(Float64, φ_data[1]   - Δφ / 2)
    φ₂_data = convert(Float64, φ_data[end] + Δφ / 2)

solves the issue.

[simone-silvestri]

But maybe the constructor for the lat lon grid should be allowed to fudge the bounds by very small amounts to ensure correctness?

Agreed, we can do something like this:

φ₂ <= 90 + tolerance  || throw(ArgumentError("The northern latitude cannot be less than -90 degrees."))
φ₂ = φ₂ > 90 ? 90 : φ₂

(and the same for the southern boundary)

[simone-silvestri]

However I am wondering how

ClimaOcean.jl/test/test_simulations.jl

Lines 8 to 30 in 69a111a

	grid = TripolarGrid(arch;
	size = (50, 50, 10),
	halo = (7, 7, 7),
	z = (-6000, 0),
	first_pole_longitude = 75,
	north_poles_latitude = 55)
	bottom_height = retrieve_bathymetry(grid;
	minimum_depth = 10,
	dir = "./",
	interpolation_passes = 20,
	major_basins = 1)
	grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true)
	free_surface = SplitExplicitFreeSurface(grid; substeps = 20)
	ocean = ocean_simulation(grid; free_surface)
	backend = JRA55NetCDFBackend(4)
	atmosphere = JRA55_prescribed_atmosphere(arch; backend)
	radiation = Radiation(arch)
	sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce()

manages to run without errors

[glwagner]

I think directly / manually fixing the rounding issue might be preferred to BigFloat, but it does seem like BigFloat will work for bathymetry at least, because we always interpolate to a new field on the model grid? For other situations where we need to interpolate on the fly, we'd need a different solution.