From a043e008a3c2709664a3523df68e9634c6559c29 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Tue, 2 Apr 2024 10:59:09 +1100 Subject: [PATCH 01/25] commit before refactor --- .gitmodules | 6 ++ grid_generators/cice_grid.py | 77 ++++++++++++++++++++++ grid_generators/esmgrids | 1 + grid_generators/ocean_model_grid_generator | 1 + grid_generators/pbs_make_cice_grids.sh | 31 +++++++++ om3utils/utils.py | 6 ++ pyproject.toml | 3 + 7 files changed, 125 insertions(+) create mode 100644 .gitmodules create mode 100644 grid_generators/cice_grid.py create mode 160000 grid_generators/esmgrids create mode 160000 grid_generators/ocean_model_grid_generator create mode 100644 grid_generators/pbs_make_cice_grids.sh create mode 100644 om3utils/utils.py diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..be93f98 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "grid_generators/esmgrids"] + path = grid_generators/esmgrids + url = https://github.com/COSIMA/esmgrids.git +[submodule "grid_generators/ocean_model_grid_generator"] + path = grid_generators/ocean_model_grid_generator + url = https://github.com/nikizadehgfdl/ocean_model_grid_generator.git diff --git a/grid_generators/cice_grid.py b/grid_generators/cice_grid.py new file mode 100644 index 0000000..5326133 --- /dev/null +++ b/grid_generators/cice_grid.py @@ -0,0 +1,77 @@ +""" +Script: make_cice_grid.py +Description: +This script generates a CICE grid from the MOM super grid provided in the input NetCDF file. + +Usage: +python make_cice_grid.py +- ocean_hgrid: Path to the MOM super grid NetCDF file. +- ocean_mask: Path to the corresponding mask NetCDF file. + +""" + + +#!/usr/bin/env python3 +# File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py + +import sys +import os +import argparse + +my_dir = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(os.path.join(my_dir, 'esmgrids')) +sys.path.append(os.path.join(my_dir, '../')) + +print(sys.path) + +from esmgrids.mom_grid import MomGrid # noqa +from esmgrids.cice_grid import CiceGrid # noqa +from om3utils.utils import md5sum + + +class cice_grid_from_mom() : + + """ + Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc + """ + + def run(ocean_hgrid, ocean_mask): + + + mom = MomGrid.fromfile(ocean_hgrid, mask_file=ocean_mask) + + cice = CiceGrid.fromgrid(mom) + + + grid_file = os.path.join('grid.nc') + mask_file = os.path.join('kmt.nc') + + cice.create_gridnc(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}" + + cice.write() + + cice.create_masknc(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 make_CICE_grid.py {ocean_hgrid} {ocean_mask}" + + cice.write_mask() + +if __name__ == '__main__': + + parser = argparse.ArgumentParser() + parser.add_argument('ocean_hgrid', help='ocean_hgrid.nc file') + parser.add_argument('ocean_mask', help='ocean_mask.nc file') + #add argument for CRS? + + args = vars(parser.parse_args()) + + + sys.exit(main(**args)) diff --git a/grid_generators/esmgrids b/grid_generators/esmgrids new file mode 160000 index 0000000..514b560 --- /dev/null +++ b/grid_generators/esmgrids @@ -0,0 +1 @@ +Subproject commit 514b560adf7b9094081408da49ecb3e19ef0353b diff --git a/grid_generators/ocean_model_grid_generator b/grid_generators/ocean_model_grid_generator new file mode 160000 index 0000000..790069b --- /dev/null +++ b/grid_generators/ocean_model_grid_generator @@ -0,0 +1 @@ +Subproject commit 790069b31f9791864ccd514a2b8f53f385a0452e diff --git a/grid_generators/pbs_make_cice_grids.sh b/grid_generators/pbs_make_cice_grids.sh new file mode 100644 index 0000000..8aa30a5 --- /dev/null +++ b/grid_generators/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 make_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 make_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 make_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 \ No newline at end of file diff --git a/om3utils/utils.py b/om3utils/utils.py new file mode 100644 index 0000000..c068699 --- /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() \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index c4d4334..62b0e75 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,6 +9,9 @@ dynamic = ["version"] dependencies = [ "f90nml", "ruamel.yaml", + "numpypi", + "xarray", + "numpy" ] [build-system] From fe4c53870c3349fd6999f5f08336f406fd71a654 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Tue, 2 Apr 2024 16:36:44 +1100 Subject: [PATCH 02/25] Working pytest suite for cice grid --- .gitmodules | 8 +- grid_generators/cice_grid.py | 77 -------- grid_generators/esmgrids | 1 - om3utils/cice_grid.py | 85 +++++++++ om3utils/esmgrids | 1 + .../pbs_make_cice_grids.sh | 0 om3utils/utils.py | 14 +- .../ocean_model_grid_generator | 0 tests/test_cice_grid.py | 177 ++++++++++++++++++ tests/utils.py | 2 +- 10 files changed, 276 insertions(+), 89 deletions(-) delete mode 100644 grid_generators/cice_grid.py delete mode 160000 grid_generators/esmgrids create mode 100644 om3utils/cice_grid.py create mode 160000 om3utils/esmgrids rename {grid_generators => om3utils}/pbs_make_cice_grids.sh (100%) rename {grid_generators => tests}/ocean_model_grid_generator (100%) create mode 100644 tests/test_cice_grid.py diff --git a/.gitmodules b/.gitmodules index be93f98..f595200 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,6 +1,6 @@ -[submodule "grid_generators/esmgrids"] - path = grid_generators/esmgrids +[submodule "om3utils/esmgrids"] + path = om3utils/esmgrids url = https://github.com/COSIMA/esmgrids.git -[submodule "grid_generators/ocean_model_grid_generator"] - path = grid_generators/ocean_model_grid_generator +[submodule "tests/ocean_model_grid_generator"] + path = tests/ocean_model_grid_generator url = https://github.com/nikizadehgfdl/ocean_model_grid_generator.git diff --git a/grid_generators/cice_grid.py b/grid_generators/cice_grid.py deleted file mode 100644 index 5326133..0000000 --- a/grid_generators/cice_grid.py +++ /dev/null @@ -1,77 +0,0 @@ -""" -Script: make_cice_grid.py -Description: -This script generates a CICE grid from the MOM super grid provided in the input NetCDF file. - -Usage: -python make_cice_grid.py -- ocean_hgrid: Path to the MOM super grid NetCDF file. -- ocean_mask: Path to the corresponding mask NetCDF file. - -""" - - -#!/usr/bin/env python3 -# File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py - -import sys -import os -import argparse - -my_dir = os.path.dirname(os.path.realpath(__file__)) -sys.path.append(os.path.join(my_dir, 'esmgrids')) -sys.path.append(os.path.join(my_dir, '../')) - -print(sys.path) - -from esmgrids.mom_grid import MomGrid # noqa -from esmgrids.cice_grid import CiceGrid # noqa -from om3utils.utils import md5sum - - -class cice_grid_from_mom() : - - """ - Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc - """ - - def run(ocean_hgrid, ocean_mask): - - - mom = MomGrid.fromfile(ocean_hgrid, mask_file=ocean_mask) - - cice = CiceGrid.fromgrid(mom) - - - grid_file = os.path.join('grid.nc') - mask_file = os.path.join('kmt.nc') - - cice.create_gridnc(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}" - - cice.write() - - cice.create_masknc(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 make_CICE_grid.py {ocean_hgrid} {ocean_mask}" - - cice.write_mask() - -if __name__ == '__main__': - - parser = argparse.ArgumentParser() - parser.add_argument('ocean_hgrid', help='ocean_hgrid.nc file') - parser.add_argument('ocean_mask', help='ocean_mask.nc file') - #add argument for CRS? - - args = vars(parser.parse_args()) - - - sys.exit(main(**args)) diff --git a/grid_generators/esmgrids b/grid_generators/esmgrids deleted file mode 160000 index 514b560..0000000 --- a/grid_generators/esmgrids +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 514b560adf7b9094081408da49ecb3e19ef0353b diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py new file mode 100644 index 0000000..e6ccf8e --- /dev/null +++ b/om3utils/cice_grid.py @@ -0,0 +1,85 @@ +""" +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. + +""" + + +#!/usr/bin/env python3 +# File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py + +import sys +# import os +import argparse + +# my_dir = os.path.dirname(os.path.realpath(__file__)) +# sys.path.append(os.path.join(my_dir, 'esmgrids')) + +# print(sys.path) + +from .esmgrids.esmgrids.mom_grid import MomGrid # noqa +from .esmgrids.esmgrids.cice_grid import CiceGrid # noqa +from .utils import Md5sum + +class cice_grid_from_mom() : + + """ + Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc + """ + + def run(self, ocean_hgrid, ocean_mask): + + + mom = MomGrid.fromfile(ocean_hgrid, mask_file=ocean_mask) + + cice = CiceGrid.fromgrid(mom) + + # grid_file = os.path.join('grid.nc') + # mask_file = os.path.join('kmt.nc') + + cice.create_gridnc('grid.nc') + + # Add versioning information + cice.grid_f.inputfile = f"{ocean_hgrid}" + cice.grid_f.inputfile_md5 = Md5sum(ocean_hgrid).sum + 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('kmt.nc') + + # Add versioning information + cice.mask_f.inputfile = f"{ocean_mask}" + cice.mask_f.inputfile_md5 = Md5sum(ocean_mask).sum + cice.mask_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.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__': + + parser = argparse.ArgumentParser() + parser.add_argument('ocean_hgrid', help='ocean_hgrid.nc file') + parser.add_argument('ocean_mask', help='ocean_mask.nc file') + #add argument for CRS? + + args = vars(parser.parse_args()) + + grid = cice_grid_from_mom() + + sys.exit(grid.run(**args)) diff --git a/om3utils/esmgrids b/om3utils/esmgrids new file mode 160000 index 0000000..5ca8760 --- /dev/null +++ b/om3utils/esmgrids @@ -0,0 +1 @@ +Subproject commit 5ca87608a34464521729bfb87607b0931bfe6da0 diff --git a/grid_generators/pbs_make_cice_grids.sh b/om3utils/pbs_make_cice_grids.sh similarity index 100% rename from grid_generators/pbs_make_cice_grids.sh rename to om3utils/pbs_make_cice_grids.sh diff --git a/om3utils/utils.py b/om3utils/utils.py index c068699..69814fe 100644 --- a/om3utils/utils.py +++ b/om3utils/utils.py @@ -1,6 +1,8 @@ -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() \ No newline at end of file +class Md5sum: + + def __init__(self,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: + self.sum = md5(file).hexdigest() diff --git a/grid_generators/ocean_model_grid_generator b/tests/ocean_model_grid_generator similarity index 100% rename from grid_generators/ocean_model_grid_generator rename to tests/ocean_model_grid_generator diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py new file mode 100644 index 0000000..9e8630c --- /dev/null +++ b/tests/test_cice_grid.py @@ -0,0 +1,177 @@ +import pytest +import xarray as xr +from numpy.testing import assert_allclose +from numpy import deg2rad + +#im not sure how micael is importing om3utils? +import sys +import os +my_dir = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(os.path.join(my_dir, '..')) + +from om3utils.cice_grid import cice_grid_from_mom +from ocean_model_grid_generator.ocean_grid_generator import main as ocean_grid_generator + +import warnings + +class MomGrid: + + """Generate a sample tripole grid to use as test data""" + + path = 'ocean_hgrid.nc' + mask_path = 'ocean_mask.nc' + + def __init__(self): + # generate an tripolar grid as test data + args = { + 'inverse_resolution': 0.25 , #4 degree grid + 'no_south_cap': True, + 'ensure_nj_even': True, + 'match_dy': ["bp", "so", "p125sc", ""], + 'gridfilename': 'ocean_hgrid.nc' + } + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category = RuntimeWarning) + ocean_grid_generator(**args) #generates ocean_hgrid.nc + + self.ds = xr.open_dataset(self.path) + + # an ocean mask with a random 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""" + + path = 'grid.nc' + kmt_path = 'kmt.nc' + def __init__(self, mom_grid): + # the data under test + cice_grid = cice_grid_from_mom() + cice_grid.run(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) + +def make_test_grid(mom_grid_ds): + # this generates the expected answers + + test_grid = xr.Dataset() + # corners of the cice grid are NE corner + test_grid['ulat']=deg2rad(mom_grid_ds.y.isel(nxp=slice(2,None,2), nyp=slice(2,None,2))) + test_grid['ulon']=deg2rad(mom_grid_ds.x.isel(nxp=slice(2,None,2), nyp=slice(2,None,2))) + + # centers of cice grid + test_grid['tlat']=deg2rad(mom_grid_ds.y.isel(nxp=slice(1,None,2), nyp=slice(1,None,2))) + test_grid['tlon']=deg2rad(mom_grid_ds.x.isel(nxp=slice(1,None,2), nyp=slice(1,None,2))) + + # length of top edge of cells + test_grid['htn']=(mom_grid_ds.dx.isel(nyp=slice(2,None,2)).coarsen(nx=2).sum()*100) + # length of right edge of cells + test_grid['hte']=mom_grid_ds.dy.isel(nxp=slice(2,None,2)).coarsen(ny=2).sum()*100 + + # angle at u point + test_grid['angle']=deg2rad(mom_grid_ds.angle_dx.isel(nyp=slice(2,None,2), nxp=slice(2,None,2))) + # angle a t points + test_grid['angleT']=deg2rad(mom_grid_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([ + mom_grid_ds.area.isel(nx=slice(1,None)), + mom_grid_ds.area.isel(nx=0) + ], dim='nx') + + top_row = xr.concat([ + mom_grid_ds.area.isel(ny=-1, nx=slice(-2,0,-1)), + mom_grid_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 + +@pytest.fixture +def mom_grid(): + return MomGrid() + +@pytest.fixture +def cice_grid(mom_grid): + return CiceGrid(mom_grid) + +@pytest.mark.filterwarnings('ignore::DeprecationWarning') +def test_cice_grid(mom_grid, cice_grid): + #to-do: run at high res? + + test_grid_ds = make_test_grid(mom_grid.ds) + + #Test1 : Are there missing vars in cice_grid? + assert(set(test_grid_ds.variables).difference(cice_grid.ds.variables) == set() ) + + #Test2 : Is the data correct + for jVar in test_grid_ds.variables: + + anom = (cice_grid.ds[jVar].values-test_grid_ds[jVar].values) + + print(f'{jVar} anom min: {anom.min()}, anom max: {anom.max()}') + + assert_allclose( + cice_grid.ds[jVar], + test_grid_ds[jVar], + rtol=1e-13, + verbose=True, + err_msg=f'{jVar} mismatch' + ) + + pass + +def test_cice_kmt(mom_grid, cice_grid): + mask = mom_grid.mask_ds.mask + kmt = cice_grid.kmt_ds.kmt + + assert_allclose( + mask, + kmt, + rtol=1e-13, + verbose=True, + err_msg=f'mask mismatch' + ) + + pass + +def test_cice_grid_attributes(cice_grid): + + #expected attributes to exist in the 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): + # todo: open with GDAL and rioxarray and confirm they find the crs? + assert hasattr(cice_grid.ds, 'crs') 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 From 4a0477e622631c768958a6de49d90c5838a60f80 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 14:27:21 +1100 Subject: [PATCH 03/25] updated tests --- om3utils/cice_grid.py | 54 +++++++++------- om3utils/esmgrids | 2 +- om3utils/pbs_make_cice_grids.sh | 6 +- om3utils/utils.py | 6 +- pyproject.toml | 1 + tests/test_cice_grid.py | 108 +++++++++++++++++++------------- 6 files changed, 101 insertions(+), 76 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index e6ccf8e..5d11111 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -8,46 +8,47 @@ - ocean_hgrid: Path to the MOM super grid NetCDF file. - ocean_mask: Path to the corresponding mask NetCDF file. -""" +Dependencies: Developed using +'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10' +""" #!/usr/bin/env python3 # File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py import sys -# import os +import os import argparse -# my_dir = os.path.dirname(os.path.realpath(__file__)) -# sys.path.append(os.path.join(my_dir, 'esmgrids')) - -# print(sys.path) +my_dir = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(os.path.join(my_dir, 'esmgrids')) -from .esmgrids.esmgrids.mom_grid import MomGrid # noqa -from .esmgrids.esmgrids.cice_grid import CiceGrid # noqa -from .utils import Md5sum +from esmgrids.mom_grid import MomGrid +from esmgrids.cice_grid import CiceGrid -class cice_grid_from_mom() : +class cice_grid_nc: """ Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc """ - def run(self, ocean_hgrid, ocean_mask): + 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) - - # grid_file = os.path.join('grid.nc') - # mask_file = os.path.join('kmt.nc') - - cice.create_gridnc('grid.nc') + + cice.create_gridnc(self.grid_file) # Add versioning information cice.grid_f.inputfile = f"{ocean_hgrid}" - cice.grid_f.inputfile_md5 = Md5sum(ocean_hgrid).sum + 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). @@ -57,12 +58,12 @@ def run(self, ocean_hgrid, ocean_mask): cice.write() - cice.create_masknc('kmt.nc') + cice.create_masknc(self.mask_file) # Add versioning information cice.mask_f.inputfile = f"{ocean_mask}" - cice.mask_f.inputfile_md5 = Md5sum(ocean_mask).sum - cice.mask_f.history_command = f"python make_CICE_grid.py {ocean_hgrid} {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') @@ -73,13 +74,20 @@ def run(self, ocean_hgrid, ocean_mask): if __name__ == '__main__': + from utils import md5sum + parser = argparse.ArgumentParser() parser.add_argument('ocean_hgrid', help='ocean_hgrid.nc file') parser.add_argument('ocean_mask', help='ocean_mask.nc file') - #add argument for CRS? + #to-do: add argument for CRS? args = vars(parser.parse_args()) - grid = cice_grid_from_mom() + grid = cice_grid_nc() + + sys.exit(grid.build_from_mom(**args)) + +else: + #for testing + from .utils import md5sum - sys.exit(grid.run(**args)) diff --git a/om3utils/esmgrids b/om3utils/esmgrids index 5ca8760..fb79a79 160000 --- a/om3utils/esmgrids +++ b/om3utils/esmgrids @@ -1 +1 @@ -Subproject commit 5ca87608a34464521729bfb87607b0931bfe6da0 +Subproject commit fb79a794b6e0b2ee62e181d6cc336a381669fc50 diff --git a/om3utils/pbs_make_cice_grids.sh b/om3utils/pbs_make_cice_grids.sh index 8aa30a5..6655a38 100644 --- a/om3utils/pbs_make_cice_grids.sh +++ b/om3utils/pbs_make_cice_grids.sh @@ -13,19 +13,19 @@ module use /g/data/hh5/public/modules module load conda/analysis3-23.10 echo "1 degree" -python3 make_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 +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 make_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 +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 make_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 +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 \ No newline at end of file diff --git a/om3utils/utils.py b/om3utils/utils.py index 69814fe..ac39867 100644 --- a/om3utils/utils.py +++ b/om3utils/utils.py @@ -1,8 +1,6 @@ -class Md5sum: - - def __init__(self,filename): +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: - self.sum = md5(file).hexdigest() + return md5(file).hexdigest() diff --git a/pyproject.toml b/pyproject.toml index 62b0e75..722d48c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,7 @@ test = [ [tool.pytest.ini_options] addopts = ["--cov=om3utils", "--cov-report=term", "--cov-report=xml"] testpaths = ["tests"] +norecursedirs = ["tests/ocean_model_grid_generator"] [tool.coverage.run] omit = ["om3utils/__init__.py", "om3utils/_version.py"] diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 9e8630c..135e1c6 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -9,26 +9,30 @@ my_dir = os.path.dirname(os.path.realpath(__file__)) sys.path.append(os.path.join(my_dir, '..')) -from om3utils.cice_grid import cice_grid_from_mom +from om3utils.cice_grid import cice_grid_nc from ocean_model_grid_generator.ocean_grid_generator import main as ocean_grid_generator import warnings +# ---------------- +# test data: + class MomGrid: """Generate a sample tripole grid to use as test data""" - path = 'ocean_hgrid.nc' - mask_path = 'ocean_mask.nc' + def __init__(self, tmp_path): + + self.path = str(tmp_path)+'/ocean_hgrid.nc' + self.mask_path = str(tmp_path)+'/ocean_mask.nc' - def __init__(self): # generate an tripolar grid as test data args = { 'inverse_resolution': 0.25 , #4 degree grid 'no_south_cap': True, 'ensure_nj_even': True, 'match_dy': ["bp", "so", "p125sc", ""], - 'gridfilename': 'ocean_hgrid.nc' + 'gridfilename': self.path } with warnings.catch_warnings(): @@ -45,51 +49,53 @@ def __init__(self): class CiceGrid: """Make the CICE grid, using script under test""" - path = 'grid.nc' - kmt_path = 'kmt.nc' - def __init__(self, mom_grid): - # the data under test - cice_grid = cice_grid_from_mom() - cice_grid.run(mom_grid.path, mom_grid.mask_path) + def __init__(self, mom_grid, tmp_path): + self.path = str(tmp_path)+'/grid.nc' + self.kmt_path = str(tmp_path)+'/kmt.nc' + cice_grid = cice_grid_nc(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) -def make_test_grid(mom_grid_ds): +@pytest.fixture +def test_grid_ds(mom_grid, tmp_path): # 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(mom_grid_ds.y.isel(nxp=slice(2,None,2), nyp=slice(2,None,2))) - test_grid['ulon']=deg2rad(mom_grid_ds.x.isel(nxp=slice(2,None,2), nyp=slice(2,None,2))) + 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(mom_grid_ds.y.isel(nxp=slice(1,None,2), nyp=slice(1,None,2))) - test_grid['tlon']=deg2rad(mom_grid_ds.x.isel(nxp=slice(1,None,2), nyp=slice(1,None,2))) + 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']=(mom_grid_ds.dx.isel(nyp=slice(2,None,2)).coarsen(nx=2).sum()*100) + 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']=mom_grid_ds.dy.isel(nxp=slice(2,None,2)).coarsen(ny=2).sum()*100 + test_grid['hte']=ds.dy.isel(nxp=slice(2,None,2)).coarsen(ny=2).sum()*100 # angle at u point - test_grid['angle']=deg2rad(mom_grid_ds.angle_dx.isel(nyp=slice(2,None,2), nxp=slice(2,None,2))) + 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(mom_grid_ds.angle_dx.isel(nyp=slice(1,None,2), nxp=slice(1,None,2))) + 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() + 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 = mom_grid.ds.area area_wrapped = xr.concat([ - mom_grid_ds.area.isel(nx=slice(1,None)), - mom_grid_ds.area.isel(nx=0) + ds.area.isel(nx=slice(1,None)), + ds.area.isel(nx=0) ], dim='nx') top_row = xr.concat([ - mom_grid_ds.area.isel(ny=-1, nx=slice(-2,0,-1)), - mom_grid_ds.area.isel(ny=-1, nx=[-1,0]) + ds.area.isel(ny=-1, nx=slice(-2,0,-1)), + ds.area.isel(ny=-1, nx=[-1,0]) ], dim='nx') area_folded = xr.concat([ @@ -101,24 +107,25 @@ def make_test_grid(mom_grid_ds): return test_grid +#pytest doesn't support class fixtures, so we need these two constructor funcs @pytest.fixture -def mom_grid(): - return MomGrid() - +def mom_grid(tmp_path): + return MomGrid(tmp_path) @pytest.fixture -def cice_grid(mom_grid): - return CiceGrid(mom_grid) +def cice_grid(mom_grid, tmp_path): + return CiceGrid(mom_grid, tmp_path) -@pytest.mark.filterwarnings('ignore::DeprecationWarning') -def test_cice_grid(mom_grid, cice_grid): - #to-do: run at high res? - - test_grid_ds = make_test_grid(mom_grid.ds) +# ---------------- +# the tests in earnest: - #Test1 : Are there missing vars in cice_grid? +@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() ) - #Test2 : Is the data correct +@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: anom = (cice_grid.ds[jVar].values-test_grid_ds[jVar].values) @@ -133,9 +140,8 @@ def test_cice_grid(mom_grid, cice_grid): err_msg=f'{jVar} mismatch' ) - pass - 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 @@ -146,12 +152,9 @@ def test_cice_kmt(mom_grid, cice_grid): verbose=True, err_msg=f'mask mismatch' ) - - pass def test_cice_grid_attributes(cice_grid): - - #expected attributes to exist in the ds + #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'}, @@ -171,7 +174,22 @@ def test_cice_grid_attributes(cice_grid): 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 = os.popen(f'md5sum {ds.inputfile} | cut -f 1 -d " "').read().strip() + 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})" + assert (cice_grid.kmt_ds.inputfile == mom_grid.mask_path), "mask inputfile attribute incorrect ({cice_grid.kmt_ds.inputfile})" \ No newline at end of file From 36b70e1c99cce4988a4ea65706717ff6c5e82bda Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 14:32:01 +1100 Subject: [PATCH 04/25] black --- om3utils/cice_grid.py | 36 ++++---- om3utils/utils.py | 10 +-- tests/test_cice_grid.py | 186 +++++++++++++++++++++------------------- 3 files changed, 121 insertions(+), 111 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 5d11111..617dddf 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -21,38 +21,37 @@ import argparse my_dir = os.path.dirname(os.path.realpath(__file__)) -sys.path.append(os.path.join(my_dir, 'esmgrids')) +sys.path.append(os.path.join(my_dir, "esmgrids")) from esmgrids.mom_grid import MomGrid from esmgrids.cice_grid import CiceGrid + class cice_grid_nc: """ 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 + 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 + # 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') + # 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"]]' @@ -60,26 +59,26 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): cice.create_masknc(self.mask_file) - # Add versioning information + # 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') + # 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__': +if __name__ == "__main__": from utils import md5sum 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? + 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? args = vars(parser.parse_args()) @@ -88,6 +87,5 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): sys.exit(grid.build_from_mom(**args)) else: - #for testing + # for testing from .utils import md5sum - diff --git a/om3utils/utils.py b/om3utils/utils.py index ac39867..157791e 100644 --- a/om3utils/utils.py +++ b/om3utils/utils.py @@ -1,6 +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() + 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/tests/test_cice_grid.py b/tests/test_cice_grid.py index 135e1c6..96473ea 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -3,60 +3,63 @@ from numpy.testing import assert_allclose from numpy import deg2rad -#im not sure how micael is importing om3utils? +# im not sure how micael is importing om3utils? import sys import os + my_dir = os.path.dirname(os.path.realpath(__file__)) -sys.path.append(os.path.join(my_dir, '..')) +sys.path.append(os.path.join(my_dir, "..")) from om3utils.cice_grid import cice_grid_nc from ocean_model_grid_generator.ocean_grid_generator import main as ocean_grid_generator -import warnings +import warnings # ---------------- # 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' + self.path = str(tmp_path) + "/ocean_hgrid.nc" + self.mask_path = str(tmp_path) + "/ocean_mask.nc" # generate an tripolar grid as test data args = { - 'inverse_resolution': 0.25 , #4 degree grid - 'no_south_cap': True, - 'ensure_nj_even': True, - 'match_dy': ["bp", "so", "p125sc", ""], - 'gridfilename': self.path + "inverse_resolution": 0.25, # 4 degree grid + "no_south_cap": True, + "ensure_nj_even": True, + "match_dy": ["bp", "so", "p125sc", ""], + "gridfilename": self.path, } with warnings.catch_warnings(): - warnings.simplefilter("ignore", category = RuntimeWarning) - ocean_grid_generator(**args) #generates ocean_hgrid.nc + warnings.simplefilter("ignore", category=RuntimeWarning) + ocean_grid_generator(**args) # generates ocean_hgrid.nc self.ds = xr.open_dataset(self.path) - # an ocean mask with a random mask + # an ocean mask with a random 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["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 = cice_grid_nc(self.path, self.kmt_path) + self.path = str(tmp_path) + "/grid.nc" + self.kmt_path = str(tmp_path) + "/kmt.nc" + cice_grid = cice_grid_nc(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.ds = xr.open_dataset(self.path, decode_cf=False) self.kmt_ds = xr.open_dataset(self.kmt_path, decode_cf=False) + @pytest.fixture def test_grid_ds(mom_grid, tmp_path): # this generates the expected answers @@ -65,131 +68,140 @@ def test_grid_ds(mom_grid, tmp_path): 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))) + 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))) + 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) + 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 + 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))) + 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))) + 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() + 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') + 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') + 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') + 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() + test_grid["uarea"] = area_folded.coarsen(ny=2).sum().coarsen(nx=2).sum() return test_grid -#pytest doesn't support class fixtures, so we need these two constructor funcs + +# 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) + # ---------------- # the tests in earnest: -@pytest.mark.filterwarnings('ignore::DeprecationWarning') + +@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() ) + # 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') + +@pytest.mark.filterwarnings("ignore::DeprecationWarning") def test_cice_grid(cice_grid, test_grid_ds): - #Test : Is the data the same as the test_grid + # Test : Is the data the same as the test_grid for jVar in test_grid_ds.variables: - - anom = (cice_grid.ds[jVar].values-test_grid_ds[jVar].values) + anom = cice_grid.ds[jVar].values - test_grid_ds[jVar].values + + print(f"{jVar} anom min: {anom.min()}, anom max: {anom.max()}") - print(f'{jVar} anom min: {anom.min()}, anom max: {anom.max()}') + assert_allclose(cice_grid.ds[jVar], test_grid_ds[jVar], rtol=1e-13, verbose=True, err_msg=f"{jVar} mismatch") - 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 + # 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=f'mask mismatch' - ) + assert_allclose(mask, kmt, rtol=1e-13, verbose=True, err_msg=f"mask mismatch") + def test_cice_grid_attributes(cice_grid): - #Test: do the expected attributes to exist in the cice ds + # 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'} + "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 ? + # 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') + 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 ? + # 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" + 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 = os.popen(f'md5sum {ds.inputfile} | cut -f 1 -d " "').read().strip() - assert (ds.inputfile_md5 == sys_md5), f"inputfile md5sum attribute incorrect, {ds.inputfile_md5} != {sys_md5}" + 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})" - assert (cice_grid.kmt_ds.inputfile == mom_grid.mask_path), "mask inputfile attribute incorrect ({cice_grid.kmt_ds.inputfile})" \ No newline at end of file + assert cice_grid.ds.inputfile == mom_grid.path, "inputfile attribute incorrect ({cice_grid.ds.inputfile})" + assert ( + cice_grid.kmt_ds.inputfile == mom_grid.mask_path + ), "mask inputfile attribute incorrect ({cice_grid.kmt_ds.inputfile})" From 5b1bd80b61b941d3bd57db3dbeecaef8e08ae25e Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 14:54:39 +1100 Subject: [PATCH 05/25] correct black version --- om3utils/cice_grid.py | 1 - tests/test_cice_grid.py | 1 - 2 files changed, 2 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 617dddf..d1e2a0c 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -28,7 +28,6 @@ class cice_grid_nc: - """ Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc """ diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 96473ea..0f5c903 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -20,7 +20,6 @@ class MomGrid: - """Generate a sample tripole grid to use as test data""" def __init__(self, tmp_path): From 34fb292966df7d91452ba0f6ca64f0b96ea80d51 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:00:28 +1100 Subject: [PATCH 06/25] rm numpypi --- pyproject.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 722d48c..a12a452 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,7 +9,6 @@ dynamic = ["version"] dependencies = [ "f90nml", "ruamel.yaml", - "numpypi", "xarray", "numpy" ] From c39dc4647837157d67c9168d01179c975cfcdf72 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:15:00 +1100 Subject: [PATCH 07/25] submodules true in CI --- .github/workflows/ci.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 748358a..9f985c1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,6 +32,8 @@ jobs: steps: - uses: actions/checkout@v3 + with: + submodules: true - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v3 with: From 7ecf9eb6a70a2b58623fe68989e80f52143c6965 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:19:46 +1100 Subject: [PATCH 08/25] exclude esmgrids from flake8 --- .github/workflows/ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9f985c1..38677eb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -45,9 +45,9 @@ 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 --exclude esmgrids # 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 + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics --exclude esmgrids - name: Test with pytest run: | python -m pytest From e4da7c2a9694024921a35ed51d0a71182c86a904 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:21:33 +1100 Subject: [PATCH 09/25] add netcdf4 to dependencies --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index a12a452..eefc915 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,8 @@ dependencies = [ "f90nml", "ruamel.yaml", "xarray", - "numpy" + "numpy", + "netcdf4" ] [build-system] From b6f1d01dbef47de63b101c335064bda3e3d870ed Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:23:06 +1100 Subject: [PATCH 10/25] add pyproj to dependencies --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index eefc915..538aa37 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,8 @@ dependencies = [ "ruamel.yaml", "xarray", "numpy", - "netcdf4" + "netcdf4", + "pyproj" ] [build-system] From c2fd7c9d089b6f46290c5bd12fba34d15648a515 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:25:27 +1100 Subject: [PATCH 11/25] add shapely to dependencies --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 538aa37..e360771 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,8 @@ dependencies = [ "xarray", "numpy", "netcdf4", - "pyproj" + "pyproj", + "shapely" ] [build-system] From b2cf75016041ba674b4e93a486123a552c5d0b63 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:28:42 +1100 Subject: [PATCH 12/25] add scipy to dependencies --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index e360771..c2f95d8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,8 @@ dependencies = [ "numpy", "netcdf4", "pyproj", - "shapely" + "shapely", + "scipy" ] [build-system] From 8cc463a7f16a654f9f3c57e6a3e620ff4dc7a757 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 15:33:30 +1100 Subject: [PATCH 13/25] switch scipy to numpypi to dependencies --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c2f95d8..f405d0f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ dependencies = [ "netcdf4", "pyproj", "shapely", - "scipy" + "numpypi@git+https://github.com/underwoo/numpypi#egg=a79f1b53ab122340ea31ee487def8d6011224c31" ] [build-system] From 224697ae4adfa2513f7111e88b5d28a3b3849cb6 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 16:21:06 +1100 Subject: [PATCH 14/25] fix dependencies and tidy up --- .gitmodules | 6 ----- om3utils/cice_grid.py | 4 ---- om3utils/esmgrids | 1 - pyproject.toml | 3 ++- tests/ocean_model_grid_generator | 1 - tests/test_cice_grid.py | 41 +++++++++++++------------------- 6 files changed, 18 insertions(+), 38 deletions(-) delete mode 160000 om3utils/esmgrids delete mode 160000 tests/ocean_model_grid_generator diff --git a/.gitmodules b/.gitmodules index f595200..e69de29 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,6 +0,0 @@ -[submodule "om3utils/esmgrids"] - path = om3utils/esmgrids - url = https://github.com/COSIMA/esmgrids.git -[submodule "tests/ocean_model_grid_generator"] - path = tests/ocean_model_grid_generator - url = https://github.com/nikizadehgfdl/ocean_model_grid_generator.git diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index d1e2a0c..839b439 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -17,12 +17,8 @@ # File based on https://github.com/COSIMA/access-om2/blob/29118914d5224152ce286e0590394b231fea632e/tools/make_cice_grid.py import sys -import os import argparse -my_dir = os.path.dirname(os.path.realpath(__file__)) -sys.path.append(os.path.join(my_dir, "esmgrids")) - from esmgrids.mom_grid import MomGrid from esmgrids.cice_grid import CiceGrid diff --git a/om3utils/esmgrids b/om3utils/esmgrids deleted file mode 160000 index fb79a79..0000000 --- a/om3utils/esmgrids +++ /dev/null @@ -1 +0,0 @@ -Subproject commit fb79a794b6e0b2ee62e181d6cc336a381669fc50 diff --git a/pyproject.toml b/pyproject.toml index f405d0f..73bfce5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,8 @@ dependencies = [ "netcdf4", "pyproj", "shapely", - "numpypi@git+https://github.com/underwoo/numpypi#egg=a79f1b53ab122340ea31ee487def8d6011224c31" + "numpypi@git+https://github.com/underwoo/numpypi@a79f1b53ab122340ea31ee487def8d6011224c31", + "esmgrids@git+https://github.com/anton-seaice/esmgrids@0a563c2916674a6a963140c2cbf1c9c6ebfee1c7", ] [build-system] diff --git a/tests/ocean_model_grid_generator b/tests/ocean_model_grid_generator deleted file mode 160000 index 790069b..0000000 --- a/tests/ocean_model_grid_generator +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 790069b31f9791864ccd514a2b8f53f385a0452e diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 0f5c903..74332b2 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -2,13 +2,7 @@ import xarray as xr from numpy.testing import assert_allclose from numpy import deg2rad - -# im not sure how micael is importing om3utils? -import sys -import os - -my_dir = os.path.dirname(os.path.realpath(__file__)) -sys.path.append(os.path.join(my_dir, "..")) +from os import popen from om3utils.cice_grid import cice_grid_nc from ocean_model_grid_generator.ocean_grid_generator import main as ocean_grid_generator @@ -39,9 +33,10 @@ def __init__(self, tmp_path): warnings.simplefilter("ignore", category=RuntimeWarning) ocean_grid_generator(**args) # generates ocean_hgrid.nc + # open ocean_hgrid.nc self.ds = xr.open_dataset(self.path) - # an ocean mask with a random mask + # 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) @@ -59,6 +54,17 @@ def __init__(self, mom_grid, tmp_path): 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, tmp_path): # this generates the expected answers @@ -101,17 +107,6 @@ def test_grid_ds(mom_grid, tmp_path): return test_grid -# 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) - - # ---------------- # the tests in earnest: @@ -126,10 +121,6 @@ def test_cice_var_list(cice_grid, test_grid_ds): 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: - anom = cice_grid.ds[jVar].values - test_grid_ds[jVar].values - - print(f"{jVar} anom min: {anom.min()}, anom max: {anom.max()}") - assert_allclose(cice_grid.ds[jVar], test_grid_ds[jVar], rtol=1e-13, verbose=True, err_msg=f"{jVar} mismatch") @@ -138,7 +129,7 @@ def test_cice_kmt(mom_grid, cice_grid): mask = mom_grid.mask_ds.mask kmt = cice_grid.kmt_ds.kmt - assert_allclose(mask, kmt, rtol=1e-13, verbose=True, err_msg=f"mask mismatch") + assert_allclose(mask, kmt, rtol=1e-13, verbose=True, err_msg="mask mismatch") def test_cice_grid_attributes(cice_grid): @@ -197,7 +188,7 @@ def test_inputs_logged(cice_grid, mom_grid): assert hasattr(ds, "inputfile"), "inputfile attribute missing" assert hasattr(ds, "inputfile_md5"), "inputfile md5sum attribute missing" - sys_md5 = os.popen(f'md5sum {ds.inputfile} | cut -f 1 -d " "').read().strip() + sys_md5 = popen(f'md5sum {ds.inputfile} | cut -f 1 -d " "').read().strip() 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})" From cc4ebff2b23db7149eb22f239fc5fdeaad4693cc Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 16:50:49 +1100 Subject: [PATCH 15/25] fix dependencies --- .github/workflows/ci.yml | 3 --- pyproject.toml | 4 +--- tests/test_cice_grid.py | 14 +++++--------- 3 files changed, 6 insertions(+), 15 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 38677eb..15d46bb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,11 +29,8 @@ jobs: strategy: matrix: python-version: ['3.10', '3.11', '3.12'] - steps: - uses: actions/checkout@v3 - with: - submodules: true - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v3 with: diff --git a/pyproject.toml b/pyproject.toml index 73bfce5..aa0cdf6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,9 +12,7 @@ dependencies = [ "xarray", "numpy", "netcdf4", - "pyproj", - "shapely", - "numpypi@git+https://github.com/underwoo/numpypi@a79f1b53ab122340ea31ee487def8d6011224c31", + "ocean_model_grid_generator@git+https://github.com/nikizadehgfdl/ocean_model_grid_generator@790069b31f9791864ccd514a2b8f53f385a0452e", "esmgrids@git+https://github.com/anton-seaice/esmgrids@0a563c2916674a6a963140c2cbf1c9c6ebfee1c7", ] diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 74332b2..a51b3f5 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -2,12 +2,9 @@ import xarray as xr from numpy.testing import assert_allclose from numpy import deg2rad -from os import popen +from subprocess import run from om3utils.cice_grid import cice_grid_nc -from ocean_model_grid_generator.ocean_grid_generator import main as ocean_grid_generator - -import warnings # ---------------- # test data: @@ -29,10 +26,8 @@ def __init__(self, tmp_path): "gridfilename": self.path, } - with warnings.catch_warnings(): - warnings.simplefilter("ignore", category=RuntimeWarning) - ocean_grid_generator(**args) # generates ocean_hgrid.nc - + run(["ocean_grid_generator.py", "-r", "0.25", "--no_south_cap", "--ensure_nj_even", "-f", self.path]) + # open ocean_hgrid.nc self.ds = xr.open_dataset(self.path) @@ -188,7 +183,8 @@ def test_inputs_logged(cice_grid, mom_grid): assert hasattr(ds, "inputfile"), "inputfile attribute missing" assert hasattr(ds, "inputfile_md5"), "inputfile md5sum attribute missing" - sys_md5 = popen(f'md5sum {ds.inputfile} | cut -f 1 -d " "').read().strip() + 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})" From dd29aa5f4b7f664ecf535a2e9556f786dd07b91d Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 16:52:52 +1100 Subject: [PATCH 16/25] black --- tests/test_cice_grid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index a51b3f5..f80819b 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -27,7 +27,7 @@ def __init__(self, tmp_path): } run(["ocean_grid_generator.py", "-r", "0.25", "--no_south_cap", "--ensure_nj_even", "-f", self.path]) - + # open ocean_hgrid.nc self.ds = xr.open_dataset(self.path) @@ -183,7 +183,7 @@ def test_inputs_logged(cice_grid, mom_grid): 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 = 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}" From 9cb25e32be5c8c8cd18c80da50dbe88425eefb17 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 16:56:11 +1100 Subject: [PATCH 17/25] pyproj dep --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index aa0cdf6..ade66f7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ dependencies = [ "numpy", "netcdf4", "ocean_model_grid_generator@git+https://github.com/nikizadehgfdl/ocean_model_grid_generator@790069b31f9791864ccd514a2b8f53f385a0452e", - "esmgrids@git+https://github.com/anton-seaice/esmgrids@0a563c2916674a6a963140c2cbf1c9c6ebfee1c7", + "esmgrids@git+https://github.com/anton-seaice/esmgrids@53e336d90cf6aad386721df574e8f9d33f820806", ] [build-system] From 97a4a87a203a5cd940eac6c6e8119349f95eab12 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Wed, 3 Apr 2024 16:58:21 +1100 Subject: [PATCH 18/25] shapely dep --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ade66f7..72acd24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ dependencies = [ "numpy", "netcdf4", "ocean_model_grid_generator@git+https://github.com/nikizadehgfdl/ocean_model_grid_generator@790069b31f9791864ccd514a2b8f53f385a0452e", - "esmgrids@git+https://github.com/anton-seaice/esmgrids@53e336d90cf6aad386721df574e8f9d33f820806", + "esmgrids@git+https://github.com/anton-seaice/esmgrids@e5791f17a699627566687021e250aac32567ecd4", ] [build-system] From 47edb5d9fa97038a11b8d8afc8f1075f7731bacd Mon Sep 17 00:00:00 2001 From: anton-climate Date: Thu, 4 Apr 2024 09:36:49 +1100 Subject: [PATCH 19/25] tidyup --- om3utils/cice_grid.py | 9 ++++----- tests/test_cice_grid.py | 28 ++++++++++++---------------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 839b439..792b8f8 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -23,7 +23,7 @@ from esmgrids.cice_grid import CiceGrid -class cice_grid_nc: +class CiceGridNc: """ Create CICE grid.nc and kmt.nc from MOM ocean_hgrid.nc and ocean_mask.nc """ @@ -68,7 +68,7 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): if __name__ == "__main__": - from utils import md5sum + from utils import md5sum # load from file parser = argparse.ArgumentParser() parser.add_argument("ocean_hgrid", help="ocean_hgrid.nc file") @@ -77,10 +77,9 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): args = vars(parser.parse_args()) - grid = cice_grid_nc() + grid = CiceGridNc() sys.exit(grid.build_from_mom(**args)) else: - # for testing - from .utils import md5sum + from .utils import md5sum #load from package diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index f80819b..0b6069a 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -4,7 +4,7 @@ from numpy import deg2rad from subprocess import run -from om3utils.cice_grid import cice_grid_nc +from om3utils.cice_grid import CiceGridNc # ---------------- # test data: @@ -18,15 +18,13 @@ def __init__(self, tmp_path): self.mask_path = str(tmp_path) + "/ocean_mask.nc" # generate an tripolar grid as test data - args = { - "inverse_resolution": 0.25, # 4 degree grid - "no_south_cap": True, - "ensure_nj_even": True, - "match_dy": ["bp", "so", "p125sc", ""], - "gridfilename": self.path, - } - - run(["ocean_grid_generator.py", "-r", "0.25", "--no_south_cap", "--ensure_nj_even", "-f", self.path]) + 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) @@ -43,7 +41,7 @@ class CiceGrid: def __init__(self, mom_grid, tmp_path): self.path = str(tmp_path) + "/grid.nc" self.kmt_path = str(tmp_path) + "/kmt.nc" - cice_grid = cice_grid_nc(self.path, self.kmt_path) + 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) @@ -61,7 +59,7 @@ def cice_grid(mom_grid, tmp_path): @pytest.fixture -def test_grid_ds(mom_grid, tmp_path): +def test_grid_ds(mom_grid): # this generates the expected answers ds = mom_grid.ds @@ -187,7 +185,5 @@ def test_inputs_logged(cice_grid, mom_grid): 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})" - assert ( - cice_grid.kmt_ds.inputfile == mom_grid.mask_path - ), "mask inputfile attribute incorrect ({cice_grid.kmt_ds.inputfile})" + 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})" From 94903cf6d60256d2f5e92399a0b66ef029e7e83f Mon Sep 17 00:00:00 2001 From: anton-climate Date: Thu, 4 Apr 2024 10:09:58 +1100 Subject: [PATCH 20/25] more tidyup --- .github/workflows/ci.yml | 4 ++-- .gitmodules | 0 om3utils/cice_grid.py | 5 +++-- pyproject.toml | 1 - 4 files changed, 5 insertions(+), 5 deletions(-) delete mode 100644 .gitmodules diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 15d46bb..9496b0d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,9 +42,9 @@ 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 --exclude esmgrids + 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 --exclude esmgrids + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - name: Test with pytest run: | python -m pytest diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index e69de29..0000000 diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 792b8f8..29c20bf 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -8,8 +8,9 @@ - ocean_hgrid: Path to the MOM super grid NetCDF file. - ocean_mask: Path to the corresponding mask NetCDF file. -Dependencies: Developed using -'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10' +Dependencies: +- 'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10'. +- or 'pyproject.toml' file """ diff --git a/pyproject.toml b/pyproject.toml index 72acd24..c7f4b43 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,7 +40,6 @@ test = [ [tool.pytest.ini_options] addopts = ["--cov=om3utils", "--cov-report=term", "--cov-report=xml"] testpaths = ["tests"] -norecursedirs = ["tests/ocean_model_grid_generator"] [tool.coverage.run] omit = ["om3utils/__init__.py", "om3utils/_version.py"] From 8be55091cb175fccb8c7bbef151d630170b3f488 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Thu, 4 Apr 2024 10:11:05 +1100 Subject: [PATCH 21/25] more tidyup --- om3utils/cice_grid.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 29c20bf..929e258 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -69,12 +69,13 @@ def build_from_mom(self, ocean_hgrid, ocean_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? + # to-do: add argument for CRS & output filenames? args = vars(parser.parse_args()) From 6b8c64498a3213e9b31603d989a5fcfb5f843132 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Thu, 4 Apr 2024 12:07:01 +1100 Subject: [PATCH 22/25] black --- om3utils/cice_grid.py | 6 +++--- tests/test_cice_grid.py | 26 +++++++++++++++++--------- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 929e258..409045d 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -69,8 +69,8 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): if __name__ == "__main__": - #command line arguments - from utils import md5sum # load from file + # command line arguments + from utils import md5sum # load from file parser = argparse.ArgumentParser() parser.add_argument("ocean_hgrid", help="ocean_hgrid.nc file") @@ -84,4 +84,4 @@ def build_from_mom(self, ocean_hgrid, ocean_mask): sys.exit(grid.build_from_mom(**args)) else: - from .utils import md5sum #load from package + from .utils import md5sum # load from package diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 0b6069a..ba6da2d 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -18,13 +18,17 @@ def __init__(self, tmp_path): 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 - ]) + 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) @@ -185,5 +189,9 @@ def test_inputs_logged(cice_grid, mom_grid): 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})" + 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})" From 29c22918b2dc77207eb8d8d25bb086e4ed679d39 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Fri, 5 Apr 2024 10:52:25 +1100 Subject: [PATCH 23/25] Create scripts folder for reference scripts --- README.md | 1 + om3utils/pbs_make_cice_grids.sh | 31 ------------------------------- scripts/README.md | 5 +++++ scripts/pbs_make_cice_grids.sh | 33 +++++++++++++++++++++++++++++++++ 4 files changed, 39 insertions(+), 31 deletions(-) delete mode 100644 om3utils/pbs_make_cice_grids.sh create mode 100644 scripts/README.md create mode 100644 scripts/pbs_make_cice_grids.sh 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/pbs_make_cice_grids.sh b/om3utils/pbs_make_cice_grids.sh deleted file mode 100644 index 6655a38..0000000 --- a/om3utils/pbs_make_cice_grids.sh +++ /dev/null @@ -1,31 +0,0 @@ -#!/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 \ No newline at end of file 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 From d836e70352316bfeacce9e6bcb3a58db3c139f8d Mon Sep 17 00:00:00 2001 From: Anton Steketee <79179784+anton-seaice@users.noreply.github.com> Date: Thu, 11 Apr 2024 10:16:21 +1000 Subject: [PATCH 24/25] Apply suggestions from code review Co-authored-by: Andrew Kiss <31054815+aekiss@users.noreply.github.com> --- tests/test_cice_grid.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index ba6da2d..012a161 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -17,7 +17,7 @@ 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 + # generate a tripolar grid as test data run( [ "ocean_grid_generator.py", @@ -84,15 +84,14 @@ def test_grid_ds(mom_grid): # 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 + # angle at 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() + 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 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") From 76665f76f9c1ffec6e639a6a9c1e1f1218f550f1 Mon Sep 17 00:00:00 2001 From: anton-climate Date: Thu, 11 Apr 2024 11:18:41 +1000 Subject: [PATCH 25/25] Bit of a refactor to hopefully improve clarity --- om3utils/cice_grid.py | 2 +- tests/test_cice_grid.py | 43 ++++++++++++++++++++++------------------- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/om3utils/cice_grid.py b/om3utils/cice_grid.py index 409045d..91d8ed3 100644 --- a/om3utils/cice_grid.py +++ b/om3utils/cice_grid.py @@ -9,7 +9,7 @@ - ocean_mask: Path to the corresponding mask NetCDF file. Dependencies: -- 'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10'. +- 'module use /g/data/hh5/public/modules ; module load conda/analysis3-23.10' + pip install esmgrids - or 'pyproject.toml' file """ diff --git a/tests/test_cice_grid.py b/tests/test_cice_grid.py index 012a161..43acb80 100644 --- a/tests/test_cice_grid.py +++ b/tests/test_cice_grid.py @@ -65,40 +65,43 @@ def cice_grid(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() - # 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))) + 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 edge of cells + # 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 edge of cells + # length of right (eastern) 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 at 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"] =ds.area.coarsen(ny=2).sum().coarsen(nx=2).sum() + 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 wrap in latitude and fold on longitude to calculate this - 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") + # 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") - area_folded = xr.concat([area_wrapped.isel(ny=slice(1, None)), top_row], 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_folded.coarsen(ny=2).sum().coarsen(nx=2).sum() + test_grid["uarea"] = area_wrapped.coarsen(ny=2).sum().coarsen(nx=2).sum() return test_grid