Review Mapping, add support for seaice variables

src/access_moppy/derivations/__init__.py

@@ -8,6 +8,14 @@
     calc_topsoil,
     extract_tilefrac,
 )
+from access_moppy.derivations.calc_ocean import (
+    calc_areacello,
+    calc_global_ave_ocean,
+)
+from access_moppy.derivations.calc_seaice import (
+    calc_seaice_extent,
+    calc_hemi_seaice,
+)
 
 custom_functions = {
     "add": lambda *args: reduce(operator.add, args),
@@ -24,6 +32,10 @@
     "calc_landcover": calc_landcover,
     "extract_tilefrac": extract_tilefrac,
     "average_tile": average_tile,
+    "calc_areacello": calc_areacello,
+    "calc_global_ave_ocean": calc_global_ave_ocean,
+    "calc_seaice_extent": calc_seaice_extent,
+    "calc_hemi_seaice": calc_hemi_seaice,
 }
 
 

src/access_moppy/derivations/calc_ocean.py

@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+# This file contains a collection of functions to calculate ocean
+# derived variables from ACCESS model output.
+#
+# To propose new calculations and/or update to existing ones see documentation:
+#
+# and open a new issue on github.
+
+import numpy as np
+import xarray as xr
+
+
+def calc_areacello(area, mask_v):
+    """
+    Calculate ocean grid cell area with proper masking.
+
+    This function applies ocean mask to the grid cell area array, ensuring that
+    only valid ocean cells are included in area calculations. Land and invalid
+    grid cells are set to zero.
+
+    For ACCESS-ESM1.6, this function is typically used with the ocean-2d-ht.nc
+    mask file, where grid cells with zero depth indicate land areas that should
+    be excluded from ocean calculations.
+
+    Parameters
+    ----------
+    area : array-like
+        Grid cell area values (m²).
+        Can be numpy masked array or regular array.
+    mask_v : array-like
+        Ocean mask array where:
+        - Valid ocean cells: unmasked/True values (depth > 0)
+        - Land or invalid cells: masked/False values (depth = 0)
+        Must have compatible dimensions with area.
+        For ACCESS-ESM1.6: typically derived from ocean-2d-ht.nc depth field.
+
+    Returns
+    -------
+    array-like
+        Ocean grid cell area with land cells set to zero.
+        - Units: m²
+        - Shape: Same as input area
+        - Type: Filled array (no masked values)
+
+    Examples
+    --------
+    Apply ocean mask to calculate valid ocean cell areas:
+
+    >>> ocean_areas = calc_areacello(grid_area, ocean_mask)
+
+    ACCESS-ESM1.6 typical usage with depth-based masking:
+
+    >>> # Load ocean depth from ACCESS-ESM1.6 auxiliary file
+    >>> depth = load_data('ocean-2d-ht.nc')['ht']
+    >>> ocean_mask = depth > 0  # Create mask where depth > 0 (ocean)
+    >>> ocean_areas = calc_areacello(grid_area, ocean_mask)
+
+    Notes
+    -----
+    - Land and invalid grid cells are filled with zeros
+    - The function preserves the original area values for valid ocean cells
+    - Output is a regular (non-masked) array for consistent downstream processing
+    - Commonly used for ocean budget calculations and spatial integration
+
+    ACCESS-ESM1.6 Implementation Details:
+    - Uses ocean-2d-ht.nc file containing ocean depth (ht variable)
+    - Grid cells with depth = 0 are identified as land areas
+    - Grid cells with depth > 0 are valid ocean areas
+    - This masking ensures accurate ocean-only area calculations
+    """
+    area.mask = mask_v.mask
+    return area.filled(0)
+
+
+def calc_global_ave_ocean(var, rho_dzt, area_t):
+    """
+    Calculate global volume-weighted average of ocean variables.
+
+    This function computes the global ocean average of a variable by weighting
+    each grid cell by its mass (density × thickness × area). This provides
+    a proper volume-weighted mean that accounts for the varying cell sizes
+    and water column thickness in the ocean model.
+
+    Parameters
+    ----------
+    var : array-like
+        Ocean variable to average (e.g., temperature, salinity).
+        Can be 3D (time, depth, lat, lon) or 2D (time, lat, lon).
+        Should have compatible dimensions with rho_dzt and area_t.
+    rho_dzt : array-like
+        Ocean cell mass per unit area (kg/m²).
+        Typically calculated as: density × cell_thickness × area
+        Must have dimensions compatible with var for broadcasting.
+    area_t : array-like
+        Ocean grid cell area (m²).
+        Must have compatible horizontal dimensions with var.
+
+    Returns
+    -------
+    array-like
+        Global volume-weighted average of the input variable.
+        - Shape: (time,) - preserves time dimension only
+        - Units: Same as input variable
+        - Type: 1D array with time series of global averages
+
+    Examples
+    --------
+    Calculate global mean ocean temperature:
+
+    >>> global_temp = calc_global_ave_ocean(temperature, rho_dzt, area_t)
+
+    Calculate global mean sea surface temperature:
+
+    >>> global_sst = calc_global_ave_ocean(sst_2d, rho_dzt_surface, area_t)
+
+    Notes
+    -----
+    - Uses mass-weighted averaging for proper volume representation
+    - Automatically handles both 3D and 2D input variables
+    - For 3D variables: averages over depth, latitude, and longitude (axes 1,2,3)
+    - For 2D variables: averages over latitude and longitude (axes 1,2)
+    - Mass weighting accounts for varying cell volumes in ocean models
+
+    ACCESS-ESM1.6 Implementation Details:
+    - rho_dzt typically from ocean model output (density × thickness)
+    - area_t from ocean grid specification
+    - Commonly used for global ocean heat content, salt content calculations
+    - Handles irregular ocean grid geometries and partial cells
+    """
+    # Calculate mass weighting field
+    mass = rho_dzt * area_t
+
+    # Debug: print shape information
+    print("Variable shape:", np.shape(var))
+
+    try:
+        # Try 3D averaging (time, depth, lat, lon) -> (time,)
+        vnew = np.average(var, axis=(1, 2, 3), weights=mass)
+    except:
+        # Fall back to 2D averaging (time, lat, lon) -> (time,)
+        # Use only surface layer of mass weights
+        vnew = np.average(var, axis=(1, 2), weights=mass[:, 0, :, :])
+
+    return vnew

src/access_moppy/derivations/calc_seaice.py

@@ -0,0 +1,175 @@
+#!/usr/bin/env python
+# This file contains a collection of functions to calculate sea ice
+# derived variables from ACCESS model output.
+#
+# To propose new calculations and/or update to existing ones see documentation:
+#
+# and open a new issue on github.
+
+import logging
+import numpy as np
+import xarray as xr
+from typing import Union, List
+
+
+def calc_hemi_seaice(invar, tarea, hemi, extent=False):
+    """
+    Calculate hemispheric sea ice properties (area, volume, or extent).
+
+    This function computes sea ice properties over either the northern or southern
+    hemisphere. It can calculate:
+    - Sea ice area/volume: by multiplying ice concentration/thickness with grid cell area
+    - Sea ice extent: by summing grid cell areas where ice concentration ≥ 15%
+
+    Parameters
+    ----------
+    invar : xarray.DataArray
+        Input sea ice variable to process:
+        - For area calculation: sea ice concentration (aice, fractional 0-1)
+        - For volume calculation: sea ice thickness (hi, meters)
+        - For extent calculation: sea ice concentration (aice, fractional 0-1)
+    tarea : xarray.DataArray
+        Grid cell area (m²). Must have compatible dimensions with invar.
+    hemi : str
+        Hemisphere to calculate over:
+        - 'north': Northern hemisphere (latitude ≥ 0°)
+        - 'south': Southern hemisphere (latitude < 0°)
+    extent : bool, default False
+        Calculation mode:
+        - False: Calculate area/volume (invar × tarea)
+        - True: Calculate extent (sum tarea where invar ≥ 0.15)
+
+    Returns
+    -------
+    xarray.DataArray
+        Hemispheric sum of the requested property:
+        - Units: m² (for area/extent) or m³ (for volume)
+        - Dimensions: Time dimension preserved, spatial dimensions summed
+
+    Examples
+    --------
+    Calculate northern hemisphere sea ice extent:
+
+    >>> extent = calc_hemi_seaice(aice, areacello, 'north', extent=True)
+
+    Calculate southern hemisphere sea ice area:
+
+    >>> area = calc_hemi_seaice(aice, areacello, 'south', extent=False)
+
+    Calculate northern hemisphere sea ice volume:
+
+    >>> volume = calc_hemi_seaice(ice_thickness, areacello, 'north', extent=False)
+
+    Notes
+    -----
+    - Sea ice extent uses the standard 15% concentration threshold
+    - Function automatically detects latitude coordinates in input data
+    - Handles coordinate dimension alignment between invar and tarea
+    - Time dimension is preserved; spatial dimensions are summed
+    """
+    # Get latitude coordinates from input data
+    lat = None
+
+    # Try to find latitude coordinate in invar first
+    for coord_name in invar.coords:
+        if 'lat' in coord_name.lower():
+            lat = invar.coords[coord_name]
+            break
+
+    # If not found in invar, try tarea
+    if lat is None:
+        for coord_name in tarea.coords:
+            if 'lat' in coord_name.lower():
+                lat = tarea.coords[coord_name]
+                break
+
+    # If still not found, raise error
+    if lat is None:
+        raise ValueError("Cannot find latitude coordinate in input data")
+
+    # Ensure latitude coordinates align with data dimensions
+    invar_dims = invar.dims[1:] if invar.dims[0] == 'time' else invar.dims
+
+    # Align latitude coordinate with data dimensions if needed
+    if any(d not in invar_dims for d in lat.dims):
+        # Create mapping from lat dims to invar dims
+        lat_mapping = {}
+        for i, d in enumerate(lat.dims):
+            if i < len(invar_dims):
+                lat_mapping[d] = invar_dims[i]
+        lat = lat.rename(lat_mapping)
+
+    # if calculating extent sum carea and aice is used as filter
+    # with volume and area invar is multiplied by carea first
+    if extent:
+        var = tarea.where(invar >= 0.15).where(invar <= 1.)
+    else:
+        var = invar * tarea
+
+    if hemi == 'north':
+        var = var.where(lat >= 0.)
+    elif hemi == 'south':
+        var = var.where(lat < 0.)
+    else:
+        logging.error(f"invalid hemisphere: {hemi}")
+        raise ValueError(f"invalid hemisphere: {hemi}")
+
+    # sum over latitude and longitude (skip time dimension)
+    # This implements CMIP7 cell_methods: "area: sum time: mean" 
+    # (spatial sum, temporal mean)
+    spatial_dims = [d for d in var.dims if d != 'time']
+    vout = var.sum(dim=spatial_dims, skipna=True)
+
+    return vout
+
+
+def calc_seaice_extent(aice, areacello, region="north"):
+    """
+    Calculate hemispheric sea ice extent with CMIP-standard units.
+
+    This function calculates sea ice extent (total area of grid cells with ≥15% 
+    ice concentration) for a specified hemisphere and converts the result to 
+    CMIP-standard units (10⁶ km²).
+
+    Parameters
+    ----------
+    aice : xarray.DataArray
+        Sea ice concentration (fractional, range 0-1).
+        Must contain latitude coordinates for hemispheric selection.
+    areacello : xarray.DataArray
+        Ocean grid cell area in m².
+        Must have compatible dimensions with aice.
+    region : str, default 'north'
+        Hemisphere to calculate extent for:
+        - 'north': Northern hemisphere (latitude ≥ 0°)
+        - 'south': Southern hemisphere (latitude < 0°)
+
+    Returns
+    -------
+    xarray.DataArray
+        Sea ice extent in 10⁶ km² (CMIP standard units).
+        Time dimension preserved if present in input.
+
+    Examples
+    --------
+    Calculate Arctic sea ice extent:
+
+    >>> arctic_extent = calc_seaice_extent(aice, areacello, region="north")
+
+    Calculate Antarctic sea ice extent:
+
+    >>> antarctic_extent = calc_seaice_extent(aice, areacello, region="south")
+
+    Notes
+    -----
+    - Uses the standard 15% concentration threshold for extent calculation
+    - Output units are 10⁶ km² as required by CMIP protocols
+    - This is a convenience wrapper around calc_hemi_seaice
+    """
+    # Calculate extent using the general function
+    extent_m2 = calc_hemi_seaice(aice, areacello, region, extent=True)
+
+    # Convert from m^2 to 1e6 km^2 for CMIP standard units
+    extent_1e6km2 = extent_m2 / 1e12  # m^2 to 1e6 km^2
+
+    return extent_1e6km2