Skip to content

Allow the use of CopernicusMetadata to initialize sea ice #644

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open

Allow the use of CopernicusMetadata to initialize sea ice #644

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
bad6fb1
first fix, doesn't solve the nan issue yet
xkykai Oct 3, 2025
fb34727
write ValueInpainting, but leaves vestigial maxiter value
xkykai Oct 3, 2025
dee9a17
remove fill_nans functionaloty
xkykai Oct 3, 2025
a8df621
update default_inpainting to return ValueInpainting for sea ice metrics
xkykai Oct 4, 2025
59d3004
refactor location function for clarity and consistency in CopernicusM…
xkykai Oct 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 26 additions & 12 deletions src/DataWrangling/Copernicus/Copernicus.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@ import ClimaOcean.DataWrangling:
metadata_filename,
inpainted_metadata_path,
reversed_vertical_axis,
available_variables
available_variables,
is_three_dimensional

using Scratch

Expand Down Expand Up @@ -129,20 +130,33 @@ end

inpainted_metadata_path(metadata::CopernicusMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))

location(::CopernicusMetadata) = (Center, Center, Center)
copernicus_location(metadata::CopernicusMetadata) = is_three_dimensional(metadata) ? (Center, Center, Center) : (Center, Center, Nothing)

location(metadata::CopernicusMetadata) = copernicus_location(metadata)
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
Comment on lines +137 to +141
Copy link
Member

@glwagner glwagner Oct 3, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you should blacklist the 2D variables rather than whitelisting 3D variables, since most variables are 3D

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here I was trying to follow the design pattern here:

ClimaOcean.jl/src/DataWrangling/metadata_field.jl

Line 53 in dac9760

if is_three_dimensional(metadata)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can achieve this by using a ! conditional because is_three_dimensional == !is_two_dimensional

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for example:

is_three_dimensional(metadata::CopernicusMetadata) =
    !(metadata.name == :ice_concentration ||
      metadata.name == :ice_thickness)


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)
Comment on lines +145 to +146
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be better to build a Flat grid for 2D variables? Why do the z_interfaces matter?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this line seems to want to take something for z, I tried returning nothing but I guess it didn't work

ClimaOcean.jl/src/DataWrangling/metadata_field.jl

Line 35 in dac9760

grid = LatitudeLongitudeGrid(arch, FT; halo, longitude, latitude, z,

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah, the size and the halo must be 2-tuples to use nothing in z.

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
19 changes: 13 additions & 6 deletions src/DataWrangling/metadata_field.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@ function Field(metadata::Metadatum, arch=CPU();
inpainting = default_inpainting(metadata),
mask = nothing,
halo = (3, 3, 3),
cache_inpainted_data = true)
cache_inpainted_data = true,
fill_nans = nothing)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not just use true or false here? Also I am confused --- shouldn't this be handled by inpainting? Inpainting is generally the process of removing invalid entries. Usually we inpaint NaNs with fancy propagation method, but just replacing them with a neutral value (like 0) could be implemented as an additional inpainting option.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose you're right, we should do it via inpainting. Currently there is no zero inpainting available, so I suppose the solution is to write a zero inpainting method and then change default_inpainting so that it deals with CopernicusMetadata differently? It seems like for ECCO no inpainting is required for sea ice

ClimaOcean.jl/src/DataWrangling/DataWrangling.jl

Line 205 in dac9760

function default_inpainting(metadata)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the purpose of the default_inpainting is to allow dispatch on metadata type


download_dataset(metadata)

Expand Down Expand Up @@ -118,7 +119,7 @@ function Field(metadata::Metadatum, arch=CPU();
# Retrieve data from file according to metadata type
data = retrieve_data(metadata)

set_metadata_field!(field, data, metadata)
set_metadata_field!(field, data, metadata; fill_nans)
fill_halo_regions!(field)

if !isnothing(inpainting)
Expand Down Expand Up @@ -172,7 +173,7 @@ struct AverageNorthSouth end
@inline mangle(i, j, k, data, ::ShiftSouth) = @inbounds data[i, j-1, k]
@inline mangle(i, j, k, data, ::AverageNorthSouth) = @inbounds (data[i, j+1, k] + data[i, j, k]) / 2

function set_metadata_field!(field, data, metadatum)
function set_metadata_field!(field, data, metadatum; fill_nans)
grid = field.grid
arch = architecture(grid)

Expand Down Expand Up @@ -206,12 +207,13 @@ function set_metadata_field!(field, data, metadatum)
end

data = on_architecture(arch, data)
Oceananigans.Utils.launch!(arch, grid, spec, _kernel, field, data, mangling, temp_units, conc_units)
Oceananigans.Utils.launch!(arch, grid, spec, _kernel, field, data, mangling, temp_units, conc_units, fill_nans)

return nothing
end

@kernel function _set_2d_metadata_field!(field, data, mangling, temp_units, conc_units)
# 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, fill_nans)
i, j = @index(Global, NTuple)
d = mangle(i, j, data, mangling)

Expand All @@ -223,13 +225,16 @@ end
elseif !isnothing(conc_units)
d = convert_concentration(d, conc_units)
end
d = ifelse(isnothing(fill_nans), d, ifelse(isnan(d), convert(FT, 0), d))

d = convert_temperature(d, temp_units)
@inbounds field[i, j, 1] = d
end

@inline nan_convert_missing(FT, ::Missing) = convert(FT, NaN)
@inline nan_convert_missing(FT, d::Number) = convert(FT, d)

@kernel function _set_3d_metadata_field!(field, data, mangling, temp_units, conc_units)
@kernel function _set_3d_metadata_field!(field, data, mangling, temp_units, conc_units, fill_nans)
i, j, k = @index(Global, NTuple)
d = mangle(i, j, k, data, mangling)

Expand All @@ -241,6 +246,8 @@ end
elseif !isnothing(conc_units)
d = convert_concentration(d, conc_units)
end
d = ifelse(isnothing(fill_nans), d, ifelse(isnan(d), convert(FT, 0), d))

@inbounds field[i, j, k] = d
end

Expand Down
Loading