From 5dcbcf3860a7d53f7d374a0715c25996669006bb Mon Sep 17 00:00:00 2001 From: Finn Heukamp <57488986+FinnHeu@users.noreply.github.com> Date: Tue, 17 Mar 2026 17:49:31 +0100 Subject: [PATCH 1/3] Remove non-finite fill values for pf.ice_ext() calculation With the fill value change from 0 to nan in fesom2.7 the pf.ice_ext() did return incorrect values. --- pyfesom2/diagnostics.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index 14db4ee..502e781 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -76,7 +76,9 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): data = data.where(data < threshhold, 1) data = data.where(data > threshhold) - ext = (data[:, hemis_mask] * mesh.lump2[hemis_mask]).sum(axis=1) + finite_mask = data.squeeze().values + + ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) da = xr.DataArray( ext, dims=["time"], coords={"time": data.time}, name=varname, attrs=attrs ) @@ -84,10 +86,11 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): else: logger.debug(data) + finite_mask = data.squeeze() i, j = np.where(data < 0.15) data[:] = 1 data[i, j] = 0 - ext = (data[:, hemis_mask] * mesh.lump2[hemis_mask]).sum(axis=1) + ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) return ext From fa259053d9e1ace09e29b7234064082d6170e85e Mon Sep 17 00:00:00 2001 From: Finn Heukamp <57488986+FinnHeu@users.noreply.github.com> Date: Wed, 25 Mar 2026 11:22:44 +0100 Subject: [PATCH 2/3] Fix finite mask calculation and handle NaN values Replace finite mask calculation to handle NaN values and update data accordingly. --- pyfesom2/diagnostics.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index 502e781..e469e7d 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -76,7 +76,7 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): data = data.where(data < threshhold, 1) data = data.where(data > threshhold) - finite_mask = data.squeeze().values + finite_mask = np.isfinite(data.squeeze().values) ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) da = xr.DataArray( @@ -86,8 +86,9 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): else: logger.debug(data) - finite_mask = data.squeeze() - i, j = np.where(data < 0.15) + finite_mask = np.isfinite(data.squeeze()) + data = np.where(~finite_mask, 0, data) # replace nan values by 0s + i, j = np.where(data < threshhold) data[:] = 1 data[i, j] = 0 ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) From 60d9e9e090060af843dc63cf94e4b38998a0d52e Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 25 Mar 2026 14:46:10 +0100 Subject: [PATCH 3/3] Fix ice_ext to properly preserve NaN land cells MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous fix still had the critical bug where .where() operations converted NaN land cells to ice (value=1), causing 10x overestimation. This commit replaces the .where() logic with xr.where() that creates a proper binary field while explicitly preserving NaN values for land. Changes: - xarray path: Use xr.where() to create binary field, then preserve NaN - numpy path: Create binary field with np.where(), then restore NaN values - Both paths now correctly exclude land cells from ice extent calculation Tested with AWI-ESM3 data - values now show realistic ~20 million km² instead of inflated ~200 million km². --- pyfesom2/diagnostics.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index e469e7d..c24d329 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -73,12 +73,13 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): hemis_mask = mesh.y2 < 0 if isinstance(data, xr.DataArray): - data = data.where(data < threshhold, 1) - data = data.where(data > threshhold) + # Create binary ice extent field: 1 where ice >= threshold, 0 where < threshold, NaN stays NaN + ice_binary = xr.where(data >= threshhold, 1, 0) + ice_binary = ice_binary.where(~np.isnan(data)) # Preserve NaN values (land) - finite_mask = np.isfinite(data.squeeze().values) + finite_mask = ~np.isnan(ice_binary.squeeze().values) - ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) + ext = (ice_binary[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) da = xr.DataArray( ext, dims=["time"], coords={"time": data.time}, name=varname, attrs=attrs ) @@ -86,12 +87,12 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): else: logger.debug(data) - finite_mask = np.isfinite(data.squeeze()) - data = np.where(~finite_mask, 0, data) # replace nan values by 0s - i, j = np.where(data < threshhold) - data[:] = 1 - data[i, j] = 0 - ext = (data[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) + # Create binary ice extent field: 1 where ice >= threshold, 0 where < threshold, NaN stays NaN + ice_binary = np.where(data >= threshhold, 1, 0).astype(float) + ice_binary[np.isnan(data)] = np.nan # Preserve NaN values (land) + + finite_mask = ~np.isnan(ice_binary.squeeze()) + ext = (ice_binary[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) return ext