diff --git a/src/DataWrangling/Copernicus/Copernicus.jl b/src/DataWrangling/Copernicus/Copernicus.jl index 2ee7f0116..8c7ac0dce 100644 --- a/src/DataWrangling/Copernicus/Copernicus.jl +++ b/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 \ No newline at end of file diff --git a/src/DataWrangling/DataWrangling.jl b/src/DataWrangling/DataWrangling.jl index e5da9795f..ca0d7ced2 100644 --- a/src/DataWrangling/DataWrangling.jl +++ b/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 diff --git a/src/DataWrangling/inpainting.jl b/src/DataWrangling/inpainting.jl index dde0f1f91..93414de91 100644 --- a/src/DataWrangling/inpainting.jl +++ b/src/DataWrangling/inpainting.jl @@ -11,6 +11,30 @@ 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) @@ -18,6 +42,9 @@ function propagating(field, mask, iter, inpainting::NearestNeighborInpainting) 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)) diff --git a/src/DataWrangling/metadata_field.jl b/src/DataWrangling/metadata_field.jl index a57b6622c..1802585d3 100644 --- a/src/DataWrangling/metadata_field.jl +++ b/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