Allow the use of CopernicusMetadata to initialize sea ice

src/DataWrangling/Copernicus/Copernicus.jl

@@ -22,7 +22,8 @@ import ClimaOcean.DataWrangling:
     metadata_filename,
     inpainted_metadata_path,
     reversed_vertical_axis,
-    available_variables
+    available_variables,
+    is_three_dimensional
 
 using Scratch
 
@@ -129,20 +130,31 @@ end
 
 inpainted_metadata_path(metadata::CopernicusMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
 
-location(::CopernicusMetadata) = (Center, Center, Center)
+location(metadata::CopernicusMetadata) = is_three_dimensional(metadata) ? (Center, Center, Center) : (Center, Center, Nothing)
 longitude_interfaces(::CopernicusMetadata) = (0, 360)
 latitude_interfaces(::CopernicusMetadata) = (-80, 90)
 
+is_three_dimensional(metadata::CopernicusMetadata) =
+    metadata.name == :temperature ||
+    metadata.name == :salinity ||
+    metadata.name == :u_velocity ||
+    metadata.name == :v_velocity
+
 function z_interfaces(metadata::CopernicusMetadata)
-    path = metadata_path(metadata)
-    ds = Dataset(path)
-    zc = - reverse(ds["depth"][:])
-    close(ds)
-    dz = zc[2] - zc[1]
-    zf = zc[1:end-1] .+ zc[2:end]
-    push!(zf, 0)
-    pushfirst!(zf, zf[1] - dz)
-    return zf
+    # Check if this is a 2D surface variable
+    if !is_three_dimensional(metadata)
+        return (0, 1)
+    else
+        path = metadata_path(metadata)
+        ds = Dataset(path)
+        zc = - reverse(ds["depth"][:])
+        close(ds)
+        dz = zc[2] - zc[1]
+        zf = zc[1:end-1] .+ zc[2:end]
+        push!(zf, 0)
+        pushfirst!(zf, zf[1] - dz)
+        return zf
+    end
 end
 
-end # module Copernicus
+end # module Copernicus

src/DataWrangling/DataWrangling.jl

@@ -206,7 +206,7 @@ function default_inpainting(metadata)
     if metadata.name in (:temperature, :salinity)
         return NearestNeighborInpainting(Inf)
     elseif metadata.name in (:sea_ice_thickness, :sea_ice_concentration)
-        return nothing
+        return ValueInpainting(0)
     else
         return NearestNeighborInpainting(5)
     end

src/DataWrangling/inpainting.jl

@@ -11,13 +11,40 @@ struct NearestNeighborInpainting{M}
     maxiter :: M
 end
 
+"""
+    ValueInpainting{V}
+
+A structure representing a simple value inpainting algorithm, where all missing values are
+substituted with a specified constant value.
+
+# Fields
+- `value :: V`: The constant value to use for inpainting missing data.
+- `maxiter`: For cache compatibility (not used in the algorithm itself). Defaults to 0.
+
+# Example
+```julia
+inpainting = ValueInpainting(0.0)  # Fill missing values with 0
+inpaint_mask!(field, mask; inpainting)
+```
+"""
+struct ValueInpainting{V, M}
+    value :: V
+    maxiter :: M
+end
+
+# Convenience constructor that sets maxiter to a default value
+ValueInpainting(value) = ValueInpainting(value, 0)
+
 propagate_horizontally!(field, ::Nothing, substituting_field=deepcopy(field); kw...) = field
 
 function propagating(field, mask, iter, inpainting::NearestNeighborInpainting)
     nans = sum(isnan, field; condition=interior(mask))
     return nans > 0 && iter < inpainting.maxiter
 end
 
+# ValueInpainting doesn't need iteration
+propagating(field, mask, iter, inpainting::ValueInpainting) = false
+
 """
     propagate_horizontally!(inpainting, field, mask [, substituting_field=deepcopy(field)])
 
@@ -57,6 +84,17 @@ function propagate_horizontally!(inpainting::NearestNeighborInpainting, field, m
     return field
 end
 
+function propagate_horizontally!(inpainting::ValueInpainting, field, mask,
+                                 substituting_field=deepcopy(field))
+    grid = field.grid
+    arch = architecture(grid)
+
+    launch!(arch, grid, size(field), _fill_with_number!, field, mask, inpainting.value)
+    fill_halo_regions!(field)
+
+    return field
+end
+
 # Maybe we can remove this propagate field in lieu of a diffusion,
 # Still we'll need to do this a couple of steps on the original grid
 @kernel function _propagate_field!(substituting_field, ::NearestNeighborInpainting, field)
@@ -102,6 +140,12 @@ end
     @inbounds field[i, j, k] *= !isnan(field[i, j, k])
 end
 
+@kernel function _fill_with_number!(field, mask, value)
+    i, j, k = @index(Global, NTuple)
+    FT_value = convert(eltype(field), value)
+    @inbounds field[i, j, k] = ifelse(mask[i, j, k], FT_value, field[i, j, k])
+end
+
 """
     inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf))
 
@@ -117,10 +161,24 @@ Arguments
 - `mask`: Boolean-valued `Field`, values where
           `mask[i, j, k] == true` are inpainted.
 
-- `inpainting`: The inpainting algorithm to use. The only option is
-                `NearestNeighborInpainting(maxiter)`, where an average
-                of the valid surrounding values is used `maxiter` times.
-                Default: `NearestNeighborInpainting(Inf)`.
+- `inpainting`: The inpainting algorithm to use. Options include:
+                * `NearestNeighborInpainting(maxiter)`: Uses an average
+                  of the valid surrounding values, repeated `maxiter` times.
+                  Default: `NearestNeighborInpainting(Inf)`.
+                * `ValueInpainting(value)`: Fills all missing values with
+                  the specified constant `value`.
+
+# Examples
+```julia
+# Use nearest neighbor inpainting with default settings
+inpaint_mask!(field, mask)
+
+# Fill missing values with zero
+inpaint_mask!(field, mask; inpainting=ValueInpainting(0))
+
+# Fill missing values with a specific temperature
+inpaint_mask!(field, mask; inpainting=ValueInpainting(273.15))
+```
 """
 function inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf))
 

src/DataWrangling/metadata_field.jl

@@ -211,6 +211,7 @@ function set_metadata_field!(field, data, metadatum)
     return nothing
 end
 
+# TODO: sort out elegant way to handle NaNs for sea ice data in GLORYS
 @kernel function _set_2d_metadata_field!(field, data, mangling, temp_units, conc_units)
     i, j = @index(Global, NTuple)
     d = mangle(i, j, data, mangling)
@@ -223,6 +224,8 @@ end
     elseif !isnothing(conc_units)
         d = convert_concentration(d, conc_units)
     end
+
+    d = convert_temperature(d, temp_units)
     @inbounds field[i, j, 1] = d
 end
 
@@ -241,6 +244,7 @@ end
     elseif !isnothing(conc_units)
         d = convert_concentration(d, conc_units)
     end
+
     @inbounds field[i, j, k] = d
 end