Addition of bathymetry modification utilities

src/Bathymetry/regrid_bathymetry.jl

@@ -517,3 +517,92 @@ function remove_minor_basins!(zb, keep_major_basins, core_size)
 
     return nothing
 end
+
+@kernel function _modify_bathymetry_kernel!(bathymetry, lons, lats, 
+                                            target_lon, target_lat, target_depth, 
+                                            rx_km, ry_km, deg_buffer, R)
+    i, j = @index(Global, NTuple)
+
+    @inbounds begin
+        # 1. Bounding box check (fast rejection to avoid expensive trig)
+        lat_diff = abs(lats[i, j] - target_lat)
+        lon_diff = abs(mod(lons[i, j] - target_lon + 180, 360) - 180)
+
+        if lat_diff <= deg_buffer && lon_diff <= deg_buffer
+
+            # 2. True physical distances in km
+            dy = R * deg2rad(lats[i, j] - target_lat)
+
+            # dx: East-West distance (accounting for convergence of meridians)
+            avg_lat = (lats[i, j] + target_lat) / 2.0
+            dλ = mod(lons[i, j] - target_lon + 180, 360) - 180
+            dx = R * cosd(avg_lat) * deg2rad(dλ)
+
+            # 3. Ellipsoidal distance (normalized)
+            dist_sq = (dx / rx_km)^2 + (dy / ry_km)^2
+
+            if dist_sq < 1.0
+                dist = sqrt(dist_sq)
+
+                # Smooth taper (Cosine) blending
+                weight = 0.5 * (1.0 + cos(π * dist))
+
+                # Blend the original bathymetry with the target depth
+                original_depth = bathymetry[i, j, 1]
+                bathymetry[i, j, 1] = (weight * target_depth) + ((1.0 - weight) * original_depth)
+            end
+        end
+    end
+end
+
+"""
+    modify_bathymetry_depth!(grid, bathymetry, lons, lats, target_lon, target_lat, target_depth; rx_km=150.0, ry_km=150.0, R=6371.0)
+
+Sets all points of `bathymetry` within `rx_km` and `ry_km` (in kilometers) of the target longitude `target_lon` and target latitude `target_lat` to `target_depth` (in meters) over the coordinates `lons` and `lats`. The radius of the planet can be set with the keyword argument `R` (default is Earth's radius of 6371 km).
+
+Arguments
+=========
+
+- `grid`: The Oceananigans grid object.
+- `bathymetry`: 2D array (or interior view) of bathymetric depths (negative values).
+- `lons`: 2D array of longitudes corresponding to `bathymetry`.
+- `lats`: 2D array of latitudes corresponding to `bathymetry`.
+- `target_lon`: Longitude of the strait center.
+- `target_lat`: Latitude of the strait center.
+- `target_depth`: Depth to set within the strait (negative value, e.g., -500 for 500 m depth).
+
+Keyword Arguments
+=================
+- `rx_km`: Radius in kilometers in the east-west direction. Default: 150.0 km.
+- `ry_km`: Radius in kilometers in the north-south direction. Default: 150.0 km. 
+- `R`: Earth's radius in kilometers. Default: 6371.0 km.
+"""
+function modify_bathymetry_depth!(grid, bathymetry, lons, lats, target_lon, target_lat, target_depth; 
+                                  rx_km=150.0, ry_km=150.0, R=6371.0)
+
+    arch = architecture(grid)
+
+    # Ensure inputs are appropriate types for precision mapping
+    FT = eltype(grid)
+    t_lon, t_lat, t_depth = FT(target_lon), FT(target_lat), FT(target_depth)
+    rx, ry, R_ft = FT(rx_km), FT(ry_km), FT(R)
+
+    # CPU Buffer calculation (broad bounding box in degrees)
+    deg_buffer = (max(rx, ry) / 111.0) * 2.0
+
+    # SAFETY CHECK for halo trap:
+    # grid.Nx is interior size, grid.Hx is halo size
+    full_size_x = grid.Nx + 2 * grid.Hx
+
+    if size(bathymetry, 1) == full_size_x && size(lons, 1) == grid.Nx
+        @warn "Mismatch detected: `bathymetry` includes halos but `lons` does not. " *
+              "The carving will be shifted! Pass `interior(bathymetry)` instead."
+    end
+
+    # Launch the kernel
+    launch!(arch, grid, :xy, _modify_bathymetry_kernel!, 
+            bathymetry, lons, lats, t_lon, t_lat, t_depth, 
+            rx, ry, deg_buffer, R_ft)
+
+    return nothing
+end

src/NumericalEarth.jl

@@ -28,6 +28,7 @@ export
     JRA55PrescribedAtmosphere,
     JRA55NetCDFBackend,
     regrid_bathymetry,
+    modify_bathymetry_depth!,
     Metadata,
     Metadatum,
     ECCOMetadatum,

test/test_bathymetry.jl

@@ -5,7 +5,8 @@ using NumericalEarth.Bathymetry: remove_minor_basins!,
                                  BathymetryRegridding,
                                  cache_filename,
                                  load_bathymetry_cache,
-                                 save_bathymetry_cache
+                                 save_bathymetry_cache,
+                                 modify_bathymetry_depth!
 using NumericalEarth.DataWrangling.ETOPO
 using Statistics
 
@@ -140,3 +141,86 @@ end
     result4 = regrid_bathymetry(grid; cache=false)
     @test parent(result1) == parent(result4)
 end
+
+
+@testset "Manual Bathymetry Carving (modify_bathymetry_depth!)" begin
+    @info "Testing manual bathymetry carving on a global 1-degree grid..."
+
+    # Fallback just in case test_architectures isn't globally defined in the test suite
+    archs = @isdefined(test_architectures) ? test_architectures : [CPU()]
+
+    for arch in archs
+        # 1. Setup a global 1-degree grid to precisely match the simulation environment
+        #    Range: -180 to 180 degrees (360 cells), -90 to 90 degrees (180 cells)
+        grid = LatitudeLongitudeGrid(arch;
+                                     size = (360, 180, 1),
+                                     longitude = (-180, 180),
+                                     latitude = (-90, 90),
+                                     z = (0, 1))
+
+        # Initialize flat bathymetry at -1000m
+        h = Field{Center, Center, Nothing}(grid)
+        set!(h, -1000.0)
+
+        # 2. Prepare 2D Coordinate Arrays using robust mathematical bounds
+        #    This completely bypasses any xnode/ynode API quirks across Oceananigans versions
+        #    and ensures the coordinate matrices perfectly map to the 1-degree grid.
+        lons_2d = zeros(grid.Nx, grid.Ny)
+        lats_2d = zeros(grid.Nx, grid.Ny)
+
+        for j in 1:grid.Ny, i in 1:grid.Nx
+            lons_2d[i, j] = -180.0 + (i - 0.5) * (360.0 / grid.Nx)
+            lats_2d[i, j] =  -90.0 + (j - 0.5) * (180.0 / grid.Ny)
+        end
+
+        # Move coordinates to the testing architecture (GPU/CPU)
+        lons_dev = on_architecture(arch, lons_2d)
+        lats_dev = on_architecture(arch, lats_2d)
+
+        # Target the exact center of the global map (near Equator / Prime Meridian)
+        target_i, target_j = 180, 90
+        target_lon = lons_2d[target_i, target_j]
+        target_lat = lats_2d[target_i, target_j]
+
+        # --- TEST 1: Basic Carving & Depth Accuracy ---
+        modify_bathymetry_depth!(grid, h, lons_dev, lats_dev, target_lon, target_lat, -100.0; 
+                                 rx_km=150.0, ry_km=150.0)
+
+        # Pulling Array safely fetches the data back to CPU
+        h_cpu = Array(interior(h, :, :, 1))
+
+        # Because we targeted the exact cell center, distance = 0, so weight = 1.0
+        center_val = h_cpu[target_i, target_j] 
+        @test center_val ≈ -100.0 atol=1.0
+
+        # Points far away should remain un-carved
+        corner_val = h_cpu[1, 1]
+        @test corner_val ≈ -1000.0
+
+        # --- TEST 2: Ellipsoidal Anisotropy (Stretching) ---
+        set!(h, -1000.0) # Reset Bathymetry
+
+        # Carve an ellipsoid: Wide in X (400km), Narrow in Y (50km)
+        modify_bathymetry_depth!(grid, h, lons_dev, lats_dev, target_lon, target_lat, -100.0; 
+                                 rx_km=400.0, ry_km=50.0)
+
+        h_cpu = Array(interior(h, :, :, 1))
+
+        # Point A: East (+2 cells at 1.0 deg/cell = 2.0 deg = ~222km). 
+        # Inside the 400km rx_km radius -> Should be modified
+        val_east = h_cpu[target_i + 2, target_j] 
+        @test val_east > -1000.0 
+
+        # Point B: North (+1 cell at 1.0 deg/cell = 1.0 deg = ~111km). 
+        # Outside the 50km ry_km radius -> Should remain unmodified
+        val_north = h_cpu[target_i, target_j + 1]
+        @test val_north ≈ -1000.0
+
+        # --- TEST 3: Smooth Tapering ---
+        # Check a point halfway inside the East-West radius (+1 cell = 1.0 deg = ~111km)
+        val_mid = h_cpu[target_i + 1, target_j] 
+
+        # It should be mathematically blended smoothly between -100 and -1000
+        @test -1000.0 < val_mid < -100.0
+    end
+end