diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 748358a..9496b0d 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/README.md b/README.md index 472929e..622947c 100644 --- a/README.md +++ b/README.md @@ -7,3 +7,4 @@ Collection of utilities aimed at simplifying the creation and handling of ACCESS-OM3 runs. It currently includes: - functions to read and write ACCESS-OM3 configuration files - functions to read and process profiling data + - scripts to create the CICE grid file from an existing MOM grid file diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py new file mode 100644 index 0000000..91d8ed3 --- /dev/null +++ b/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: 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' + pip install esmgrids +- 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 diff --git a/om3utils/utils.py b/om3utils/utils.py new file mode 100644 index 0000000..157791e --- /dev/null +++ b/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() diff --git a/pyproject.toml b/pyproject.toml index c4d4334..c7f4b43 100644 --- a/pyproject.toml +++ b/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] diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 0000000..01f3ef0 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,5 @@ +# Scripts + +This folder contains scripts which are the reference implementation of how to use software within OM3 utils. + +i.e. run `qsub pbs_make_cice_grids.sh` to create the cice grid for ACCESS-OM3 at all three standard resolutions. After sanity checking, these are then moved by hand to the desired folder to run the model. \ No newline at end of file diff --git a/scripts/pbs_make_cice_grids.sh b/scripts/pbs_make_cice_grids.sh new file mode 100644 index 0000000..d5e46f8 --- /dev/null +++ b/scripts/pbs_make_cice_grids.sh @@ -0,0 +1,33 @@ +#!/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 + +run='python3 ../om3utils/cice_grid.py' + +umask 0003 + +module purge +module use /g/data/hh5/public/modules +module load conda/analysis3-23.10 + +echo "1 degree" +$run /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" +$run /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" +$run /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 \ No newline at end of file diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py new file mode 100644 index 0000000..43acb80 --- /dev/null +++ b/tests/test_cice_grid.py @@ -0,0 +1,199 @@ +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 a 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 + # In simple terms the MOM supergrid has four cells for each model grid cell. The MOM supergrid includes all edges (left and right) but CICE only uses right/east edges. (e.g. For points/edges of first cell: 0,0 is SW corner, 1,1 is middle of cell, 2,2, is NE corner/edges) + + ds = mom_grid.ds + + # u points are at top-right (NE) corner + u_points = ds.isel(nxp=slice(2, None, 2), nyp=slice(2, None, 2)) + + # t points are in middle of cell + t_points = ds.isel(nxp=slice(1, None, 2), nyp=slice(1, None, 2)) + + test_grid = xr.Dataset() + + test_grid["ulat"] = deg2rad(u_points.y) + test_grid["ulon"] = deg2rad(u_points.x) + test_grid["tlat"] = deg2rad(t_points.y) + test_grid["tlon"] = deg2rad(t_points.x) + + test_grid["angle"] = deg2rad(u_points.angle_dx) # angle at u point + test_grid["angleT"] = deg2rad(t_points.angle_dx) + + # length of top (northern) edge of cells + test_grid["htn"] = ds.dx.isel(nyp=slice(2, None, 2)).coarsen(nx=2).sum() * 100 + # length of right (eastern) edge of cells + test_grid["hte"] = ds.dy.isel(nxp=slice(2, None, 2)).coarsen(ny=2).sum() * 100 + + # area of cells + test_grid["tarea"] = ds.area.coarsen(ny=2).sum().coarsen(nx=2).sum() + + # uarea is area of a cell centred around the u point + # we need to fold on longitude and wrap in latitude to calculate this + # drop the bottom row, new top row is reverse of current top row + area_folded = xr.concat([ds.area.isel(ny=slice(1, None)), ds.area.isel(ny=-1, nx=slice(-1, None, -1))], dim="ny") + + # drop the first column, make the new last column the first column + area_wrapped = xr.concat([area_folded.isel(nx=slice(1, None)), area_folded.isel(nx=0)], dim="nx") + + test_grid["uarea"] = area_wrapped.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})" diff --git a/tests/utils.py b/tests/utils.py index 1d2fb1c..c3a8246 100644 --- a/tests/utils.py +++ b/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