From bfb95b0e04d77365df7233abf18fdf6c2e7bafd7 Mon Sep 17 00:00:00 2001 From: Mukulika Pahari Date: Mon, 11 May 2026 15:09:30 +0200 Subject: [PATCH 1/3] update ice_area to handle NaNs --- pyfesom2/diagnostics.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index c24d329..2010073 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -140,8 +140,8 @@ def ice_vol(data, mesh, hemisphere="N", attrs={}): return vol -def ice_area(data, mesh, hemisphere="N", attrs={}): - """ Compute sea ice volume. +def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): + """ Compute sea ice area. Parameters ---------- @@ -175,13 +175,26 @@ def ice_area(data, mesh, hemisphere="N", attrs={}): hemis_mask = mesh.y2 < 0 if isinstance(data, xr.DataArray): - vol = (data[:, hemis_mask] * mesh.lump2[hemis_mask]).sum(axis=1) + # Create ice area field: data where ice >= threshold, 0 where < threshold, NaN stays NaN + ice_field = xr.where(data >= threshhold, data, 0) + ice_field = ice_field.where(~np.isnan(data)) # Preserve NaN values (land) + + finite_mask = ~np.isnan(ice_field.squeeze().values) + + area = (ice_field[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) da = xr.DataArray( - vol, dims=["time"], coords={"time": data.time}, name=varname, attrs=attrs + area, dims=["time"], coords={"time": data.time}, name=varname, attrs=attrs ) return da + else: - area = (data[:, hemis_mask] * mesh.lump2[hemis_mask]).sum(axis=1) + # logger.debug(data) + # Create ice area field: data where ice >= threshold, 0 where < threshold, NaN stays NaN + ice_field = np.where(data >= threshhold, data, 0).astype(float) + ice_field[np.isnan(data)] = np.nan # Preserve NaN values (land) + + finite_mask = ~np.isnan(ice_field.squeeze()) + area = (ice_field[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) return area From f10f759f4364bd3d4ebf5735c9760748e9b089e6 Mon Sep 17 00:00:00 2001 From: Paul Gierz Date: Tue, 12 May 2026 12:19:39 +0200 Subject: [PATCH 2/3] Drop stray debug comment and trailing whitespace in ice_area --- pyfesom2/diagnostics.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index 2010073..251a7a9 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -180,7 +180,7 @@ def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): ice_field = ice_field.where(~np.isnan(data)) # Preserve NaN values (land) finite_mask = ~np.isnan(ice_field.squeeze().values) - + area = (ice_field[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) da = xr.DataArray( area, dims=["time"], coords={"time": data.time}, name=varname, attrs=attrs @@ -188,11 +188,10 @@ def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): return da else: - # logger.debug(data) # Create ice area field: data where ice >= threshold, 0 where < threshold, NaN stays NaN ice_field = np.where(data >= threshhold, data, 0).astype(float) ice_field[np.isnan(data)] = np.nan # Preserve NaN values (land) - + finite_mask = ~np.isnan(ice_field.squeeze()) area = (ice_field[:, hemis_mask & finite_mask] * mesh.lump2[hemis_mask & finite_mask]).sum(axis=1) return area From b908b0d2be7e773115bf363ff22c47ca559827ca Mon Sep 17 00:00:00 2001 From: Mukulika Pahari Date: Tue, 12 May 2026 13:23:17 +0200 Subject: [PATCH 3/3] change default threshold and fix spelling --- pyfesom2/diagnostics.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyfesom2/diagnostics.py b/pyfesom2/diagnostics.py index 251a7a9..17094b2 100644 --- a/pyfesom2/diagnostics.py +++ b/pyfesom2/diagnostics.py @@ -37,7 +37,7 @@ def add_timedim(data, date="1970-01-01"): return data -def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): +def ice_ext(data, mesh, hemisphere="N", threshold=0.15, attrs={}): """ Compute sea ice extent. Parameters @@ -50,7 +50,7 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): FESOM2 mesh object. hemisphere: str can be ether 'N' or 'S' - threshhold: float + threshold: float default is 0.15 attrs: dict dictionary of attributes that will be put in the resulting @@ -74,7 +74,7 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): if isinstance(data, xr.DataArray): # 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 = xr.where(data >= threshold, 1, 0) ice_binary = ice_binary.where(~np.isnan(data)) # Preserve NaN values (land) finite_mask = ~np.isnan(ice_binary.squeeze().values) @@ -88,7 +88,7 @@ def ice_ext(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): else: logger.debug(data) # 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.where(data >= threshold, 1, 0).astype(float) ice_binary[np.isnan(data)] = np.nan # Preserve NaN values (land) finite_mask = ~np.isnan(ice_binary.squeeze()) @@ -140,7 +140,7 @@ def ice_vol(data, mesh, hemisphere="N", attrs={}): return vol -def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): +def ice_area(data, mesh, hemisphere="N", threshold=0., attrs={}): """ Compute sea ice area. Parameters @@ -176,7 +176,7 @@ def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): if isinstance(data, xr.DataArray): # Create ice area field: data where ice >= threshold, 0 where < threshold, NaN stays NaN - ice_field = xr.where(data >= threshhold, data, 0) + ice_field = xr.where(data >= threshold, data, 0) ice_field = ice_field.where(~np.isnan(data)) # Preserve NaN values (land) finite_mask = ~np.isnan(ice_field.squeeze().values) @@ -189,7 +189,7 @@ def ice_area(data, mesh, hemisphere="N", threshhold=0.15, attrs={}): else: # Create ice area field: data where ice >= threshold, 0 where < threshold, NaN stays NaN - ice_field = np.where(data >= threshhold, data, 0).astype(float) + ice_field = np.where(data >= threshold, data, 0).astype(float) ice_field[np.isnan(data)] = np.nan # Preserve NaN values (land) finite_mask = ~np.isnan(ice_field.squeeze())