Cice grid

.github/workflows/ci.yml

@@ -29,7 +29,6 @@ jobs:
     strategy:
       matrix:
         python-version: ['3.10', '3.11', '3.12']
-
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python ${{ matrix.python-version }}
@@ -43,7 +42,7 @@ jobs:
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics  
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest

om3utils/cice_grid.py

@@ -0,0 +1,87 @@
+"""
+Script: cice_grid.py
+Description: 
+This script generates a CICE grid from the MOM super grid provided in the input NetCDF file.
+
+Usage:
+python cice_grid.py <ocean_hgrid> <ocean_hgrid>
+- ocean_hgrid: Path to the MOM super grid NetCDF file.
+- ocean_mask: Path to the corresponding mask NetCDF file.
+
+Dependencies: 
+- 'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10'.
+- or 'pyproject.toml' file
+
+"""
+
+#!/usr/bin/env python3
+# File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py
+
+import sys
+import argparse
+
+from esmgrids.mom_grid import MomGrid
+from esmgrids.cice_grid import CiceGrid
+
+
+class CiceGridNc:
+    """
+    Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc
+    """
+
+    def __init__(self, grid_file="grid.nc", mask_file="kmt.nc"):
+        self.grid_file = grid_file
+        self.mask_file = mask_file
+        return
+
+    def build_from_mom(self, ocean_hgrid, ocean_mask):
+        mom = MomGrid.fromfile(ocean_hgrid, mask_file=ocean_mask)
+
+        cice = CiceGrid.fromgrid(mom)
+
+        cice.create_gridnc(self.grid_file)
+
+        # Add versioning information
+        cice.grid_f.inputfile = f"{ocean_hgrid}"
+        cice.grid_f.inputfile_md5 = md5sum(ocean_hgrid)
+        cice.grid_f.history_command = f"python make_CICE_grid.py {ocean_hgrid} {ocean_mask}"
+
+        # Add the typical crs (i.e. WGS84/EPSG4326 , but in radians).
+        crs = cice.grid_f.createVariable("crs", "S1")
+        crs.grid_mapping_name = "tripolar_latitude_longitude"
+        crs.crs_wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["radians",1,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
+
+        cice.write()
+
+        cice.create_masknc(self.mask_file)
+
+        # Add versioning information
+        cice.mask_f.inputfile = f"{ocean_mask}"
+        cice.mask_f.inputfile_md5 = md5sum(ocean_mask)
+        cice.mask_f.history_command = f"python cice_grid.py {ocean_hgrid} {ocean_mask}"
+
+        # Add the typical crs (i.e. WGS84/EPSG4326 , but in radians).
+        crs = cice.mask_f.createVariable("crs", "S1")
+        crs.grid_mapping_name = "tripolar_latitude_longitude"
+        crs.crs_wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["radians",1,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
+
+        cice.write_mask()
+
+
+if __name__ == "__main__":
+    # command line arguments
+    from utils import md5sum  # load from file
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("ocean_hgrid", help="ocean_hgrid.nc file")
+    parser.add_argument("ocean_mask", help="ocean_mask.nc file")
+    # to-do: add argument for CRS & output filenames?
+
+    args = vars(parser.parse_args())
+
+    grid = CiceGridNc()
+
+    sys.exit(grid.build_from_mom(**args))
+
+else:
+    from .utils import md5sum  # load from package

om3utils/pbs_make_cice_grids.sh

@@ -0,0 +1,31 @@
+#!/bin/bash
+
+#PBS -W umask=0022
+#PBS -l mem=24gb
+#PBS -l storage=gdata/ik11+gdata/tm70+gdata/hh5
+#PBS -l wd
+#PBS -j oe
+
+umask 0003
+
+module purge
+module use /g/data/hh5/public/modules
+module load conda/analysis3-23.10
+
+echo "1 degree"
+python3 cice_grid.py /g/data/ik11/inputs/access-om2/input_20201102/mom_1deg/ocean_hgrid.nc /g/data/ik11/inputs/access-om2/input_20201102/mom_1deg/ocean_mask.nc
+
+mkdir 1deg
+mv grid.nc kmt.nc 1deg
+
+echo "0.25 deg"
+python3 cice_grid.py /g/data/ik11/inputs/access-om2/input_20230515_025deg_topog/mom_025deg/ocean_hgrid.nc /g/data/ik11/inputs/access-om2/input_20230515_025deg_topog/mom_025deg/ocean_mask.nc
+
+mkdir 025deg
+mv grid.nc kmt.nc 025deg
+
+echo "01 deg"
+python3 cice_grid.py /g/data/ik11/inputs/access-om2/input_20201102/mom_01deg/ocean_hgrid.nc /g/data/ik11/inputs/access-om2/input_20201102/mom_01deg/ocean_mask.nc
+
+mkdir 01deg
+mv grid.nc kmt.nc 01deg

om3utils/utils.py

@@ -0,0 +1,6 @@
+def md5sum(filename):
+    from hashlib import md5
+    from mmap import mmap, ACCESS_READ
+
+    with open(filename) as file, mmap(file.fileno(), 0, access=ACCESS_READ) as file:
+        return md5(file).hexdigest()

pyproject.toml

@@ -9,6 +9,11 @@ dynamic = ["version"]
 dependencies = [
     "f90nml",
     "ruamel.yaml",
+    "xarray",
+    "numpy",
+    "netcdf4",
+    "ocean_model_grid_generator@git+https://github.com/nikizadehgfdl/ocean_model_grid_generator@790069b31f9791864ccd514a2b8f53f385a0452e",
+    "esmgrids@git+https://github.com/anton-seaice/esmgrids@e5791f17a699627566687021e250aac32567ecd4",
 ]
 
 [build-system]

tests/test_cice_grid.py

@@ -0,0 +1,197 @@
+import pytest
+import xarray as xr
+from numpy.testing import assert_allclose
+from numpy import deg2rad
+from subprocess import run
+
+from om3utils.cice_grid import CiceGridNc
+
+# ----------------
+# test data:
+
+
+class MomGrid:
+    """Generate a sample tripole grid to use as test data"""
+
+    def __init__(self, tmp_path):
+        self.path = str(tmp_path) + "/ocean_hgrid.nc"
+        self.mask_path = str(tmp_path) + "/ocean_mask.nc"
+
+        # generate an tripolar grid as test data
+        run(
+            [
+                "ocean_grid_generator.py",
+                "-r",
+                "0.25",  # 4 degree grid
+                "--no_south_cap",
+                "--ensure_nj_even",
+                "-f",
+                self.path,
+            ]
+        )
+
+        # open ocean_hgrid.nc
+        self.ds = xr.open_dataset(self.path)
+
+        # an ocean mask with a arbitrary mask
+        self.mask_ds = xr.Dataset()
+        self.mask_ds["mask"] = (self.ds.area.coarsen(ny=2).sum().coarsen(nx=2).sum()) > 5e9
+        self.mask_ds.to_netcdf(self.mask_path)
+
+
+class CiceGrid:
+    """Make the CICE grid, using script under test"""
+
+    def __init__(self, mom_grid, tmp_path):
+        self.path = str(tmp_path) + "/grid.nc"
+        self.kmt_path = str(tmp_path) + "/kmt.nc"
+        cice_grid = CiceGridNc(self.path, self.kmt_path)
+        cice_grid.build_from_mom(mom_grid.path, mom_grid.mask_path)
+        self.ds = xr.open_dataset(self.path, decode_cf=False)
+        self.kmt_ds = xr.open_dataset(self.kmt_path, decode_cf=False)
+
+
+# pytest doesn't support class fixtures, so we need these two constructor funcs
+@pytest.fixture
+def mom_grid(tmp_path):
+    return MomGrid(tmp_path)
+
+
+@pytest.fixture
+def cice_grid(mom_grid, tmp_path):
+    return CiceGrid(mom_grid, tmp_path)
+
+
+@pytest.fixture
+def test_grid_ds(mom_grid):
+    # this generates the expected answers
+
+    ds = mom_grid.ds
+
+    test_grid = xr.Dataset()
+    # corners of the cice grid are NE corner
+    test_grid["ulat"] = deg2rad(ds.y.isel(nxp=slice(2, None, 2), nyp=slice(2, None, 2)))
+    test_grid["ulon"] = deg2rad(ds.x.isel(nxp=slice(2, None, 2), nyp=slice(2, None, 2)))
+
+    # centers of cice grid
+    test_grid["tlat"] = deg2rad(ds.y.isel(nxp=slice(1, None, 2), nyp=slice(1, None, 2)))
+    test_grid["tlon"] = deg2rad(ds.x.isel(nxp=slice(1, None, 2), nyp=slice(1, None, 2)))
+
+    # length of top edge of cells
+    test_grid["htn"] = ds.dx.isel(nyp=slice(2, None, 2)).coarsen(nx=2).sum() * 100
+    # length of right edge of cells
+    test_grid["hte"] = ds.dy.isel(nxp=slice(2, None, 2)).coarsen(ny=2).sum() * 100
+
+    # angle at u point
+    test_grid["angle"] = deg2rad(ds.angle_dx.isel(nyp=slice(2, None, 2), nxp=slice(2, None, 2)))
+    # angle a t points
+    test_grid["angleT"] = deg2rad(ds.angle_dx.isel(nyp=slice(1, None, 2), nxp=slice(1, None, 2)))
+
+    # area of cells
+    test_grid["tarea"] = mom_grid.ds.area.coarsen(ny=2).sum().coarsen(nx=2).sum()
+
+    # uarea is area of a cell centred around the u point
+    # we need to wrap in latitude and fold on longitude to calculate this
+    area_wrapped = mom_grid.ds.area
+    area_wrapped = xr.concat([ds.area.isel(nx=slice(1, None)), ds.area.isel(nx=0)], dim="nx")
+
+    top_row = xr.concat([ds.area.isel(ny=-1, nx=slice(-2, 0, -1)), ds.area.isel(ny=-1, nx=[-1, 0])], dim="nx")
+
+    area_folded = xr.concat([area_wrapped.isel(ny=slice(1, None)), top_row], dim="ny")
+
+    test_grid["uarea"] = area_folded.coarsen(ny=2).sum().coarsen(nx=2).sum()
+
+    return test_grid
+
+
+# ----------------
+# the tests in earnest:
+
+
+@pytest.mark.filterwarnings("ignore::DeprecationWarning")
+def test_cice_var_list(cice_grid, test_grid_ds):
+    # Test : Are there missing vars in cice_grid?
+    assert set(test_grid_ds.variables).difference(cice_grid.ds.variables) == set()
+
+
+@pytest.mark.filterwarnings("ignore::DeprecationWarning")
+def test_cice_grid(cice_grid, test_grid_ds):
+    # Test : Is the data the same as the test_grid
+    for jVar in test_grid_ds.variables:
+        assert_allclose(cice_grid.ds[jVar], test_grid_ds[jVar], rtol=1e-13, verbose=True, err_msg=f"{jVar} mismatch")
+
+
+def test_cice_kmt(mom_grid, cice_grid):
+    # Test : does the mask match
+    mask = mom_grid.mask_ds.mask
+    kmt = cice_grid.kmt_ds.kmt
+
+    assert_allclose(mask, kmt, rtol=1e-13, verbose=True, err_msg="mask mismatch")
+
+
+def test_cice_grid_attributes(cice_grid):
+    # Test: do the expected attributes to exist in the cice ds
+    cf_attributes = {
+        "ulat": {"standard_name": "latitude", "units": "radians"},
+        "ulon": {"standard_name": "longitude", "units": "radians"},
+        "tlat": {"standard_name": "latitude", "units": "radians"},
+        "tlon": {"standard_name": "longitude", "units": "radians"},
+        "uarea": {
+            "standard_name": "cell_area",
+            "units": "m^2",
+            "grid_mapping": "crs",
+            "coordinates": "ulat ulon",
+        },
+        "tarea": {
+            "standard_name": "cell_area",
+            "units": "m^2",
+            "grid_mapping": "crs",
+            "coordinates": "tlat tlon",
+        },
+        "angle": {
+            "standard_name": "angle_of_rotation_from_east_to_x",
+            "units": "radians",
+            "grid_mapping": "crs",
+            "coordinates": "ulat ulon",
+        },
+        "angleT": {
+            "standard_name": "angle_of_rotation_from_east_to_x",
+            "units": "radians",
+            "grid_mapping": "crs",
+            "coordinates": "tlat tlon",
+        },
+        "htn": {"units": "cm", "coordinates": "ulat tlon", "grid_mapping": "crs"},
+        "hte": {"units": "cm", "coordinates": "tlat ulon", "grid_mapping": "crs"},
+    }
+
+    for iVar in cf_attributes.keys():
+        print(cice_grid.ds[iVar])
+
+        for jAttr in cf_attributes[iVar].keys():
+            assert cice_grid.ds[iVar].attrs[jAttr] == cf_attributes[iVar][jAttr]
+
+
+def test_crs_exist(cice_grid):
+    # Test: has the crs been added ?
+    # todo: open with GDAL and rioxarray and confirm they find the crs?
+    assert hasattr(cice_grid.ds, "crs")
+    assert hasattr(cice_grid.kmt_ds, "crs")
+
+
+def test_inputs_logged(cice_grid, mom_grid):
+    # Test: have the source data been logged ?
+
+    for ds in [cice_grid.ds, cice_grid.kmt_ds]:
+        assert hasattr(ds, "inputfile"), "inputfile attribute missing"
+        assert hasattr(ds, "inputfile_md5"), "inputfile md5sum attribute missing"
+
+        sys_md5 = run(["md5sum", ds.inputfile], capture_output=True, text=True)
+        sys_md5 = sys_md5.stdout.split(" ")[0]
+        assert ds.inputfile_md5 == sys_md5, f"inputfile md5sum attribute incorrect, {ds.inputfile_md5} != {sys_md5}"
+
+    assert (
+        cice_grid.ds.inputfile == mom_grid.path
+    ), "inputfile attribute incorrect ({cice_grid.ds.inputfile} != {mom_grid.path})"
+    assert (
+        cice_grid.kmt_ds.inputfile == mom_grid.mask_path
+    ), "mask inputfile attribute incorrect ({cice_grid.kmt_ds.inputfile} != {mom_grid.mask_path})"

tests/utils.py

@@ -1,5 +1,5 @@
 class MockFile:
-    """Class for testing parsers that require a file.
+    """Class for testing parsers that require a text file.
 
     Usage:
     @pytest.fixture