Skip to content
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.

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

Fix derivation for sea ice extent (siextent) #2648

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 1 addition & 1 deletion esmvalcore/cmor/tables/custom/CMOR_siextent.dat
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ modeling_realm: seaIce
! Variable attributes:
!----------------------------------
standard_name:
units: m2
units: 1
cell_methods: area: mean where sea time: mean
cell_measures: area: areacello
long_name: Sea Ice Extent
Expand Down
44 changes: 28 additions & 16 deletions esmvalcore/preprocessor/_derive/siextent.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
"""Derivation of variable `sithick`."""
"""Create mask for derivation of variable `siextent`."""

import logging

Expand All @@ -14,27 +14,33 @@


class DerivedVariable(DerivedVariableBase):
"""Derivation of variable `siextent`."""
"""Create mask for derivation of variable `siextent`."""

@staticmethod
def required(project):
"""Declare the variables needed for derivation."""
required = [
{"short_name": "sic", "optional": "true"},
{"short_name": "siconca", "optional": "true"},
]
"""Declare the variable needed for derivation."""
# 'sic' is sufficient as there is already an entry
# in the mapping table esmvalcore/cmor/variable_alt_names.yml
if project == "CMIP5":
required = [{"short_name": "sic"}]
elif project == "CMIP6":
required = [
{"short_name": "sic", "optional": True},
{"short_name": "siconca", "optional": True},
]

return required

@staticmethod
def calculate(cubes):
"""Compute sea ice extent.

Returns an array of ones in every grid point where
the sea ice area fraction has values > 15 .
the sea ice area fraction has values > 15% .

Use in combination with the preprocessor
`area_statistics(operator='sum')` to weigh by the area and
compute global or regional sea ice extent values.
compute global or regional sea ice extent values (in m2).

Arguments
---------
Expand All @@ -48,16 +54,22 @@ def calculate(cubes):
sic = cubes.extract_cube(Constraint(name="sic"))
except iris.exceptions.ConstraintMismatchError:
try:
sic = cubes.extract_cube(Constraint(name="siconca"))
Copy link
Contributor

Choose a reason for hiding this comment

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

If I remember correctly, siconca was added because some models were missing siconc, but I am not sure if that is still the case.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is right. In any case, I think it does not hurt to also try "siconca", so I would prefer to keep this here if that's fine.

Copy link
Contributor

Choose a reason for hiding this comment

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

Well in that case, it looks like it needs to be re-added because tests are failing due to siconca not being required anymore.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I justed updated that part a bit: 0c8beb9
I cannot get this preprocessor working for CMIP5 data, though, when adding {"short_name": "siconca", "optional": "true"},. Alternatively, I can remove the "siconca" part. Any advice would be highly appreciated...

Copy link
Contributor

@sloosvel sloosvel Feb 21, 2025

Choose a reason for hiding this comment

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

I tested with the latest commit 1eaa6f0 loading CMIP5 sic data, CMIP6 sic data and CMIP6 data that only has siconca available and it worked finding the data that is needed for each project:

datasets:

        - {dataset: GISS-E2-1-H, grid: gr} # CMIP6 sic data 
        - {dataset: GISS-E2-1-H, exp: piControl, grid: gn, timerange: '3180/3180'} # CMIP6 siconca data
        - {dataset: GISS-E2-H-CC, project: CMIP5, ensemble: r1i1p1, mip: OImon} # CMIP5 sic data

diagnostics:

  test:
    variables:
      siextent:
        project: CMIP6
        mip: SImon
        timerange: '2000/2000'
        derive: true
        exp: historical
        ensemble: r1i1p1f1
    scripts: null

Let me know if it works for you

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thank you for looking into this, @sloosvel! I tried this and I don't think we found the optimum solution yet. If I add a CMIP6 dataset that provides both, siconc and siconca (e.g. MPI-ESM1-2-LR), I run into a shape error, e.g.

ValueError: Chunks do not add up to shape. Got chunks=((96,), (192,)), shape=(220, 256)

Also, our current solution does not support to process any observationally-based data (e.g. projects OBS, OBS6, ana4mips, obs4MIPs, native5, etc.). Here are some examples for observationally-based datasets that I tried:

  - {dataset: ESACCI-SEAICE, project: OBS6, tier: 2, type: sat, version: L4-SICONC-RE-SSMI-12.5kmEASE2-fv3.0-NH,
     supplementary_variables: [{short_name: areacello, mip: Ofx}]}
  - {dataset: HadISST, project: OBS, tier: 2, type: reanaly, version: '1', mip: OImon}
  - {dataset: CFSR, project: ana4mips, tier: 1, type: reanalysis, mip: OImon}

So maybe checking for if project == 'CMIP6' or project == 'OBS6' is enough and a plain else for all other cases in the required function? But then, there is still the shape problem.

Do you have an idea what we could do? I didn't expect this to be so complicated...

except iris.exceptions.ConstraintMismatchError as exc:
raise RecipeError(
"Derivation of siextent failed due to missing variables "
"sic and siconca."
) from exc
sic = cubes.extract_cube(Constraint(name="siconc"))
except iris.exceptions.ConstraintMismatchError:
try:
sic = cubes.extract_cube(Constraint(name="siconca"))
except iris.exceptions.ConstraintMismatchError as exc:
raise RecipeError(
"Derivation of siextent failed due to missing variables "
"sic and siconc/siconca."
) from exc

ones = da.ones_like(sic)
siextent_data = da.ma.masked_where(sic.lazy_data() < 15.0, ones)
siextent = sic.copy(siextent_data)
siextent.units = "m2"
# unit is 1 as this is just a mask that has to be used with
# preprocessor area_statistics(operator='sum') to obtain the
# sea ice extent (m2)
siextent.units = "1"

return siextent
6 changes: 3 additions & 3 deletions tests/unit/preprocessor/_derive/test_siextent.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ def test_siextent_calculation_sic(cubes_sic):
"""Test function ``calculate`` when sic is available."""
derived_var = siextent.DerivedVariable()
out_cube = derived_var.calculate(cubes_sic)
assert out_cube.units == cf_units.Unit("m2")
assert out_cube.units == cf_units.Unit("1")
out_data = out_cube.data
expected = np.ma.ones_like(cubes_sic[0].data)
expected.mask = True
Expand All @@ -75,7 +75,7 @@ def test_siextent_calculation_siconca(cubes_siconca):
"""Test function ``calculate`` when siconca is available."""
derived_var = siextent.DerivedVariable()
out_cube = derived_var.calculate(cubes_siconca)
assert out_cube.units == cf_units.Unit("m2")
assert out_cube.units == cf_units.Unit("1")
out_data = out_cube.data
expected = np.ma.ones_like(cubes_siconca[0].data)
expected.mask = True
Expand All @@ -88,7 +88,7 @@ def test_siextent_calculation(cubes):
"""Test function ``calculate`` when sic and siconca are available."""
derived_var = siextent.DerivedVariable()
out_cube = derived_var.calculate(cubes)
assert out_cube.units == cf_units.Unit("m2")
assert out_cube.units == cf_units.Unit("1")
out_data = out_cube.data
expected = np.ma.ones_like(cubes[0].data)
expected.mask = True
Expand Down