From e435a0ca81d8ce81d24c98d554c785d0ba46c035 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 24 Mar 2026 11:24:24 +0100 Subject: [PATCH 01/11] Add ocean basin labeling and masking from ClimaOcean Port the refactored labeling framework from ClimaOcean while preserving the existing bathymetry caching infrastructure. Adds Barrier type for basin separation, label_ocean_basins with periodic/tripolar support, OceanBasinMask with predefined seeds and convenience functions for Atlantic, Indian, Pacific, Southern, and Arctic oceans. Co-Authored-By: Claude Opus 4.6 (1M context) --- src/Bathymetry/Bathymetry.jl | 8 +- src/Bathymetry/label_ocean_basins.jl | 226 ++++++++++++++++ src/Bathymetry/ocean_basin_mask.jl | 375 +++++++++++++++++++++++++++ src/Bathymetry/regrid_bathymetry.jl | 93 ++----- src/NumericalEarth.jl | 8 + test/test_bathymetry.jl | 113 +++++++- 6 files changed, 755 insertions(+), 68 deletions(-) create mode 100644 src/Bathymetry/label_ocean_basins.jl create mode 100644 src/Bathymetry/ocean_basin_mask.jl diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index 40eed858d..1aad80397 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -1,6 +1,10 @@ module Bathymetry export regrid_bathymetry, ORCAGrid +export OceanBasinMask +export atlantic_ocean_mask, indian_ocean_mask, southern_ocean_mask, pacific_ocean_mask, arctic_ocean_mask +export label_ocean_basins +export Barrier using Downloads using ImageMorphology @@ -12,7 +16,7 @@ using Oceananigans.BoundaryConditions using Oceananigans.DistributedComputations using Oceananigans.DistributedComputations: DistributedGrid, reconstruct_global_grid, all_reduce, @root using Oceananigans.Fields: interpolate! -using Oceananigans.Grids: x_domain, y_domain, topology +using Oceananigans.Grids: x_domain, y_domain, topology, AbstractGrid using Oceananigans.Utils: launch! using OffsetArrays using NCDatasets @@ -22,7 +26,9 @@ using Scratch using ..DataWrangling: Metadatum, native_grid, metadata_path, download_dataset using ..DataWrangling.ETOPO: ETOPO2022 +include("label_ocean_basins.jl") include("regrid_bathymetry.jl") +include("ocean_basin_mask.jl") include("orca_grid.jl") end # module diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl new file mode 100644 index 000000000..1416bbce8 --- /dev/null +++ b/src/Bathymetry/label_ocean_basins.jl @@ -0,0 +1,226 @@ +using Oceananigans.OrthogonalSphericalShellGrids: TripolarGridOfSomeKind +using Oceananigans.Fields: convert_to_0_360 + +##### +##### Barrier type for separating ocean basins +##### + +""" + Barrier{W, E, S, N} + +A rectangular geographic region used to separate ocean basins during labeling. +All cells within this region are temporarily marked as land. + +Fields +====== +- `west`: Western longitude limit (degrees) +- `east`: Eastern longitude limit (degrees) +- `south`: Southern latitude limit (degrees) +- `north`: Northern latitude limit (degrees) + +Constructors +============ +- `Barrier(west, east, south, north)`: Create a barrier with explicit bounds +- `Barrier(; west, east, south, north)`: Keyword argument version +""" +struct Barrier{W, E, S, N} + west :: W + east :: E + south :: S + north :: N +end + +Barrier(; west, east, south, north) = Barrier(west, east, south, north) + +""" + Barrier(longitude, south, north; width=2.0) + +Create a meridional (north-south) barrier at a given longitude. +Useful for closing straits like Cape Agulhas or Indonesian passages. +""" +Barrier(longitude, south, north; width=2.0) = Barrier(longitude - width/2, longitude + width/2, south, north) + +##### +##### Barrier application functions +##### + +""" + apply_barrier!(zb_data, grid, barrier::Barrier) + +Mark all cells within the barrier region as land (z = 0). +""" +apply_barrier!(zb, grid, barrier::Barrier) = + launch!(architecture(grid), grid, :xy, _apply_barrier!, zb, grid, barrier) + +apply_barrier!(zb, grid, barriers::Nothing) = zb + +function apply_barrier!(zb, grid, barriers::AbstractVector) + for barrier in barriers + apply_barrier!(zb, grid, barrier) + end + return zb +end + +@kernel function _apply_barrier!(zb, grid, barrier::Barrier) + i, j = @index(Global, NTuple) + + in_lon = if isnothing(barrier.west) || (barrier.east - barrier.west >= 360) + true + else + bw = convert_to_0_360(barrier.west) + be = convert_to_0_360(barrier.east) + λ = λnode(i, j, 1, grid, Center(), Center(), Center()) + λ = convert_to_0_360(λ) + (bw <= λ <= be) + end + + φ = φnode(i, j, 1, grid, Center(), Center(), Center()) + in_lat = barrier.south <= φ <= barrier.north + + @inbounds zb[i, j, 1] = ifelse(in_lon & in_lat, zero(grid), zb[i, j, 1]) +end + +# Since the strel algorithm in `remove_major_basins` does not recognize periodic boundaries, +# before removing connected regions, we extend the longitude direction if it is periodic. +# An extension of half the domain is enough. +function maybe_extend_longitude(zb_cpu::Field, ::Periodic) + Nx = size(zb_cpu, 1) + nx = Nx ÷ 2 + + zb_data = zb_cpu.data[1:Nx, :, 1] + zb_parent = zb_data.parent + + # Add information on the LHS and to the RHS + zb_parent = vcat(zb_parent[nx:Nx, :], zb_parent, zb_parent[1:nx, :]) + + # Update offsets: the original domain starts at parent index (Nx - nx + 2) + # For offset index 1 to map there, we need offset = -(nx + 1) + yoffsets = zb_cpu.data.offsets[2] + xoffsets = - nx - 1 + + return OffsetArray(zb_parent, xoffsets, yoffsets) +end + +maybe_extend_longitude(zb_cpu::Field, tx) = interior(zb_cpu, :, :, 1) + +""" + label_ocean_basins(zb_field, TX, core_size) + +Creates a matrix with a unique label for each connected basin. Useful for inpainting the bathymetry and +computing the masks for oceanic basins. + +Handles periodic boundary extension internally and returns labels for the core region only. +""" +function label_ocean_basins(zb_field, TX, size) + zb = maybe_extend_longitude(zb_field, TX()) # Outputs a 2D AbstractArray + + water = zb .< 0 + + connectivity = ImageMorphology.strel(water) + labels = ImageMorphology.label_components(connectivity) + + Nx, Ny = size[1], size[2] + core_labels = labels[1:Nx, 1:Ny] + + # Enforce periodicity: merge labels that should be connected across + # the periodic boundaries. This handles cases where a barrier (e.g., blocking + # the Southern Ocean) prevents the extended domain from connecting basins + # that are actually periodic neighbors. + enforce_periodic_labels!(core_labels, TX()) + enforce_tripolar_labels!(core_labels, zb_field.grid) + + return core_labels +end + +""" + enforce_periodic_labels!(labels, ::Periodic) + +Merge labels that should be connected due to periodic boundary conditions. +For each latitude, if the westernmost and easternmost cells are both water +(non-zero labels), they are periodic neighbors and must have the same label. +""" +function enforce_periodic_labels!(labels, ::Periodic) + Nx, Ny = Base.size(labels) + + for j in 1:Ny + label_west = labels[1, j] + label_east = labels[Nx, j] + + # Both cells are water and have different labels: merge them + if label_west != 0 && label_east != 0 && label_west != label_east + # Replace all occurrences of label_east with label_west + replace!(labels, label_east => label_west) + end + end + + return labels +end + +# No-op for non-periodic domains +enforce_periodic_labels!(labels, tx) = labels + +""" + enforce_tripolar_labels!(labels, ::TripolarGridOfSomeKind) + +Merge labels that should be connected due to zipper boundary conditions. +For each longitude, if cells reflected across the fold (Nx÷2) are water +(non-zero labels), they are periodic neighbors and must have the same label. +""" +function enforce_tripolar_labels!(labels, ::TripolarGridOfSomeKind) + Nx, Ny = Base.size(labels) + + for i in 1:Nx÷2 + label_west = labels[i, Ny] + label_east = labels[Nx-i+1, Ny] + + # Both cells are water and have different labels: merge them + if label_west != 0 && label_east != 0 && label_west != label_east + # Replace all occurrences of label_east with label_west + replace!(labels, label_east => label_west) + end + end + + return labels +end + +# No-op for non-tripolar grids +enforce_tripolar_labels!(labels, grid) = labels + +# Utilities to label ocean basins passing only the grid +function label_ocean_basins(grid::AbstractGrid; barriers=nothing) + @warn "The grid is not immersed, there is only one ocean basin!" + Nx, Ny, Nz = size(grid) + return zeros(Int, Nx, Ny) +end + +""" + label_ocean_basins(grid::ImmersedBoundaryGrid; barriers=nothing) + +Label connected ocean basins in an ImmersedBoundaryGrid. + +Keyword Arguments +================= +- `barriers`: Collection of barriers to apply before labeling. Barriers temporarily + mark certain cells as land, allowing separation of connected ocean basins + (e.g., separating Atlantic from Pacific via the Southern Ocean). +""" +function label_ocean_basins(grid::ImmersedBoundaryGrid; barriers=nothing) + + # The labelling algorithm works only on CPUs + cpu_grid = on_architecture(CPU(), grid) + + TX = topology(cpu_grid, 1) + zb = cpu_grid.immersed_boundary.bottom_height + + # If barriers are specified, apply them to a copy of the bathymetry + if !isnothing(barriers) + # Create a temporary field with the modified bathymetry + zb_modified = Field{Center, Center, Nothing}(cpu_grid) + parent(zb_modified) .= parent(zb) + apply_barrier!(zb_modified, cpu_grid, barriers) + else + zb_modified = zb + end + + return label_ocean_basins(zb_modified, TX, size(cpu_grid)) +end diff --git a/src/Bathymetry/ocean_basin_mask.jl b/src/Bathymetry/ocean_basin_mask.jl new file mode 100644 index 000000000..451610e31 --- /dev/null +++ b/src/Bathymetry/ocean_basin_mask.jl @@ -0,0 +1,375 @@ +using Oceananigans.Grids: λnode, φnode, on_architecture +using ImageMorphology + +##### +##### OceanBasinMask struct +##### + +""" + OceanBasinMask{M, G} + +A mask identifying cells belonging to a specific ocean basin. +The mask is stored as a 2D Field{Center, Center, Nothing} with values: +- 1.0: cell belongs to the basin +- 0.0: cell does not belong to the basin +""" +struct OceanBasinMask{M, G, S} + mask :: M + grid :: G + seed_points :: S +end + +Base.summary(obm::OceanBasinMask) = "OceanBasinMask" + +function Base.show(io::IO, obm::OceanBasinMask) + print(io, summary(obm), " on ", summary(obm.grid)) +end + +# Forward getindex to mask +Base.getindex(obm::OceanBasinMask, i, j, k) = obm.mask[i, j, k] + +##### +##### Additive basin masks +##### + +""" + Base.:+(mask1::OceanBasinMask, mask2::OceanBasinMask) + +Combine two ocean basin masks, returning a new mask that is the union of both. +Cells that belong to either basin will have value 1.0. +""" +function Base.:+(mask1::OceanBasinMask, mask2::OceanBasinMask) + grid = mask1.grid + combined_mask = Field{Center, Center, Nothing}(grid) + + # Union: cell is in combined mask if it's in either mask + parent(combined_mask) .= max.(parent(mask1.mask), parent(mask2.mask)) + fill_halo_regions!(combined_mask) + + # Combine seed points + combined_seeds = (mask1.seed_points..., mask2.seed_points...) + + return OceanBasinMask(combined_mask, grid, combined_seeds) +end + +##### +##### Connected component utilities +##### + +""" + find_label_at_point(labels, grid, λs, φs; radius = 2) + +Find the connected component label at the given longitude/latitude seed point. +Returns the label value, or 0 if the point is on land or outside the domain. +Checks in a circular cap of radius `radius` degrees around the seed point +""" +function find_label_at_point(labels, grid, λs, φs; radius = 2) + Nx, Ny, _ = size(grid) + + # Find grid cell containing the seed point + for j in 1:Ny + for i in 1:Nx + λ = λnode(i, j, 1, grid, Center(), Center(), Center()) + φ = φnode(i, j, 1, grid, Center(), Center(), Center()) + + λ = convert_to_0_360(λ) + + Δλ = if isnothing(λs) + zero(λ) + else + λs = convert_to_0_360(λs) + λ - λs + end + + # Check if this cell contains the seed point (within a certain radius) + if Δλ^2 + (φ - φs)^2 < radius^2 + return labels[i, j] + end + end + end + + return 0 # Seed point not found +end + +""" + create_basin_mask_from_label(grid, labels, basin_label) + +Create a mask field for all cells with the given label. +""" +function create_basin_mask_from_label(grid, labels, basin_label) + arch = architecture(grid) + Nx, Ny, _ = size(grid) + FT = eltype(grid) + + mask = Field{Center, Center, Nothing}(grid) + + launch!(CPU(), grid, :xy, _compute_basin_mask!, + mask, grid, labels, basin_label) + + fill_halo_regions!(mask) + + return mask +end + +@kernel function _compute_basin_mask!(mask, grid, labels, basin_label) + i, j = @index(Global, NTuple) + correct_basin = @inbounds labels[i, j] == basin_label + @inbounds mask[i, j, 1] = correct_basin +end + +##### +##### Some useful Basin seeds and barriers +##### + +const SOUTHERN_OCEAN_SEPARATION_BARRIER = Barrier(-180.0, 180.0, -56.0, -54.0) + +const ATLANTIC_OCEAN_BARRIERS = [ + Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) + Barrier(289.0, -90.0, -30.0) # Drake passage (meridional barrier) +] + +const INDIAN_OCEAN_BARRIERS = [ + Barrier(141.0, -90.0, -3.0), # Indonesian side (meridional) + Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) + Barrier(105.0, 141.0, -4.0, -3.0), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) +] + +const SOUTHERN_OCEAN_BARRIERS = [SOUTHERN_OCEAN_SEPARATION_BARRIER] + +const PACIFIC_OCEAN_BARRIERS = [ + Barrier(141.0, -90.0, -3.0), # Indonesian side (meridional) + Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) + Barrier(105.0, 141.0, -4.0, -3.0), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) +] + +# Seed points for Atlantic Ocean (definitely in the Atlantic) +const ATLANTIC_SEED_POINTS = [ + (-30.0, 0.0), # Central equatorial Atlantic + (-40.0, 30.0), # North Atlantic + (-25.0, -20.0), # South Atlantic +] + +# Seed points for Indian Ocean +const INDIAN_SEED_POINTS = [ + (70.0, -10.0), # Central Indian Ocean + (60.0, 10.0), # Arabian Sea region + (90.0, -20.0), # Eastern Indian Ocean +] + +# Seed points for Southern Ocean +const SOUTHERN_SEED_POINTS = [ + (0.0, -60.0), # South Atlantic sector + (90.0, -60.0), # Indian Ocean sector + (180.0, -60.0), # Pacific sector (date line) + (-90.0, -60.0), # South Pacific sector +] + +# Seed points for Pacific Ocean +const PACIFIC_SEED_POINTS = [ + (180.0, 0.0), # Central equatorial Pacific (dateline) + (-150.0, 20.0), # North Pacific (Hawaii region) + (-120.0, -20.0), # South Pacific + # Same values but from 0 to 360 + (180.0, 0.0), # Central equatorial Pacific + (-150.0 + 360, 20.0), # North Pacific + (-120.0 + 360, -20.0), # South Pacific +] + +##### +##### OceanBasinMask +##### + +add_barrier(v::AbstractVector, b::Barrier) = [v..., b] +add_barrier(v::Barrier, b::Barrier) = [v, b] +add_barrier(::Nothing, b::Barrier) = b + +""" + OceanBasinMask(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) + +Create a mask identifying cells in the Ocean basin using connected component analysis. +The algorithm uses labels restricted by south, north, east, and west boundaries with a fill +algorithm that retrieves the basin corresponding to labels associated with known `seed_points`. +This ensures the mask exactly follows coastlines. + +Arguments +========= +- `grid`: The ocean grid (typically an ImmersedBoundaryGrid) + +Keyword Arguments +================= +- `south_boundary`: Southern latitude limit. Default: nothing +- `north_boundary`: Northern latitude limit. Default: nothing +- `seed_points`: Known (λ, φ) points of the ocean basin to retrieve. Default: [(0, 0)] +- `barriers`: Collection of barriers to apply before labeling. Barriers temporarily + separate connected ocean basins (e.g., separating Atlantic from Pacific). +""" +function OceanBasinMask(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) + + # The computations are 2D and require serial algorithms, so + # we perform the computation on the CPU then move the output + # to the GPU if the initial grid was a GPU grid + cpu_grid = Oceananigans.on_architecture(CPU(), grid) + + # Enforce north and south boundaries + if !isnothing(south_boundary) + barriers = add_barrier(barriers, Barrier(nothing, nothing, -90, south_boundary)) + end + + if !isnothing(north_boundary) + barriers = add_barrier(barriers, Barrier(nothing, nothing, north_boundary, 90)) + end + + # Compute connected component labels for all ocean cells + # Barriers are applied to temporarily separate connected basins + labels = label_ocean_basins(cpu_grid; barriers) + + # Find the Basin label using seed points + basin_label = 0 + for (λs, φs) in seed_points + label = find_label_at_point(labels, cpu_grid, λs, φs) + if label > 0 + basin_label = label + break + end + end + + if basin_label == 0 + @warn "Could not find the Ocean basin in grid. Returning empty mask." + mask = Field{Center, Center, Nothing}(grid) + return OceanBasinMask(mask, grid, seed_points) + end + + # Create mask from label with latitude bounds + mask = create_basin_mask_from_label(cpu_grid, labels, basin_label) + mask = Oceananigans.on_architecture(architecture(grid), mask) + + return OceanBasinMask(mask, grid, seed_points) +end + +##### +##### Convenience functions for specific ocean basins +##### + +""" + atlantic_ocean_mask(grid; include_southern_ocean=false, kw...) + +Create a mask for the Atlantic Ocean with predefined barriers and seed points. + +Keyword Arguments +================= +- `include_southern_ocean`: If `true`, extends the Atlantic basin into the Southern Ocean + sector below the standard separation latitude (~55°S). Default: `false`. +- `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) +- `north_boundary`: Northern latitude limit. Default: 65.0 +- Other keyword arguments are passed to `OceanBasinMask`. +""" +function atlantic_ocean_mask(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 65.0, + barriers = ATLANTIC_OCEAN_BARRIERS, + seed_points = ATLANTIC_SEED_POINTS, + kw...) + + if !include_southern_ocean + barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] + end + + return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +end + +""" + indian_ocean_mask(grid; include_southern_ocean=false, kw...) + +Create a mask for the Indian Ocean with predefined barriers and seed points. + +Keyword Arguments +================= +- `include_southern_ocean`: If `true`, extends the Indian basin into the Southern Ocean + sector below the standard separation latitude (~55°S). Default: `false`. +- `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) +- `north_boundary`: Northern latitude limit. Default: 30.0 +- Other keyword arguments are passed to `OceanBasinMask`. +""" +function indian_ocean_mask(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 30.0, + barriers = INDIAN_OCEAN_BARRIERS, + seed_points = INDIAN_SEED_POINTS, + kw...) + + if !include_southern_ocean + barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] + end + + return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +end + +""" + southern_ocean_mask(grid; kw...) + +Create a mask for the Southern Ocean with predefined barriers and seed points. +Default boundaries: south=-90.0, north=-35.0 +""" +function southern_ocean_mask(grid; + south_boundary = -90.0, + north_boundary = -35.0, + barriers = SOUTHERN_OCEAN_BARRIERS, + seed_points = SOUTHERN_SEED_POINTS, + kw...) + return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +end + +""" + pacific_ocean_mask(grid; include_southern_ocean=false, kw...) + +Create a mask for the Pacific Ocean with predefined barriers and seed points. + +Keyword Arguments +================= +- `include_southern_ocean`: If `true`, extends the Pacific basin into the Southern Ocean + sector below the standard separation latitude (~55°S). Default: `false`. +- `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) +- `north_boundary`: Northern latitude limit. Default: 65.0 +- Other keyword arguments are passed to `OceanBasinMask`. +""" +function pacific_ocean_mask(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 65.0, + barriers = PACIFIC_OCEAN_BARRIERS, + seed_points = PACIFIC_SEED_POINTS, + kw...) + + if !include_southern_ocean + barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] + end + + return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +end + +""" + arctic_ocean_mask(grid; kw...) + +Create a mask for the Arctic Ocean with predefined seed points. +Default boundaries: south=65.0, north=91.0 +""" +function arctic_ocean_mask(grid; + include_southern_ocean = true, + south_boundary = 65.0, + north_boundary = 91.0, + barriers = nothing, + seed_points = [(nothing, 90.0)], + kw...) + + return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +end diff --git a/src/Bathymetry/regrid_bathymetry.jl b/src/Bathymetry/regrid_bathymetry.jl index 430cad426..6227b78f0 100644 --- a/src/Bathymetry/regrid_bathymetry.jl +++ b/src/Bathymetry/regrid_bathymetry.jl @@ -413,49 +413,6 @@ Arguments """ function remove_minor_basins!(zb::Field, keep_major_basins) - zb_cpu = on_architecture(CPU(), zb) - TX = topology(zb_cpu.grid, 1) - - core_size = Nx, Ny, _ = size(zb_cpu.grid) - zb_data = maybe_extend_longitude(zb_cpu, TX()) # Outputs a 2D AbstractArray - - # For periodic grids, `zb_data` may include an extended longitude domain so that - # ImageMorphology can identify connected components correctly across periodic - # boundaries. However, we want `keep_major_basins` to refer to the number of - # major basins in the *physical* model domain only. - # - # To achieve this, we rank basins by their area within the core region - # `1:Nx, 1:Ny` of `zb_data`, rather than within the whole extended array. - remove_minor_basins!(zb_data, keep_major_basins, core_size) - set!(zb, zb_data[1:Nx, 1:Ny]) - - return zb -end - -maybe_extend_longitude(zb_cpu, tx) = interior(zb_cpu, :, :, 1) - -# Since the strel algorithm in `remove_major_basins` does not recognize periodic boundaries, -# before removing connected regions, we extend the longitude direction if it is periodic. -# An extension of half the domain is enough. -function maybe_extend_longitude(zb_cpu, ::Periodic) - Nx = size(zb_cpu, 1) - nx = Nx ÷ 2 - - zb_data = zb_cpu.data[1:Nx, :, 1] - zb_parent = zb_data.parent - - # Add information on the LHS and to the RHS - zb_parent = vcat(zb_parent[nx:Nx, :], zb_parent, zb_parent[1:nx, :]) - - # Update offsets - yoffsets = zb_cpu.data.offsets[2] - xoffsets = - nx - - return OffsetArray(zb_parent, xoffsets, yoffsets) -end - -function remove_minor_basins!(zb, keep_major_basins, core_size) - if !isfinite(keep_major_basins) throw(ArgumentError("`keep_major_basins` must be a finite number!")) end @@ -464,56 +421,60 @@ function remove_minor_basins!(zb, keep_major_basins, core_size) throw(ArgumentError("keep_major_basins must be larger than 0.")) end - water = zb .< 0 + zb_cpu = on_architecture(CPU(), zb) + TX, TY, _ = topology(zb_cpu.grid) - connectivity = ImageMorphology.strel(water) - labels = ImageMorphology.label_components(connectivity) + Nx = Base.length(Center(), TX(), zb_cpu.grid.Nx) + Ny = Base.length(Center(), TY(), zb_cpu.grid.Ny) + # Get labels for the core region (extension is handled internally by label_ocean_basins) + labels = label_ocean_basins(zb_cpu, TX, (Nx, Ny)) nlabels = maximum(labels) - # Rank labels by the number of elements they occupy within the *core* region. - - Nx, Ny = core_size[1], core_size[2] - core_view = labels[1:Nx, 1:Ny] + if nlabels == 0 + return zb # No basins found + end + # Rank labels by the number of elements they occupy total_elements = zeros(nlabels) label_elements = zeros(Int, nlabels) for e in 1:nlabels - cnt = count(==(e), core_view) + cnt = count(==(e), labels) total_elements[e] = cnt label_elements[e] = e end - # Ignore labels that do not intersect the core at all. - # These correspond to basins that live only in the periodic extension. + # Find valid basins (those with at least one cell) valid = findall(>(0), total_elements) - mm_basins = [] # major basins indexes + major_basins = Int[] # indices of major basins to keep m = 1 - # We add basin indexes until we reach the specified number (m == keep_major_basins) or - # we run out of basins to keep -> isempty(valid) + # Keep the largest basins up to keep_major_basins while (m <= keep_major_basins) && !isempty(valid) - # Among the remaining valid labels, find the one with the largest core area. _, idx = findmax(total_elements[valid]) next_label = label_elements[valid[idx]] - push!(mm_basins, next_label) + push!(major_basins, next_label) deleteat!(valid, idx) m += 1 end - labels = map(Float64, labels) + # Modify the bathymetry: set minor basin cells to 0 (land) + # Work on interior view which directly modifies the underlying data + zb_data = interior(zb_cpu, :, :, 1) - for ℓ = 1:maximum(labels) - remove_basin = all(ℓ != m for m in mm_basins) - if remove_basin - labels[labels .== ℓ] .= NaN # Regions to remove + for j in 1:Ny, i in 1:Nx + label = labels[i, j] + if label > 0 && !(label in major_basins) + zb_data[i, j] = 0 # Flatten this cell (make it land) end end - # Flatten minor basins, corresponding to regions where `labels == NaN` - zb[isnan.(labels)] .= 0 + # If original field was on a different architecture, copy back + if zb !== zb_cpu + set!(zb, zb_cpu) + end - return nothing + return zb end diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 8b23f354a..8acc0c07c 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -28,6 +28,14 @@ export JRA55PrescribedAtmosphere, JRA55NetCDFBackend, regrid_bathymetry, + label_ocean_basins, + OceanBasinMask, + Barrier, + atlantic_ocean_mask, + indian_ocean_mask, + southern_ocean_mask, + pacific_ocean_mask, + arctic_ocean_mask, Metadata, Metadatum, ECCOMetadatum, diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 73faf02c2..803adca22 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -5,7 +5,13 @@ using NumericalEarth.Bathymetry: remove_minor_basins!, BathymetryRegridding, cache_filename, load_bathymetry_cache, - save_bathymetry_cache + save_bathymetry_cache, + label_ocean_basins, + find_label_at_point, + OceanBasinMask, + atlantic_ocean_mask, + Barrier, + ATLANTIC_OCEAN_BARRIERS using NumericalEarth.DataWrangling.ETOPO using Statistics @@ -140,3 +146,108 @@ end result4 = regrid_bathymetry(grid; cache=false) @test parent(result1) == parent(result4) end + +@testset "Barrier geometry" begin + @info "Testing barrier geometry utilities..." + + # Test Barrier construction with explicit bounds + barrier = Barrier(-10.0, 10.0, -5.0, 5.0) + @test barrier.west == -10.0 + @test barrier.east == 10.0 + @test barrier.south == -5.0 + @test barrier.north == 5.0 + + # Test Barrier with keyword arguments + barrier_kw = Barrier(west=-10.0, east=10.0, south=-5.0, north=5.0) + @test barrier_kw.west == -10.0 + @test barrier_kw.east == 10.0 + + # Test meridional barrier constructor (3 args + width) + meridional = Barrier(20.0, -36.0, -30.0) # longitude, south, north + @test meridional.west == 19.0 # 20 - 2/2 + @test meridional.east == 21.0 # 20 + 2/2 + @test meridional.south == -36.0 + @test meridional.north == -30.0 + + # Test meridional barrier with custom width + meridional_wide = Barrier(20.0, -36.0, -30.0; width=4.0) + @test meridional_wide.west == 18.0 + @test meridional_wide.east == 22.0 + + # Test LatitudeBand + band = Barrier(-180.0, 180.0, -60.0, -55.0) + @test band.west == -180.0 + @test band.east == 180.0 + @test band.south == -60.0 + @test band.north == -55.0 +end + +@testset "Ocean basin labeling with barriers" begin + @info "Testing ocean basin labeling with barriers..." + + for arch in test_architectures + # Create a global grid + grid = LatitudeLongitudeGrid(arch; + size = (90, 45, 10), + longitude = (-180, 180), + latitude = (-90, 90), + z = (-6000, 0)) + + # Regrid real bathymetry onto this grid + bottom_height = regrid_bathymetry(grid) + ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) + + # Test labeling without barriers - all oceans should have the same label + # (they're connected via the Southern Ocean) + labels_no_barrier = label_ocean_basins(ibg) + + # Find labels at Atlantic and Pacific seed points + atlantic_label = find_label_at_point(labels_no_barrier, ibg, -30.0, 0.0) + pacific_label = find_label_at_point(labels_no_barrier, ibg, -170.0, 0.0) + + # Without barriers, Atlantic and Pacific should have the same label + # (connected via Southern Ocean) + @test atlantic_label == pacific_label + @test atlantic_label > 0 # Should find a valid basin + + # Test labeling with barriers - oceans should be separated + labels_with_barrier = label_ocean_basins(ibg; barriers=ATLANTIC_OCEAN_BARRIERS) + + # With barriers, Atlantic should still be found + atlantic_label_with_barrier = find_label_at_point(labels_with_barrier, ibg, -30.0, 0.0) + @test atlantic_label_with_barrier > 0 + end +end + +@testset "OceanBasinMask creation" begin + @info "Testing OceanBasinMask creation..." + + for arch in test_architectures + # Create a global grid at 1° resolution (needed to properly resolve + # Central America and separate Atlantic from Pacific) + grid = LatitudeLongitudeGrid(arch; + size = (360, 180, 10), + longitude = (-180, 180), + latitude = (-90, 90), + z = (-6000, 0)) + + bottom_height = regrid_bathymetry(grid) + ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) + + # Test atlantic_ocean_mask creation + atlantic = atlantic_ocean_mask(ibg) + @test atlantic isa OceanBasinMask + @test sum(atlantic.mask) > 0 # Should have some ocean cells + + mask = on_architecture(CPU(), atlantic.mask) + + # Test that the mask is properly bounded + # Atlantic mask should not include cells in the Pacific + # (seed point at -170°, 0° should be 0) + pacific_point_i = findfirst(i -> -175 < i < -165, range(-180, 180, length=360)) + equator_j = 90 # equator for 180 latitude points + if !isnothing(pacific_point_i) + @test mask[pacific_point_i, equator_j, 1] == 0 + end + end +end From 23fa1bbe3f5b8ba49aa6b911fe4a0ed25b091793 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 24 Mar 2026 11:24:50 +0100 Subject: [PATCH 02/11] Update to current ClimaOcean preparing for CliMA/ClimaOcean.jl#757 Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 2 +- src/Oceans/ocean_simulation.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 211dec35f..3eeb1d304 100644 --- a/Project.toml +++ b/Project.toml @@ -79,7 +79,7 @@ SeawaterPolynomials = "0.3.5" SpeedyWeather = "0.17.4" StaticArrays = "1" Statistics = "<0.0.1, 1" -Thermodynamics = "0.15.3" +Thermodynamics = "0.15.3, 1.0" WorldOceanAtlasTools = "0.6" XESMF = "0.1.6" ZipFile = "0.10" diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl index 88627bb1b..a87c6115c 100644 --- a/src/Oceans/ocean_simulation.jl +++ b/src/Oceans/ocean_simulation.jl @@ -183,6 +183,7 @@ defaults on a per-field basis. function ocean_simulation(grid; Δt = estimate_maximum_Δt(grid), closure = default_ocean_closure(), + clock = Clock(grid), tracers = (:T, :S), free_surface = default_free_surface(grid), reference_density = 1020, @@ -303,6 +304,7 @@ function ocean_simulation(grid; end ocean_model = HydrostaticFreeSurfaceModel(grid; + clock, buoyancy, closure, biogeochemistry, From 9218be2a5091a7b68ecdbef6b4732cdc80391617 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 26 Mar 2026 09:24:51 +0100 Subject: [PATCH 03/11] address the review points --- .gitignore | 1 + src/Bathymetry/Bathymetry.jl | 2 +- src/Bathymetry/label_ocean_basins.jl | 40 ++--- src/Bathymetry/ocean_basin_mask.jl | 143 ++++++++++-------- .../sea_ice_ocean_fluxes.jl | 4 +- src/NumericalEarth.jl | 2 +- test/test_bathymetry.jl | 8 +- 7 files changed, 110 insertions(+), 90 deletions(-) diff --git a/.gitignore b/.gitignore index edb030721..bf2ce4e33 100644 --- a/.gitignore +++ b/.gitignore @@ -63,4 +63,5 @@ docs/src/literated/ .CondaPkg CondaPkg.toml +docs/plans/* !docs/CondaPkg.toml diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index 1aad80397..ba0379f35 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -1,7 +1,7 @@ module Bathymetry export regrid_bathymetry, ORCAGrid -export OceanBasinMask +export BasinMask export atlantic_ocean_mask, indian_ocean_mask, southern_ocean_mask, pacific_ocean_mask, arctic_ocean_mask export label_ocean_basins export Barrier diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl index 1416bbce8..e311c8b33 100644 --- a/src/Bathymetry/label_ocean_basins.jl +++ b/src/Bathymetry/label_ocean_basins.jl @@ -2,14 +2,20 @@ using Oceananigans.OrthogonalSphericalShellGrids: TripolarGridOfSomeKind using Oceananigans.Fields: convert_to_0_360 ##### -##### Barrier type for separating ocean basins +##### Barrier type for separating connected regions ##### """ Barrier{W, E, S, N} -A rectangular geographic region used to separate ocean basins during labeling. -All cells within this region are temporarily marked as land. +A rectangular geographic region used to separate connected water regions during labeling. +All cells within this region are temporarily marked as land, preventing +the flood-fill from crossing that area. + +Note: `Barrier` differs from `BoundingBox` (in `DataWrangling`) in purpose: +a `BoundingBox` selects a spatial subset of a dataset for downloading/regridding, +whereas a `Barrier` temporarily blocks water cells so that the connected-component +labeling algorithm treats the two sides as separate basins. Fields ====== @@ -80,28 +86,28 @@ end @inbounds zb[i, j, 1] = ifelse(in_lon & in_lat, zero(grid), zb[i, j, 1]) end -# Since the strel algorithm in `remove_major_basins` does not recognize periodic boundaries, -# before removing connected regions, we extend the longitude direction if it is periodic. -# An extension of half the domain is enough. -function maybe_extend_longitude(zb_cpu::Field, ::Periodic) - Nx = size(zb_cpu, 1) - nx = Nx ÷ 2 +# Since the strel algorithm in `label_components` does not recognize periodic boundaries, +# before labeling connected regions we copy half the longitude extent on each side so that +# water cells near the boundary are correctly identified as connected. +function copy_periodic_longitude(zb_cpu::Field, ::Periodic) + Nλ = size(zb_cpu, 1) + half = Nλ ÷ 2 - zb_data = zb_cpu.data[1:Nx, :, 1] + zb_data = zb_cpu.data[1:Nλ, :, 1] zb_parent = zb_data.parent - # Add information on the LHS and to the RHS - zb_parent = vcat(zb_parent[nx:Nx, :], zb_parent, zb_parent[1:nx, :]) + # Concatenate a copy of the eastern half on the left and the western half on the right. + # This is an O(Nλ × Nφ) CPU allocation, acceptable for the serial labeling step. + zb_parent = vcat(zb_parent[half:Nλ, :], zb_parent, zb_parent[1:half, :]) - # Update offsets: the original domain starts at parent index (Nx - nx + 2) - # For offset index 1 to map there, we need offset = -(nx + 1) + # Update offsets so that index 1 maps back to the original domain start yoffsets = zb_cpu.data.offsets[2] - xoffsets = - nx - 1 + xoffsets = - half - 1 return OffsetArray(zb_parent, xoffsets, yoffsets) end -maybe_extend_longitude(zb_cpu::Field, tx) = interior(zb_cpu, :, :, 1) +copy_periodic_longitude(zb_cpu::Field, tx) = interior(zb_cpu, :, :, 1) """ label_ocean_basins(zb_field, TX, core_size) @@ -112,7 +118,7 @@ computing the masks for oceanic basins. Handles periodic boundary extension internally and returns labels for the core region only. """ function label_ocean_basins(zb_field, TX, size) - zb = maybe_extend_longitude(zb_field, TX()) # Outputs a 2D AbstractArray + zb = copy_periodic_longitude(zb_field, TX()) # Outputs a 2D AbstractArray water = zb .< 0 diff --git a/src/Bathymetry/ocean_basin_mask.jl b/src/Bathymetry/ocean_basin_mask.jl index 451610e31..4b291f80c 100644 --- a/src/Bathymetry/ocean_basin_mask.jl +++ b/src/Bathymetry/ocean_basin_mask.jl @@ -2,54 +2,59 @@ using Oceananigans.Grids: λnode, φnode, on_architecture using ImageMorphology ##### -##### OceanBasinMask struct +##### BasinMask struct ##### """ - OceanBasinMask{M, G} + BasinMask{M, G, S} -A mask identifying cells belonging to a specific ocean basin. -The mask is stored as a 2D Field{Center, Center, Nothing} with values: -- 1.0: cell belongs to the basin -- 0.0: cell does not belong to the basin +A mask identifying cells belonging to a specific basin (connected water region). + +The mask is stored as a 2D `Field{Center, Center, Nothing}` with `Bool` values: +- `true`: cell belongs to the basin +- `false`: cell does not belong to the basin + +`seed_points` is a collection of `(λ, φ)` coordinate pairs used to identify which +connected component corresponds to the desired basin. The algorithm finds the +label at each seed point and builds the mask from that label. """ -struct OceanBasinMask{M, G, S} +struct BasinMask{M, G, S} mask :: M grid :: G seed_points :: S end -Base.summary(obm::OceanBasinMask) = "OceanBasinMask" +Base.summary(bm::BasinMask) = "BasinMask" -function Base.show(io::IO, obm::OceanBasinMask) - print(io, summary(obm), " on ", summary(obm.grid)) +function Base.show(io::IO, bm::BasinMask) + print(io, summary(bm), " on ", summary(bm.grid)) end # Forward getindex to mask -Base.getindex(obm::OceanBasinMask, i, j, k) = obm.mask[i, j, k] +Base.getindex(bm::BasinMask, i, j, k) = bm.mask[i, j, k] ##### ##### Additive basin masks ##### """ - Base.:+(mask1::OceanBasinMask, mask2::OceanBasinMask) + Base.:+(mask1::BasinMask, mask2::BasinMask) -Combine two ocean basin masks, returning a new mask that is the union of both. -Cells that belong to either basin will have value 1.0. +Combine two basin masks, returning a new mask that is the union of both. +Cells that belong to either basin will be `true`. """ -function Base.:+(mask1::OceanBasinMask, mask2::OceanBasinMask) +function Base.:+(mask1::BasinMask, mask2::BasinMask) grid = mask1.grid - combined_mask = Field{Center, Center, Nothing}(grid) + combined_mask = Field{Center, Center, Nothing}(grid, Bool) # Union: cell is in combined mask if it's in either mask - parent(combined_mask) .= max.(parent(mask1.mask), parent(mask2.mask)) + parent(combined_mask) .= parent(mask1.mask) .| parent(mask2.mask) fill_halo_regions!(combined_mask) # Combine seed points combined_seeds = (mask1.seed_points..., mask2.seed_points...) - return OceanBasinMask(combined_mask, grid, combined_seeds) + return BasinMask(combined_mask, grid, combined_seeds) end ##### @@ -97,11 +102,9 @@ end Create a mask field for all cells with the given label. """ function create_basin_mask_from_label(grid, labels, basin_label) - arch = architecture(grid) Nx, Ny, _ = size(grid) - FT = eltype(grid) - mask = Field{Center, Center, Nothing}(grid) + mask = Field{Center, Center, Nothing}(grid, Bool) launch!(CPU(), grid, :xy, _compute_basin_mask!, mask, grid, labels, basin_label) @@ -176,7 +179,7 @@ const PACIFIC_SEED_POINTS = [ ] ##### -##### OceanBasinMask +##### BasinMask constructor ##### add_barrier(v::AbstractVector, b::Barrier) = [v..., b] @@ -184,34 +187,44 @@ add_barrier(v::Barrier, b::Barrier) = [v, b] add_barrier(::Nothing, b::Barrier) = b """ - OceanBasinMask(grid; - south_boundary = nothing, - north_boundary = nothing, - seed_points = [(0, 0)], - barriers = nothing) + BasinMask(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) + +Create a `Bool` mask for a single connected water basin. -Create a mask identifying cells in the Ocean basin using connected component analysis. -The algorithm uses labels restricted by south, north, east, and west boundaries with a fill -algorithm that retrieves the basin corresponding to labels associated with known `seed_points`. -This ensures the mask exactly follows coastlines. +The algorithm first labels every connected water region (using `ImageMorphology.label_components`), +optionally splitting regions with `barriers`. It then retrieves the basin whose label +matches the first `seed_point` that falls on water. + +Multiple `seed_points` are tried in order as fallbacks (the first hit wins); +multiple `barriers` are applied simultaneously before labeling so that each +barrier independently blocks water connectivity. Arguments ========= -- `grid`: The ocean grid (typically an ImmersedBoundaryGrid) +- `grid`: An `ImmersedBoundaryGrid` whose immersed boundary defines the coastlines. Keyword Arguments ================= -- `south_boundary`: Southern latitude limit. Default: nothing -- `north_boundary`: Northern latitude limit. Default: nothing -- `seed_points`: Known (λ, φ) points of the ocean basin to retrieve. Default: [(0, 0)] -- `barriers`: Collection of barriers to apply before labeling. Barriers temporarily - separate connected ocean basins (e.g., separating Atlantic from Pacific). +- `south_boundary`: Southern latitude limit — cells south of this become land. Default: `nothing`. +- `north_boundary`: Northern latitude limit — cells north of this become land. Default: `nothing`. +- `seed_points`: `(λ, φ)` pairs identifying the target basin. The first seed that lands + on water determines which connected component becomes `true` in the mask. + Multiple seeds are tried as fallbacks for grids where a single point + may fall on land. Default: `[(0, 0)]`. +- `barriers`: `Barrier`s (or a `Vector` of them) applied before labeling. + Each barrier temporarily marks its rectangle as land, preventing the + flood-fill from crossing it (e.g., closing Drake Passage separates the + Atlantic from the Pacific). Default: `nothing`. """ -function OceanBasinMask(grid; - south_boundary = nothing, - north_boundary = nothing, - seed_points = [(0, 0)], - barriers = nothing) +function BasinMask(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) # The computations are 2D and require serial algorithms, so # we perform the computation on the CPU then move the output @@ -242,26 +255,26 @@ function OceanBasinMask(grid; end if basin_label == 0 - @warn "Could not find the Ocean basin in grid. Returning empty mask." - mask = Field{Center, Center, Nothing}(grid) - return OceanBasinMask(mask, grid, seed_points) + @warn "Could not find the basin in grid. Returning empty mask." + mask = Field{Center, Center, Nothing}(grid, Bool) + return BasinMask(mask, grid, seed_points) end # Create mask from label with latitude bounds mask = create_basin_mask_from_label(cpu_grid, labels, basin_label) mask = Oceananigans.on_architecture(architecture(grid), mask) - return OceanBasinMask(mask, grid, seed_points) + return BasinMask(mask, grid, seed_points) end ##### -##### Convenience functions for specific ocean basins +##### Convenience functions for Earth's ocean basins ##### """ atlantic_ocean_mask(grid; include_southern_ocean=false, kw...) -Create a mask for the Atlantic Ocean with predefined barriers and seed points. +Create a mask for Earth's Atlantic Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -269,7 +282,7 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 65.0 -- Other keyword arguments are passed to `OceanBasinMask`. +- Other keyword arguments are passed to `BasinMask`. """ function atlantic_ocean_mask(grid; include_southern_ocean = true, @@ -283,13 +296,13 @@ function atlantic_ocean_mask(grid; barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ indian_ocean_mask(grid; include_southern_ocean=false, kw...) -Create a mask for the Indian Ocean with predefined barriers and seed points. +Create a mask for Earth's Indian Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -297,7 +310,7 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 30.0 -- Other keyword arguments are passed to `OceanBasinMask`. +- Other keyword arguments are passed to `BasinMask`. """ function indian_ocean_mask(grid; include_southern_ocean = true, @@ -311,28 +324,28 @@ function indian_ocean_mask(grid; barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ southern_ocean_mask(grid; kw...) -Create a mask for the Southern Ocean with predefined barriers and seed points. +Create a mask for Earth's Southern Ocean with predefined barriers and seed points. Default boundaries: south=-90.0, north=-35.0 """ function southern_ocean_mask(grid; - south_boundary = -90.0, - north_boundary = -35.0, - barriers = SOUTHERN_OCEAN_BARRIERS, - seed_points = SOUTHERN_SEED_POINTS, - kw...) - return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + south_boundary = -90.0, + north_boundary = -35.0, + barriers = SOUTHERN_OCEAN_BARRIERS, + seed_points = SOUTHERN_SEED_POINTS, + kw...) + return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ pacific_ocean_mask(grid; include_southern_ocean=false, kw...) -Create a mask for the Pacific Ocean with predefined barriers and seed points. +Create a mask for Earth's Pacific Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -340,7 +353,7 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 65.0 -- Other keyword arguments are passed to `OceanBasinMask`. +- Other keyword arguments are passed to `BasinMask`. """ function pacific_ocean_mask(grid; include_southern_ocean = true, @@ -354,13 +367,13 @@ function pacific_ocean_mask(grid; barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ arctic_ocean_mask(grid; kw...) -Create a mask for the Arctic Ocean with predefined seed points. +Create a mask for Earth's Arctic Ocean with predefined seed points. Default boundaries: south=65.0, north=91.0 """ function arctic_ocean_mask(grid; @@ -371,5 +384,5 @@ function arctic_ocean_mask(grid; seed_points = [(nothing, 90.0)], kw...) - return OceanBasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index ddc5f7bbf..161072987 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -198,8 +198,8 @@ end # ============================================= # Returns interfacial heat flux, melt rate qᵐ, and interface T, S 𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_interface_heat_flux(flux_formulation, - ocean_surface_state, ice_state, - liquidus, ocean_properties, ℰ, u★) + ocean_surface_state, ice_state, + liquidus, ocean_properties, ℰ, u★) # Store interface values and heat flux @inbounds T★[i, j, 1] = Tᵦ diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 8acc0c07c..2b1a88b41 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -29,7 +29,7 @@ export JRA55NetCDFBackend, regrid_bathymetry, label_ocean_basins, - OceanBasinMask, + BasinMask, Barrier, atlantic_ocean_mask, indian_ocean_mask, diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 803adca22..8244d9d6b 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -8,7 +8,7 @@ using NumericalEarth.Bathymetry: remove_minor_basins!, save_bathymetry_cache, label_ocean_basins, find_label_at_point, - OceanBasinMask, + BasinMask, atlantic_ocean_mask, Barrier, ATLANTIC_OCEAN_BARRIERS @@ -219,8 +219,8 @@ end end end -@testset "OceanBasinMask creation" begin - @info "Testing OceanBasinMask creation..." +@testset "BasinMask creation" begin + @info "Testing BasinMask creation..." for arch in test_architectures # Create a global grid at 1° resolution (needed to properly resolve @@ -236,7 +236,7 @@ end # Test atlantic_ocean_mask creation atlantic = atlantic_ocean_mask(ibg) - @test atlantic isa OceanBasinMask + @test atlantic isa BasinMask @test sum(atlantic.mask) > 0 # Should have some ocean cells mask = on_architecture(CPU(), atlantic.mask) From 6ae8282794081e4579d6b45467fd9290001c1233 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 1 Apr 2026 16:23:49 +0200 Subject: [PATCH 04/11] Update src/Bathymetry/label_ocean_basins.jl Co-authored-by: Gregory L. Wagner --- src/Bathymetry/label_ocean_basins.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl index e311c8b33..6fdb29510 100644 --- a/src/Bathymetry/label_ocean_basins.jl +++ b/src/Bathymetry/label_ocean_basins.jl @@ -91,13 +91,12 @@ end # water cells near the boundary are correctly identified as connected. function copy_periodic_longitude(zb_cpu::Field, ::Periodic) Nλ = size(zb_cpu, 1) - half = Nλ ÷ 2 - zb_data = zb_cpu.data[1:Nλ, :, 1] zb_parent = zb_data.parent # Concatenate a copy of the eastern half on the left and the western half on the right. # This is an O(Nλ × Nφ) CPU allocation, acceptable for the serial labeling step. + half = Nλ ÷ 2 zb_parent = vcat(zb_parent[half:Nλ, :], zb_parent, zb_parent[1:half, :]) # Update offsets so that index 1 maps back to the original domain start From 99d7a473c88402ea1018146896ea120cd31101b0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 15:54:08 +0200 Subject: [PATCH 05/11] start changing stuff --- src/Bathymetry/Bathymetry.jl | 2 +- src/Bathymetry/label_ocean_basins.jl | 82 ++++++++++------------------ src/Bathymetry/ocean_basin_mask.jl | 37 +++++++------ src/NumericalEarth.jl | 2 +- test/test_bathymetry.jl | 52 +++++++----------- 5 files changed, 70 insertions(+), 105 deletions(-) diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index ba0379f35..114241ddc 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -4,7 +4,7 @@ export regrid_bathymetry, ORCAGrid export BasinMask export atlantic_ocean_mask, indian_ocean_mask, southern_ocean_mask, pacific_ocean_mask, arctic_ocean_mask export label_ocean_basins -export Barrier +export meridional_barrier using Downloads using ImageMorphology diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl index 6fdb29510..a930bb9d7 100644 --- a/src/Bathymetry/label_ocean_basins.jl +++ b/src/Bathymetry/label_ocean_basins.jl @@ -1,62 +1,31 @@ using Oceananigans.OrthogonalSphericalShellGrids: TripolarGridOfSomeKind using Oceananigans.Fields: convert_to_0_360 +using ..DataWrangling: BoundingBox ##### -##### Barrier type for separating connected regions +##### Barrier application functions ##### -""" - Barrier{W, E, S, N} - -A rectangular geographic region used to separate connected water regions during labeling. -All cells within this region are temporarily marked as land, preventing -the flood-fill from crossing that area. - -Note: `Barrier` differs from `BoundingBox` (in `DataWrangling`) in purpose: -a `BoundingBox` selects a spatial subset of a dataset for downloading/regridding, -whereas a `Barrier` temporarily blocks water cells so that the connected-component -labeling algorithm treats the two sides as separate basins. - -Fields -====== -- `west`: Western longitude limit (degrees) -- `east`: Eastern longitude limit (degrees) -- `south`: Southern latitude limit (degrees) -- `north`: Northern latitude limit (degrees) - -Constructors -============ -- `Barrier(west, east, south, north)`: Create a barrier with explicit bounds -- `Barrier(; west, east, south, north)`: Keyword argument version -""" -struct Barrier{W, E, S, N} - west :: W - east :: E - south :: S - north :: N -end - -Barrier(; west, east, south, north) = Barrier(west, east, south, north) +# A barrier is a `BoundingBox` whose horizontal extent defines a rectangular +# region to be temporarily marked as land during connected-component labeling. +# The `z` field of the `BoundingBox` is ignored: the flood-fill operates on +# the 2D bathymetry slice, and barriers act at every depth. """ - Barrier(longitude, south, north; width=2.0) + meridional_barrier(longitude, south, north; width=2.0) -Create a meridional (north-south) barrier at a given longitude. -Useful for closing straits like Cape Agulhas or Indonesian passages. +Create a narrow meridional `BoundingBox` centered at `longitude` with meridional +extent `[south, north]` and zonal width `width` degrees. Useful for closing +straits like Cape Agulhas or Indonesian passages during basin labeling. """ -Barrier(longitude, south, north; width=2.0) = Barrier(longitude - width/2, longitude + width/2, south, north) - -##### -##### Barrier application functions -##### +meridional_barrier(longitude, south, north; width=2.0) = BoundingBox(longitude=(longitude - width/2, longitude + width/2), latitude=(south, north)) """ - apply_barrier!(zb_data, grid, barrier::Barrier) + apply_barrier!(zb, grid, barrier::BoundingBox) -Mark all cells within the barrier region as land (z = 0). +Mark all cells within the barrier's horizontal region as land (z = 0). """ -apply_barrier!(zb, grid, barrier::Barrier) = - launch!(architecture(grid), grid, :xy, _apply_barrier!, zb, grid, barrier) +apply_barrier!(zb, grid, barrier::BoundingBox) = launch!(architecture(grid), grid, :xy, _apply_barrier!, zb, grid, barrier) apply_barrier!(zb, grid, barriers::Nothing) = zb @@ -67,21 +36,25 @@ function apply_barrier!(zb, grid, barriers::AbstractVector) return zb end -@kernel function _apply_barrier!(zb, grid, barrier::Barrier) +@kernel function _apply_barrier!(zb, grid, barrier::BoundingBox) i, j = @index(Global, NTuple) - in_lon = if isnothing(barrier.west) || (barrier.east - barrier.west >= 360) + in_lon = if isnothing(barrier.longitude) || (barrier.longitude[2] - barrier.longitude[1] >= 360) true else - bw = convert_to_0_360(barrier.west) - be = convert_to_0_360(barrier.east) + bw = convert_to_0_360(barrier.longitude[1]) + be = convert_to_0_360(barrier.longitude[2]) λ = λnode(i, j, 1, grid, Center(), Center(), Center()) λ = convert_to_0_360(λ) (bw <= λ <= be) end - φ = φnode(i, j, 1, grid, Center(), Center(), Center()) - in_lat = barrier.south <= φ <= barrier.north + in_lat = if isnothing(barrier.latitude) + true + else + φ = φnode(i, j, 1, grid, Center(), Center(), Center()) + barrier.latitude[1] <= φ <= barrier.latitude[2] + end @inbounds zb[i, j, 1] = ifelse(in_lon & in_lat, zero(grid), zb[i, j, 1]) end @@ -205,9 +178,10 @@ Label connected ocean basins in an ImmersedBoundaryGrid. Keyword Arguments ================= -- `barriers`: Collection of barriers to apply before labeling. Barriers temporarily - mark certain cells as land, allowing separation of connected ocean basins - (e.g., separating Atlantic from Pacific via the Southern Ocean). +- `barriers`: Collection of `BoundingBox`es applied before labeling. Each barrier + temporarily marks its horizontal rectangle as land, allowing separation + of connected ocean basins (e.g., separating Atlantic from Pacific via + the Southern Ocean). """ function label_ocean_basins(grid::ImmersedBoundaryGrid; barriers=nothing) diff --git a/src/Bathymetry/ocean_basin_mask.jl b/src/Bathymetry/ocean_basin_mask.jl index 4b291f80c..fbcd518d3 100644 --- a/src/Bathymetry/ocean_basin_mask.jl +++ b/src/Bathymetry/ocean_basin_mask.jl @@ -124,25 +124,25 @@ end ##### Some useful Basin seeds and barriers ##### -const SOUTHERN_OCEAN_SEPARATION_BARRIER = Barrier(-180.0, 180.0, -56.0, -54.0) +const SOUTHERN_OCEAN_SEPARATION_BARRIER = BoundingBox(longitude=(-180.0, 180.0), latitude=(-56.0, -54.0)) const ATLANTIC_OCEAN_BARRIERS = [ - Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) - Barrier(289.0, -90.0, -30.0) # Drake passage (meridional barrier) + meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas + meridional_barrier(289.0, -90.0, -30.0), # Drake Passage ] const INDIAN_OCEAN_BARRIERS = [ - Barrier(141.0, -90.0, -3.0), # Indonesian side (meridional) - Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) - Barrier(105.0, 141.0, -4.0, -3.0), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) + meridional_barrier(141.0, -90.0, -3.0), # Indonesian side + meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas + BoundingBox(longitude=(105.0, 141.0), latitude=(-4.0, -3.0)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) ] const SOUTHERN_OCEAN_BARRIERS = [SOUTHERN_OCEAN_SEPARATION_BARRIER] const PACIFIC_OCEAN_BARRIERS = [ - Barrier(141.0, -90.0, -3.0), # Indonesian side (meridional) - Barrier(20.0, -90.0, -30.0), # Cape Agulhas (meridional barrier) - Barrier(105.0, 141.0, -4.0, -3.0), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) + meridional_barrier(141.0, -90.0, -3.0), # Indonesian side + meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas + BoundingBox(longitude=(105.0, 141.0), latitude=(-4.0, -3.0)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) ] # Seed points for Atlantic Ocean (definitely in the Atlantic) @@ -182,9 +182,9 @@ const PACIFIC_SEED_POINTS = [ ##### BasinMask constructor ##### -add_barrier(v::AbstractVector, b::Barrier) = [v..., b] -add_barrier(v::Barrier, b::Barrier) = [v, b] -add_barrier(::Nothing, b::Barrier) = b +add_barrier(v::AbstractVector, b::BoundingBox) = [v..., b] +add_barrier(v::BoundingBox, b::BoundingBox) = [v, b] +add_barrier(::Nothing, b::BoundingBox) = b """ BasinMask(grid; @@ -215,10 +215,11 @@ Keyword Arguments on water determines which connected component becomes `true` in the mask. Multiple seeds are tried as fallbacks for grids where a single point may fall on land. Default: `[(0, 0)]`. -- `barriers`: `Barrier`s (or a `Vector` of them) applied before labeling. - Each barrier temporarily marks its rectangle as land, preventing the - flood-fill from crossing it (e.g., closing Drake Passage separates the - Atlantic from the Pacific). Default: `nothing`. +- `barriers`: A `BoundingBox` (or a `Vector` of them) applied before labeling. + Each barrier temporarily marks its horizontal rectangle as land, + preventing the flood-fill from crossing it (e.g., closing Drake Passage + separates the Atlantic from the Pacific). The `z` field of the + `BoundingBox` is ignored. Default: `nothing`. """ function BasinMask(grid; south_boundary = nothing, @@ -233,11 +234,11 @@ function BasinMask(grid; # Enforce north and south boundaries if !isnothing(south_boundary) - barriers = add_barrier(barriers, Barrier(nothing, nothing, -90, south_boundary)) + barriers = add_barrier(barriers, BoundingBox(longitude=nothing, latitude=(-90, south_boundary))) end if !isnothing(north_boundary) - barriers = add_barrier(barriers, Barrier(nothing, nothing, north_boundary, 90)) + barriers = add_barrier(barriers, BoundingBox(longitude=nothing, latitude=(north_boundary, 90))) end # Compute connected component labels for all ocean cells diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 2b1a88b41..c052ef3f7 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -30,7 +30,7 @@ export regrid_bathymetry, label_ocean_basins, BasinMask, - Barrier, + meridional_barrier, atlantic_ocean_mask, indian_ocean_mask, southern_ocean_mask, diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 8244d9d6b..cb727c571 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -10,8 +10,9 @@ using NumericalEarth.Bathymetry: remove_minor_basins!, find_label_at_point, BasinMask, atlantic_ocean_mask, - Barrier, + meridional_barrier, ATLANTIC_OCEAN_BARRIERS +using NumericalEarth.DataWrangling: BoundingBox using NumericalEarth.DataWrangling.ETOPO using Statistics @@ -150,36 +151,25 @@ end @testset "Barrier geometry" begin @info "Testing barrier geometry utilities..." - # Test Barrier construction with explicit bounds - barrier = Barrier(-10.0, 10.0, -5.0, 5.0) - @test barrier.west == -10.0 - @test barrier.east == 10.0 - @test barrier.south == -5.0 - @test barrier.north == 5.0 - - # Test Barrier with keyword arguments - barrier_kw = Barrier(west=-10.0, east=10.0, south=-5.0, north=5.0) - @test barrier_kw.west == -10.0 - @test barrier_kw.east == 10.0 - - # Test meridional barrier constructor (3 args + width) - meridional = Barrier(20.0, -36.0, -30.0) # longitude, south, north - @test meridional.west == 19.0 # 20 - 2/2 - @test meridional.east == 21.0 # 20 + 2/2 - @test meridional.south == -36.0 - @test meridional.north == -30.0 - - # Test meridional barrier with custom width - meridional_wide = Barrier(20.0, -36.0, -30.0; width=4.0) - @test meridional_wide.west == 18.0 - @test meridional_wide.east == 22.0 - - # Test LatitudeBand - band = Barrier(-180.0, 180.0, -60.0, -55.0) - @test band.west == -180.0 - @test band.east == 180.0 - @test band.south == -60.0 - @test band.north == -55.0 + # A barrier is a `BoundingBox` whose horizontal extent defines the rectangle + # to be treated as land during connected-component labeling. + barrier = BoundingBox(longitude=(-10.0, 10.0), latitude=(-5.0, 5.0)) + @test barrier.longitude == (-10.0, 10.0) + @test barrier.latitude == (-5.0, 5.0) + + # Test meridional_barrier constructor (longitude, south, north) + meridional = meridional_barrier(20.0, -36.0, -30.0) + @test meridional.longitude == (19.0, 21.0) # 20 ± 2/2 + @test meridional.latitude == (-36.0, -30.0) + + # Test meridional_barrier with custom width + meridional_wide = meridional_barrier(20.0, -36.0, -30.0; width=4.0) + @test meridional_wide.longitude == (18.0, 22.0) + + # Latitude band covering the full longitude range + band = BoundingBox(longitude=(-180.0, 180.0), latitude=(-60.0, -55.0)) + @test band.longitude == (-180.0, 180.0) + @test band.latitude == (-60.0, -55.0) end @testset "Ocean basin labeling with barriers" begin From 5a6f64eaa65cbbb8add4c60eddfcf96ec16a9435 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 16:02:58 +0200 Subject: [PATCH 06/11] BasinMask -> Basin --- src/Bathymetry/Bathymetry.jl | 6 +- .../{ocean_basin_mask.jl => ocean_basin.jl} | 203 +++++++++--------- src/NumericalEarth.jl | 12 +- test/test_bathymetry.jl | 14 +- 4 files changed, 120 insertions(+), 115 deletions(-) rename src/Bathymetry/{ocean_basin_mask.jl => ocean_basin.jl} (65%) diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index 114241ddc..6fbfc45e6 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -1,8 +1,8 @@ module Bathymetry export regrid_bathymetry, ORCAGrid -export BasinMask -export atlantic_ocean_mask, indian_ocean_mask, southern_ocean_mask, pacific_ocean_mask, arctic_ocean_mask +export Basin +export atlantic_ocean_basin, indian_ocean_basin, southern_ocean_basin, pacific_ocean_basin, arctic_ocean_basin export label_ocean_basins export meridional_barrier @@ -28,7 +28,7 @@ using ..DataWrangling.ETOPO: ETOPO2022 include("label_ocean_basins.jl") include("regrid_bathymetry.jl") -include("ocean_basin_mask.jl") +include("ocean_basin.jl") include("orca_grid.jl") end # module diff --git a/src/Bathymetry/ocean_basin_mask.jl b/src/Bathymetry/ocean_basin.jl similarity index 65% rename from src/Bathymetry/ocean_basin_mask.jl rename to src/Bathymetry/ocean_basin.jl index fbcd518d3..7d09bd7d7 100644 --- a/src/Bathymetry/ocean_basin_mask.jl +++ b/src/Bathymetry/ocean_basin.jl @@ -2,59 +2,62 @@ using Oceananigans.Grids: λnode, φnode, on_architecture using ImageMorphology ##### -##### BasinMask struct +##### Basin struct ##### """ - BasinMask{M, G, S} - -A mask identifying cells belonging to a specific basin (connected water region). - -The mask is stored as a 2D `Field{Center, Center, Nothing}` with `Bool` values: -- `true`: cell belongs to the basin -- `false`: cell does not belong to the basin - -`seed_points` is a collection of `(λ, φ)` coordinate pairs used to identify which -connected component corresponds to the desired basin. The algorithm finds the -label at each seed point and builds the mask from that label. + Basin{M, G, S} + +A connected water region (ocean basin) identified on a grid, together with the +boolean mask that labels cells belonging to it and the seed points used to pick +the connected component. + +Fields +====== +- `mask`: a 2D `Field{Center, Center, Nothing}` with `Bool` values — `true` for + cells belonging to the basin, `false` otherwise. +- `grid`: the grid on which the basin is defined. +- `seed_points`: `(λ, φ)` coordinate pairs used to identify which connected + component corresponds to the desired basin. The algorithm finds + the label at each seed point and builds the mask from that label. """ -struct BasinMask{M, G, S} +struct Basin{M, G, S} mask :: M grid :: G seed_points :: S end -Base.summary(bm::BasinMask) = "BasinMask" +Base.summary(basin::Basin) = "Basin" -function Base.show(io::IO, bm::BasinMask) - print(io, summary(bm), " on ", summary(bm.grid)) +function Base.show(io::IO, basin::Basin) + print(io, summary(basin), " on ", summary(basin.grid)) end -# Forward getindex to mask -Base.getindex(bm::BasinMask, i, j, k) = bm.mask[i, j, k] +# Forward getindex to the underlying mask +Base.getindex(basin::Basin, i, j, k) = basin.mask[i, j, k] ##### -##### Additive basin masks +##### Additive basins ##### """ - Base.:+(mask1::BasinMask, mask2::BasinMask) + Base.:+(basin1::Basin, basin2::Basin) -Combine two basin masks, returning a new mask that is the union of both. -Cells that belong to either basin will be `true`. +Combine two basins into a new one whose mask is the union of both. Cells that +belong to either input basin are `true` in the combined mask. """ -function Base.:+(mask1::BasinMask, mask2::BasinMask) - grid = mask1.grid +function Base.:+(basin1::Basin, basin2::Basin) + grid = basin1.grid combined_mask = Field{Center, Center, Nothing}(grid, Bool) - # Union: cell is in combined mask if it's in either mask - parent(combined_mask) .= parent(mask1.mask) .| parent(mask2.mask) + # Union: a cell is in the combined mask if it is in either input mask + parent(combined_mask) .= parent(basin1.mask) .| parent(basin2.mask) fill_halo_regions!(combined_mask) # Combine seed points - combined_seeds = (mask1.seed_points..., mask2.seed_points...) + combined_seeds = (basin1.seed_points..., basin2.seed_points...) - return BasinMask(combined_mask, grid, combined_seeds) + return Basin(combined_mask, grid, combined_seeds) end ##### @@ -179,7 +182,7 @@ const PACIFIC_SEED_POINTS = [ ] ##### -##### BasinMask constructor +##### Basin constructor ##### add_barrier(v::AbstractVector, b::BoundingBox) = [v..., b] @@ -187,17 +190,19 @@ add_barrier(v::BoundingBox, b::BoundingBox) = [v, b] add_barrier(::Nothing, b::BoundingBox) = b """ - BasinMask(grid; - south_boundary = nothing, - north_boundary = nothing, - seed_points = [(0, 0)], - barriers = nothing) + Basin(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) -Create a `Bool` mask for a single connected water basin. +Build a `Basin` — a single connected water region on `grid` together with its +boolean mask. -The algorithm first labels every connected water region (using `ImageMorphology.label_components`), -optionally splitting regions with `barriers`. It then retrieves the basin whose label -matches the first `seed_point` that falls on water. +The algorithm first labels every connected water region (using +`ImageMorphology.label_components`), optionally splitting regions with +`barriers`. It then retrieves the basin whose label matches the first +`seed_point` that falls on water. Multiple `seed_points` are tried in order as fallbacks (the first hit wins); multiple `barriers` are applied simultaneously before labeling so that each @@ -221,11 +226,11 @@ Keyword Arguments separates the Atlantic from the Pacific). The `z` field of the `BoundingBox` is ignored. Default: `nothing`. """ -function BasinMask(grid; - south_boundary = nothing, - north_boundary = nothing, - seed_points = [(0, 0)], - barriers = nothing) +function Basin(grid; + south_boundary = nothing, + north_boundary = nothing, + seed_points = [(0, 0)], + barriers = nothing) # The computations are 2D and require serial algorithms, so # we perform the computation on the CPU then move the output @@ -245,7 +250,7 @@ function BasinMask(grid; # Barriers are applied to temporarily separate connected basins labels = label_ocean_basins(cpu_grid; barriers) - # Find the Basin label using seed points + # Find the basin label using seed points basin_label = 0 for (λs, φs) in seed_points label = find_label_at_point(labels, cpu_grid, λs, φs) @@ -258,14 +263,14 @@ function BasinMask(grid; if basin_label == 0 @warn "Could not find the basin in grid. Returning empty mask." mask = Field{Center, Center, Nothing}(grid, Bool) - return BasinMask(mask, grid, seed_points) + return Basin(mask, grid, seed_points) end # Create mask from label with latitude bounds mask = create_basin_mask_from_label(cpu_grid, labels, basin_label) mask = Oceananigans.on_architecture(architecture(grid), mask) - return BasinMask(mask, grid, seed_points) + return Basin(mask, grid, seed_points) end ##### @@ -273,9 +278,9 @@ end ##### """ - atlantic_ocean_mask(grid; include_southern_ocean=false, kw...) + atlantic_ocean_basin(grid; include_southern_ocean=false, kw...) -Create a mask for Earth's Atlantic Ocean with predefined barriers and seed points. +Build a `Basin` for Earth's Atlantic Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -283,27 +288,27 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 65.0 -- Other keyword arguments are passed to `BasinMask`. +- Other keyword arguments are passed to `Basin`. """ -function atlantic_ocean_mask(grid; - include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 65.0, - barriers = ATLANTIC_OCEAN_BARRIERS, - seed_points = ATLANTIC_SEED_POINTS, - kw...) +function atlantic_ocean_basin(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 65.0, + barriers = ATLANTIC_OCEAN_BARRIERS, + seed_points = ATLANTIC_SEED_POINTS, + kw...) if !include_southern_ocean barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ - indian_ocean_mask(grid; include_southern_ocean=false, kw...) + indian_ocean_basin(grid; include_southern_ocean=false, kw...) -Create a mask for Earth's Indian Ocean with predefined barriers and seed points. +Build a `Basin` for Earth's Indian Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -311,42 +316,42 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 30.0 -- Other keyword arguments are passed to `BasinMask`. +- Other keyword arguments are passed to `Basin`. """ -function indian_ocean_mask(grid; - include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 30.0, - barriers = INDIAN_OCEAN_BARRIERS, - seed_points = INDIAN_SEED_POINTS, - kw...) +function indian_ocean_basin(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 30.0, + barriers = INDIAN_OCEAN_BARRIERS, + seed_points = INDIAN_SEED_POINTS, + kw...) if !include_southern_ocean barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ - southern_ocean_mask(grid; kw...) + southern_ocean_basin(grid; kw...) -Create a mask for Earth's Southern Ocean with predefined barriers and seed points. +Build a `Basin` for Earth's Southern Ocean with predefined barriers and seed points. Default boundaries: south=-90.0, north=-35.0 """ -function southern_ocean_mask(grid; - south_boundary = -90.0, - north_boundary = -35.0, - barriers = SOUTHERN_OCEAN_BARRIERS, - seed_points = SOUTHERN_SEED_POINTS, - kw...) - return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +function southern_ocean_basin(grid; + south_boundary = -90.0, + north_boundary = -35.0, + barriers = SOUTHERN_OCEAN_BARRIERS, + seed_points = SOUTHERN_SEED_POINTS, + kw...) + return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ - pacific_ocean_mask(grid; include_southern_ocean=false, kw...) + pacific_ocean_basin(grid; include_southern_ocean=false, kw...) -Create a mask for Earth's Pacific Ocean with predefined barriers and seed points. +Build a `Basin` for Earth's Pacific Ocean with predefined barriers and seed points. Keyword Arguments ================= @@ -354,36 +359,36 @@ Keyword Arguments sector below the standard separation latitude (~55°S). Default: `false`. - `south_boundary`: Southern latitude limit. Default: -50.0 (or -90.0 if `include_southern_ocean=true`) - `north_boundary`: Northern latitude limit. Default: 65.0 -- Other keyword arguments are passed to `BasinMask`. +- Other keyword arguments are passed to `Basin`. """ -function pacific_ocean_mask(grid; - include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 65.0, - barriers = PACIFIC_OCEAN_BARRIERS, - seed_points = PACIFIC_SEED_POINTS, - kw...) +function pacific_ocean_basin(grid; + include_southern_ocean = true, + south_boundary = include_southern_ocean ? -90.0 : -50.0, + north_boundary = 65.0, + barriers = PACIFIC_OCEAN_BARRIERS, + seed_points = PACIFIC_SEED_POINTS, + kw...) if !include_southern_ocean barriers = [barriers..., SOUTHERN_OCEAN_SEPARATION_BARRIER] end - return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) + return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end """ - arctic_ocean_mask(grid; kw...) + arctic_ocean_basin(grid; kw...) -Create a mask for Earth's Arctic Ocean with predefined seed points. +Build a `Basin` for Earth's Arctic Ocean with predefined seed points. Default boundaries: south=65.0, north=91.0 """ -function arctic_ocean_mask(grid; - include_southern_ocean = true, - south_boundary = 65.0, - north_boundary = 91.0, - barriers = nothing, - seed_points = [(nothing, 90.0)], - kw...) - - return BasinMask(grid; south_boundary, north_boundary, barriers, seed_points, kw...) +function arctic_ocean_basin(grid; + include_southern_ocean = true, + south_boundary = 65.0, + north_boundary = 91.0, + barriers = nothing, + seed_points = [(nothing, 90.0)], + kw...) + + return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) end diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index c052ef3f7..8d9d8cb85 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -29,13 +29,13 @@ export JRA55NetCDFBackend, regrid_bathymetry, label_ocean_basins, - BasinMask, + Basin, meridional_barrier, - atlantic_ocean_mask, - indian_ocean_mask, - southern_ocean_mask, - pacific_ocean_mask, - arctic_ocean_mask, + atlantic_ocean_basin, + indian_ocean_basin, + southern_ocean_basin, + pacific_ocean_basin, + arctic_ocean_basin, Metadata, Metadatum, ECCOMetadatum, diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index cb727c571..3795bdcba 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -8,8 +8,8 @@ using NumericalEarth.Bathymetry: remove_minor_basins!, save_bathymetry_cache, label_ocean_basins, find_label_at_point, - BasinMask, - atlantic_ocean_mask, + Basin, + atlantic_ocean_basin, meridional_barrier, ATLANTIC_OCEAN_BARRIERS using NumericalEarth.DataWrangling: BoundingBox @@ -209,8 +209,8 @@ end end end -@testset "BasinMask creation" begin - @info "Testing BasinMask creation..." +@testset "Basin creation" begin + @info "Testing Basin creation..." for arch in test_architectures # Create a global grid at 1° resolution (needed to properly resolve @@ -224,9 +224,9 @@ end bottom_height = regrid_bathymetry(grid) ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) - # Test atlantic_ocean_mask creation - atlantic = atlantic_ocean_mask(ibg) - @test atlantic isa BasinMask + # Test atlantic_ocean_basin creation + atlantic = atlantic_ocean_basin(ibg) + @test atlantic isa Basin @test sum(atlantic.mask) > 0 # Should have some ocean cells mask = on_architecture(CPU(), atlantic.mask) From a71782c2572a5b291b2058f8f44e2f559d6123ec Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 16:08:04 +0200 Subject: [PATCH 07/11] use an or --- src/Bathymetry/ocean_basin.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Bathymetry/ocean_basin.jl b/src/Bathymetry/ocean_basin.jl index 7d09bd7d7..17e8392c8 100644 --- a/src/Bathymetry/ocean_basin.jl +++ b/src/Bathymetry/ocean_basin.jl @@ -37,16 +37,17 @@ end Base.getindex(basin::Basin, i, j, k) = basin.mask[i, j, k] ##### -##### Additive basins +##### Basin union ##### """ - Base.:+(basin1::Basin, basin2::Basin) + Base.:|(basin1::Basin, basin2::Basin) Combine two basins into a new one whose mask is the union of both. Cells that -belong to either input basin are `true` in the combined mask. +belong to either input basin are `true` in the combined mask. This mirrors the +elementwise `.|` applied to the underlying boolean masks. """ -function Base.:+(basin1::Basin, basin2::Basin) +function Base.:|(basin1::Basin, basin2::Basin) grid = basin1.grid combined_mask = Field{Center, Center, Nothing}(grid, Bool) From cebf2550879381fd76f3cf9f68c42c64c9fa6641 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 16:29:49 +0200 Subject: [PATCH 08/11] floats to integers --- src/Bathymetry/label_ocean_basins.jl | 4 +- src/Bathymetry/ocean_basin.jl | 74 ++++++++++++++-------------- 2 files changed, 39 insertions(+), 39 deletions(-) diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl index a930bb9d7..eb3144fc5 100644 --- a/src/Bathymetry/label_ocean_basins.jl +++ b/src/Bathymetry/label_ocean_basins.jl @@ -12,13 +12,13 @@ using ..DataWrangling: BoundingBox # the 2D bathymetry slice, and barriers act at every depth. """ - meridional_barrier(longitude, south, north; width=2.0) + meridional_barrier(longitude, south, north; width=2) Create a narrow meridional `BoundingBox` centered at `longitude` with meridional extent `[south, north]` and zonal width `width` degrees. Useful for closing straits like Cape Agulhas or Indonesian passages during basin labeling. """ -meridional_barrier(longitude, south, north; width=2.0) = BoundingBox(longitude=(longitude - width/2, longitude + width/2), latitude=(south, north)) +meridional_barrier(longitude, south, north; width=2) = BoundingBox(longitude=(longitude - width/2, longitude + width/2), latitude=(south, north)) """ apply_barrier!(zb, grid, barrier::BoundingBox) diff --git a/src/Bathymetry/ocean_basin.jl b/src/Bathymetry/ocean_basin.jl index 17e8392c8..e4c9eca74 100644 --- a/src/Bathymetry/ocean_basin.jl +++ b/src/Bathymetry/ocean_basin.jl @@ -128,58 +128,58 @@ end ##### Some useful Basin seeds and barriers ##### -const SOUTHERN_OCEAN_SEPARATION_BARRIER = BoundingBox(longitude=(-180.0, 180.0), latitude=(-56.0, -54.0)) +const SOUTHERN_OCEAN_SEPARATION_BARRIER = BoundingBox(longitude=(-180, 180), latitude=(-56, -54)) const ATLANTIC_OCEAN_BARRIERS = [ - meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas - meridional_barrier(289.0, -90.0, -30.0), # Drake Passage + meridional_barrier(20, -90, -30), # Cape Agulhas + meridional_barrier(289, -90, -30), # Drake Passage ] const INDIAN_OCEAN_BARRIERS = [ - meridional_barrier(141.0, -90.0, -3.0), # Indonesian side - meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas - BoundingBox(longitude=(105.0, 141.0), latitude=(-4.0, -3.0)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) + meridional_barrier(141, -90, -3), # Indonesian side + meridional_barrier(20, -90, -30), # Cape Agulhas + BoundingBox(longitude=(105, 141), latitude=(-4, -3)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) ] const SOUTHERN_OCEAN_BARRIERS = [SOUTHERN_OCEAN_SEPARATION_BARRIER] const PACIFIC_OCEAN_BARRIERS = [ - meridional_barrier(141.0, -90.0, -3.0), # Indonesian side - meridional_barrier(20.0, -90.0, -30.0), # Cape Agulhas - BoundingBox(longitude=(105.0, 141.0), latitude=(-4.0, -3.0)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) + meridional_barrier(141, -90, -3), # Indonesian side + meridional_barrier(20, -90, -30), # Cape Agulhas + BoundingBox(longitude=(105, 141), latitude=(-4, -3)), # Indonesian/Asian seas (zonal barrier at 3.5ᵒ S) ] # Seed points for Atlantic Ocean (definitely in the Atlantic) const ATLANTIC_SEED_POINTS = [ - (-30.0, 0.0), # Central equatorial Atlantic - (-40.0, 30.0), # North Atlantic - (-25.0, -20.0), # South Atlantic + (-30, 0), # Central equatorial Atlantic + (-40, 30), # North Atlantic + (-25, -20), # South Atlantic ] # Seed points for Indian Ocean const INDIAN_SEED_POINTS = [ - (70.0, -10.0), # Central Indian Ocean - (60.0, 10.0), # Arabian Sea region - (90.0, -20.0), # Eastern Indian Ocean + (70, -10), # Central Indian Ocean + (60, 10), # Arabian Sea region + (90, -20), # Eastern Indian Ocean ] # Seed points for Southern Ocean const SOUTHERN_SEED_POINTS = [ - (0.0, -60.0), # South Atlantic sector - (90.0, -60.0), # Indian Ocean sector - (180.0, -60.0), # Pacific sector (date line) - (-90.0, -60.0), # South Pacific sector + (0, -60), # South Atlantic sector + (90, -60), # Indian Ocean sector + (180, -60), # Pacific sector (date line) + (-90, -60), # South Pacific sector ] # Seed points for Pacific Ocean const PACIFIC_SEED_POINTS = [ - (180.0, 0.0), # Central equatorial Pacific (dateline) - (-150.0, 20.0), # North Pacific (Hawaii region) - (-120.0, -20.0), # South Pacific - # Same values but from 0 to 360 - (180.0, 0.0), # Central equatorial Pacific - (-150.0 + 360, 20.0), # North Pacific - (-120.0 + 360, -20.0), # South Pacific + (180, 0), # Central equatorial Pacific (dateline) + (-150, 20), # North Pacific (Hawaii region) + (-120, -20), # South Pacific + # Same values but mapped to [0, 360) + (180, 0), # Central equatorial Pacific + (-150 + 360, 20), # North Pacific + (-120 + 360, -20), # South Pacific ] ##### @@ -293,8 +293,8 @@ Keyword Arguments """ function atlantic_ocean_basin(grid; include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 65.0, + south_boundary = include_southern_ocean ? -90 : -50, + north_boundary = 65, barriers = ATLANTIC_OCEAN_BARRIERS, seed_points = ATLANTIC_SEED_POINTS, kw...) @@ -321,8 +321,8 @@ Keyword Arguments """ function indian_ocean_basin(grid; include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 30.0, + south_boundary = include_southern_ocean ? -90 : -50, + north_boundary = 30, barriers = INDIAN_OCEAN_BARRIERS, seed_points = INDIAN_SEED_POINTS, kw...) @@ -341,8 +341,8 @@ Build a `Basin` for Earth's Southern Ocean with predefined barriers and seed poi Default boundaries: south=-90.0, north=-35.0 """ function southern_ocean_basin(grid; - south_boundary = -90.0, - north_boundary = -35.0, + south_boundary = -90, + north_boundary = -35, barriers = SOUTHERN_OCEAN_BARRIERS, seed_points = SOUTHERN_SEED_POINTS, kw...) @@ -364,8 +364,8 @@ Keyword Arguments """ function pacific_ocean_basin(grid; include_southern_ocean = true, - south_boundary = include_southern_ocean ? -90.0 : -50.0, - north_boundary = 65.0, + south_boundary = include_southern_ocean ? -90 : -50, + north_boundary = 65, barriers = PACIFIC_OCEAN_BARRIERS, seed_points = PACIFIC_SEED_POINTS, kw...) @@ -385,10 +385,10 @@ Default boundaries: south=65.0, north=91.0 """ function arctic_ocean_basin(grid; include_southern_ocean = true, - south_boundary = 65.0, - north_boundary = 91.0, + south_boundary = 65, + north_boundary = 91, barriers = nothing, - seed_points = [(nothing, 90.0)], + seed_points = [(nothing, 90)], kw...) return Basin(grid; south_boundary, north_boundary, barriers, seed_points, kw...) From 71d149a8da736be13e08722d269d8af21bce4a30 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 16:33:29 +0200 Subject: [PATCH 09/11] test compared to integers --- test/test_bathymetry.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 3795bdcba..57022def8 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -153,23 +153,23 @@ end # A barrier is a `BoundingBox` whose horizontal extent defines the rectangle # to be treated as land during connected-component labeling. - barrier = BoundingBox(longitude=(-10.0, 10.0), latitude=(-5.0, 5.0)) - @test barrier.longitude == (-10.0, 10.0) - @test barrier.latitude == (-5.0, 5.0) + barrier = BoundingBox(longitude=(-10, 10), latitude=(-5, 5)) + @test barrier.longitude == (-10, 10) + @test barrier.latitude == (-5, 5) # Test meridional_barrier constructor (longitude, south, north) - meridional = meridional_barrier(20.0, -36.0, -30.0) - @test meridional.longitude == (19.0, 21.0) # 20 ± 2/2 - @test meridional.latitude == (-36.0, -30.0) + meridional = meridional_barrier(20, -36, -30) + @test meridional.longitude == (19, 21) # 20 ± 2/2 + @test meridional.latitude == (-36, -30) # Test meridional_barrier with custom width - meridional_wide = meridional_barrier(20.0, -36.0, -30.0; width=4.0) - @test meridional_wide.longitude == (18.0, 22.0) + meridional_wide = meridional_barrier(20, -36, -30; width=4) + @test meridional_wide.longitude == (18, 22) # Latitude band covering the full longitude range - band = BoundingBox(longitude=(-180.0, 180.0), latitude=(-60.0, -55.0)) - @test band.longitude == (-180.0, 180.0) - @test band.latitude == (-60.0, -55.0) + band = BoundingBox(longitude=(-180, 180), latitude=(-60, -55)) + @test band.longitude == (-180, 180) + @test band.latitude == (-60, -55) end @testset "Ocean basin labeling with barriers" begin From d1a51823976888f14c4c1be73de95cbc2ff095b0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 16:39:25 +0200 Subject: [PATCH 10/11] some fixes --- src/Bathymetry/label_ocean_basins.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/Bathymetry/label_ocean_basins.jl b/src/Bathymetry/label_ocean_basins.jl index eb3144fc5..58d0bb385 100644 --- a/src/Bathymetry/label_ocean_basins.jl +++ b/src/Bathymetry/label_ocean_basins.jl @@ -64,19 +64,21 @@ end # water cells near the boundary are correctly identified as connected. function copy_periodic_longitude(zb_cpu::Field, ::Periodic) Nλ = size(zb_cpu, 1) + half = Nλ ÷ 2 + + # 2D slice of the bathymetry at the bottom level. zb_data = zb_cpu.data[1:Nλ, :, 1] - zb_parent = zb_data.parent + zb_parent = parent(zb_data) # Concatenate a copy of the eastern half on the left and the western half on the right. # This is an O(Nλ × Nφ) CPU allocation, acceptable for the serial labeling step. - half = Nλ ÷ 2 - zb_parent = vcat(zb_parent[half:Nλ, :], zb_parent, zb_parent[1:half, :]) + concatenated_zb = vcat(zb_parent[half:Nλ, :], zb_parent, zb_parent[1:half, :]) - # Update offsets so that index 1 maps back to the original domain start - yoffsets = zb_cpu.data.offsets[2] - xoffsets = - half - 1 + # Offsets so that index 1 maps back to the original domain start + xoffset = - half - 1 + yoffset = zb_cpu.data.offsets[2] - return OffsetArray(zb_parent, xoffsets, yoffsets) + return OffsetArray(concatenated_zb, xoffset, yoffset) end copy_periodic_longitude(zb_cpu::Field, tx) = interior(zb_cpu, :, :, 1) @@ -191,10 +193,11 @@ function label_ocean_basins(grid::ImmersedBoundaryGrid; barriers=nothing) TX = topology(cpu_grid, 1) zb = cpu_grid.immersed_boundary.bottom_height - # If barriers are specified, apply them to a copy of the bathymetry + # If barriers are specified, apply them to a copy of the bathymetry. + # Using `similar(zb)` preserves the field type, which allows for future + # bathymetry representations beyond plain 2D `GridFittedBottom` fields. if !isnothing(barriers) - # Create a temporary field with the modified bathymetry - zb_modified = Field{Center, Center, Nothing}(cpu_grid) + zb_modified = similar(zb) parent(zb_modified) .= parent(zb) apply_barrier!(zb_modified, cpu_grid, barriers) else From dcf747ea28f212219358b2d103cde91b38c3b4fd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 19:42:56 +0200 Subject: [PATCH 11/11] fix bathymetry tests --- test/test_bathymetry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 57022def8..24ad9dd94 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -227,7 +227,7 @@ end # Test atlantic_ocean_basin creation atlantic = atlantic_ocean_basin(ibg) @test atlantic isa Basin - @test sum(atlantic.mask) > 0 # Should have some ocean cells + @test sum(interior(atlantic.mask)) > 0 # Should have some ocean cells mask = on_architecture(CPU(), atlantic.mask)