regrid_bathymetry does not support RotatedLatitudeLongitudeGrid #666
Description
regrid_bathymetry does not support RotatedLatitudeLongitudeGrid
Issue Description
The regrid_bathymetry function fails when used with a RotatedLatitudeLongitudeGrid (which is an OrthogonalSphericalShellGrid type). This prevents users from loading realistic bathymetry for rotated grids, which are useful for high-latitude ocean simulations to avoid pole singularities.
Minimal Reproducible Example
using Oceananigans
using ClimaOcean
using ClimaOcean.OrthogonalSphericalShellGrids: RotatedLatitudeLongitudeGrid
# Create a rotated grid for Arctic region (60°N - 85°N)
grid = RotatedLatitudeLongitudeGrid(
size = (180, 25, 30),
longitude = (-180, 180),
latitude = (60, 85),
z = (-2000, 0),
north_pole = (0, 90)
)
# Attempt to load bathymetry - this fails
bottom_height = regrid_bathymetry(grid,
interpolation_passes = 5,
minimum_depth = 10
)Error Message
ERROR: unknown target grid type
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] interpolate_bathymetry_in_passes(native_z::Field{…}, target_grid::OrthogonalSphericalShellGrid{…}; passes::Int64)
@ ClimaOcean.Bathymetry ~/.julia/packages/ClimaOcean/i4kd3/src/Bathymetry.jl:207
[3] _regrid_bathymetry(target_grid::OrthogonalSphericalShellGrid{…}, metadata::Metadatum{…}; height_above_water::Nothing, minimum_depth::Int64, interpolation_passes::Int64, major_basins::Int64)
@ ClimaOcean.Bathymetry ~/.julia/packages/ClimaOcean/i4kd3/src/Bathymetry.jl:125
Root Cause
In src/Bathymetry.jl, lines 202-207, the interpolate_bathymetry_in_passes function performs a type check that doesn't include OrthogonalSphericalShellGrid:
gridtype = target_grid isa TripolarGrid ? "TripolarGrid" :
target_grid isa LatitudeLongitudeGrid ? "LatitudeLongitudeGrid" :
target_grid isa RectilinearGrid ? "RectilinearGrid" :
error("unknown target grid type")RotatedLatitudeLongitudeGrid is defined as an OrthogonalSphericalShellGrid (see OrthogonalSphericalShellGrids/rotated_latitude_longitude_grid.jl), which is not in this type check list. However, the underlying interpolation machinery (interpolate! function) should work fine with this grid type since it's sufficiently generic.
Proposed Solution
Add OrthogonalSphericalShellGrid to the type check in src/Bathymetry.jl:
gridtype = target_grid isa TripolarGrid ? "TripolarGrid" :
target_grid isa LatitudeLongitudeGrid ? "LatitudeLongitudeGrid" :
target_grid isa RectilinearGrid ? "RectilinearGrid" :
target_grid isa OrthogonalSphericalShellGrid ? "RotatedLatitudeLongitudeGrid" :
error("unknown target grid type")Note: The gridtype variable is only used for logging messages, not for any actual logic, so this change should be safe and backward compatible.
Workaround
Users can work around this issue by overloading the function in their own code:
using Oceananigans.Architectures: architecture
using Oceananigans.Grids: topology, x_domain, y_domain
using Oceananigans.Fields: interpolate!, Field
import ClimaOcean.Bathymetry: interpolate_bathymetry_in_passes
function interpolate_bathymetry_in_passes(native_z, target_grid::OrthogonalSphericalShellGrid;
passes = 10)
Nλt, Nφt = Nt = size(target_grid)
Nλn, Nφn = Nn = size(native_z)
latitude = y_domain(native_z.grid)
longitude = x_domain(native_z.grid)
ΔNλ = floor((Nλn - Nλt) / passes)
ΔNφ = floor((Nφn - Nφt) / passes)
Nλ = [Nλn - ΔNλ * pass for pass in 1:passes-1]
Nφ = [Nφn - ΔNφ * pass for pass in 1:passes-1]
Nλ = Int[Nλ..., Nλt]
Nφ = Int[Nφ..., Nφt]
old_z = native_z
TX, TY = topology(target_grid)
Hx, Hy, Hz = Oceananigans.halo_size(native_z.grid)
@info "Interpolation passes of bathymetry size $(size(old_z)) onto a RotatedLatitudeLongitudeGrid target grid of size $Nt:"
for pass = 1:passes - 1
new_size = (Nλ[pass], Nφ[pass], 1)
@info " pass $pass to size $new_size"
new_grid = LatitudeLongitudeGrid(architecture(target_grid), Float32,
size = new_size,
latitude = (latitude[1], latitude[2]),
longitude = (longitude[1], longitude[2]),
z = (0, 1),
topology = (TX, TY, Bounded),
halo = (Hx, Hy, Hz))
new_z = Field{Center, Center, Nothing}(new_grid)
interpolate!(new_z, old_z)
old_z = new_z
end
new_size = (Nλ[passes], Nφ[passes], 1)
@info " pass $passes to size $new_size"
target_z = Field{Center, Center, Nothing}(target_grid)
interpolate!(target_z, old_z)
return target_z
endEnvironment
- Julia version: [v1.10.10]
- ClimaOcean version: [0.8.4]
- Oceananigans version: [0.97.6, 0.98]
- OS: [Linux]
Use Case
This feature is needed for Arctic and Antarctic ocean simulations where rotated grids help avoid numerical issues at the poles while still using realistic bathymetry data. The rotated grid approach is particularly important for high-latitude regions (e.g., 60°N - 85°N) where standard latitude-longitude grids suffer from coordinate singularities.