Add mom_to_cice metadata and tests

README.md

@@ -1,6 +1,7 @@
 
 [![Code Health](https://landscape.io/github/DoublePrecision/esmgrids/master/landscape.svg?style=flat)](https://landscape.io/github/DoublePrecision/esmgrids/master)
-[![Build Status](https://travis-ci.org/DoublePrecision/esmgrids.svg?branch=master)](https://travis-ci.org/DoublePrecision/esmgrids)
+
+![Code Health](https://github.com/COSIMA/esmgrids/actions/workflows/ci.yml/badge.svg)
 
 # esmgrids
 
@@ -19,9 +20,8 @@ Grids currently supported:
 ## To run the tests
 
 ```
-conda env create -n grids python3
-source activate grids
-python -m pytest
+pip install '.[tests]'
+python -m pytest -m "not broken"
 ```
 
 Warning: this will download a rather large tarball of test inputs.

esmgrids/_version.py

@@ -0,0 +1,16 @@
+# file generated by setuptools_scm
+# don't change, don't track in version control
+TYPE_CHECKING = False
+if TYPE_CHECKING:
+    from typing import Tuple, Union
+    VERSION_TUPLE = Tuple[Union[int, str], ...]
+else:
+    VERSION_TUPLE = object
+
+version: str
+__version__: str
+__version_tuple__: VERSION_TUPLE
+version_tuple: VERSION_TUPLE
+
+__version__ = version = '0.1.dev90+gb53e07d.d20240412'
+__version_tuple__ = version_tuple = (0, 1, 'dev90', 'gb53e07d.d20240412')

esmgrids/base_grid.py

@@ -1,7 +1,7 @@
 import numpy as np
 import netCDF4 as nc
 
-from .util import calc_area_of_polygons
+from esmgrids.util import calc_area_of_polygons
 
 
 class BaseGrid(object):
@@ -49,7 +49,7 @@ def __init__(self, **kwargs):
         self.description = kwargs.get("description", "")
         self.calc_areas = kwargs.get("calc_areas", True)
 
-        if len(x_t.shape) == 1:
+        if x_t is not None and len(x_t.shape) == 1:
             # Tile x_t
             self.x_t = np.tile(x_t, (y_t.shape[0], 1))
             self.y_t = np.tile(y_t, (x_t.shape[0], 1))

esmgrids/cice_grid.py

@@ -1,31 +1,18 @@
 import numpy as np
 import netCDF4 as nc
 
-from .base_grid import BaseGrid
+from esmgrids.base_grid import BaseGrid
+from esmgrids.mom_grid import MomGrid
+from esmgrids.util import *
 
-from os import environ
-from datetime import datetime
-
-from .util import *
-
-
-def create_nc(filename):
-
-    f = nc.Dataset(filename, 'w')
-
-    f.timeGenerated = f"{datetime.now()}"
-    f.created_by = f"{environ.get('USER')}"
-    if is_git_repo():
-        git_url, _ , git_hash = git_info()
-        f.history = f"Created using commit {git_hash} of {git_url}"
-
-    return f 
 
 class CiceGrid(BaseGrid):
 
     def __init__(self, **kwargs):
         self.type = "Arakawa B"
         self.full_name = "CICE"
+        # self.grid_file = grid_file
+        # self.mask_file = mask_file
 
         super(CiceGrid, self).__init__(**kwargs)
 
@@ -85,30 +72,32 @@ def fromfile(cls, h_grid_def, mask_file=None, description="CICE tripolar"):
     def create_gridnc(self, grid_filename):
         self.grid_f = create_nc(grid_filename)
         return True
-    
+
     def create_masknc(self, mask_filename):
         self.mask_f = create_nc(mask_filename)
         return True
-    
+
     def create_2d_grid_var(self, name):
-        #set chunksizes based on OM2 config
-        #To-do: load these from a configuration file?
-        if (self.num_lon_points==360): #1deg
-            chunksizes=(150,180)
-        elif (self.num_lon_points==1440): #0.25deg
-            chunksizes=(540,720)
-        elif (self.num_lon_points==3600): #0.01deg
-            chunksizes=(270,360)
-        else :
-            chunksizes=None
+        # set chunksizes based on OM2 config
+        # To-do: load these from a configuration file?
+        if self.num_lon_points == 360:  # 1deg
+            chunksizes = (150, 180)
+        elif self.num_lon_points == 1440:  # 0.25deg
+            chunksizes = (540, 720)
+        elif self.num_lon_points == 3600:  # 0.01deg
+            chunksizes = (270, 360)
+        else:
+            chunksizes = None
 
         return self.grid_f.createVariable(
-            name, 'f8', dimensions=('ny', 'nx'),
-            compression='zlib', complevel=1,
-            chunksizes=chunksizes
+            name,
+            "f8",
+            dimensions=("ny", "nx"),
+            compression="zlib",
+            complevel=1,
+            chunksizes=chunksizes,
         )
 
-
     def write(self):
         """
         Write out CICE grid to netcdf.
@@ -117,69 +106,68 @@ def write(self):
         f = self.grid_f
 
         # Create dimensions.
-        f.createDimension('nx', self.num_lon_points)
-        #nx is the grid_longitude but doesn't have a value other than its index
-        f.createDimension('ny', self.num_lat_points)
-        #ny is the grid_latitude but doesn't have a value other than its index
+        f.createDimension("nx", self.num_lon_points)
+        # nx is the grid_longitude but doesn't have a value other than its index
+        f.createDimension("ny", self.num_lat_points)
+        # ny is the grid_latitude but doesn't have a value other than its index
 
         # Make all CICE grid variables.
         # names are based on https://cfconventions.org/Data/cf-standard-names/current/build/cf-standard-name-table.html
         f.Conventions = "CF-1.6"
-        
-        ulat = self.create_2d_grid_var('ulat')
+
+        ulat = self.create_2d_grid_var("ulat")
         ulat.units = "radians"
         ulat.long_name = "Latitude of U points"
         ulat.standard_name = "latitude"
-        ulon = self.create_2d_grid_var('ulon')
+        ulon = self.create_2d_grid_var("ulon")
         ulon.units = "radians"
         ulon.long_name = "Longitude of U points"
         ulon.standard_name = "longitude"
-        tlat = self.create_2d_grid_var('tlat')
+        tlat = self.create_2d_grid_var("tlat")
         tlat.units = "radians"
         tlat.long_name = "Latitude of T points"
         tlat.standard_name = "latitude"
-        tlon = self.create_2d_grid_var('tlon')
+        tlon = self.create_2d_grid_var("tlon")
         tlon.units = "radians"
         tlon.long_name = "Longitude of T points"
         tlon.standard_name = "longitude"
-        
-        htn = self.create_2d_grid_var('htn')
+
+        htn = self.create_2d_grid_var("htn")
         htn.units = "cm"
         htn.long_name = "Width of T cells on North side."
         htn.coordinates = "ulat tlon"
         htn.grid_mapping = "crs"
-        hte = self.create_2d_grid_var('hte')
+        hte = self.create_2d_grid_var("hte")
         hte.units = "cm"
         hte.long_name = "Width of T cells on East side."
         hte.coordinates = "tlat ulon"
         hte.grid_mapping = "crs"
-        
-        angle = self.create_2d_grid_var('angle')
+
+        angle = self.create_2d_grid_var("angle")
         angle.units = "radians"
         angle.long_name = "Rotation angle of U cells."
         angle.standard_name = "angle_of_rotation_from_east_to_x"
         angle.coordinates = "ulat ulon"
         angle.grid_mapping = "crs"
-        angleT = self.create_2d_grid_var('angleT')
+        angleT = self.create_2d_grid_var("angleT")
         angleT.units = "radians"
         angleT.long_name = "Rotation angle of T cells."
         angleT.standard_name = "angle_of_rotation_from_east_to_x"
         angleT.coordinates = "tlat tlon"
         angleT.grid_mapping = "crs"
-        
-        area_t = self.create_2d_grid_var('tarea')
+
+        area_t = self.create_2d_grid_var("tarea")
         area_t.units = "m^2"
         area_t.long_name = "Area of T cells."
         area_t.standard_name = "cell_area"
         area_t.coordinates = "tlat tlon"
         area_t.grid_mapping = "crs"
-        area_u = self.create_2d_grid_var('uarea')
+        area_u = self.create_2d_grid_var("uarea")
         area_u.units = "m^2"
         area_u.long_name = "Area of U cells."
         area_u.standard_name = "cell_area"
         area_u.coordinates = "ulat ulon"
         area_u.grid_mapping = "crs"
-
 
         area_t[:] = self.area_t[:]
         area_u[:] = self.area_u[:]
@@ -200,22 +188,78 @@ def write(self):
         f.close()
 
     def write_mask(self):
-
         """
         Write out CICE mask/kmt to netcdf.
         """
 
         f = self.mask_f
-        f.createDimension('nx', self.num_lon_points)
-        f.createDimension('ny', self.num_lat_points)
-        mask = f.createVariable(
-            'kmt', 'f8', dimensions=('ny', 'nx'),
-            compression='zlib', complevel=1
-        )
-
+        f.createDimension("nx", self.num_lon_points)
+        f.createDimension("ny", self.num_lat_points)
+        mask = f.createVariable("kmt", "f8", dimensions=("ny", "nx"), compression="zlib", complevel=1)
+
         mask.grid_mapping = "crs"
         mask.standard_name = "sea_binary_mask"
-        
+
         # CICE uses 0 as masked, whereas internally we use 1 as masked.
-        mask[:] = (1 - self.mask_t)
+        mask[:] = 1 - self.mask_t
         f.close()
+
+def cice_from_mom(ocean_hgrid, ocean_mask, grid_file="grid.nc", mask_file="kmt.nc"):
+
+    mom = MomGrid.fromfile(ocean_hgrid, mask_file=ocean_mask)
+
+    cice = CiceGrid.fromgrid(mom)
+
+    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}"
+
+    # 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(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}"
+
+    # 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__":
+    import argparse
+    import sys
+    # command line arguments
+    # import os, sys
+    # sys.path.append(os.path.dirname(__file__))
+
+    # print(os.path.dirname(__file__))
+
+    # from base_grid import BaseGrid
+    # from mom_grid import MomGrid
+    # from util import *  # 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())
+
+    sys.exit(cice_from_mom(**args))
+
+# else:
+#     # for testing

esmgrids/mom_grid.py

@@ -1,7 +1,7 @@
 import numpy as np
 import netCDF4 as nc
 
-from .base_grid import BaseGrid
+from esmgrids.base_grid import BaseGrid
 
 
 class MomGrid(BaseGrid):
@@ -85,8 +85,8 @@ def fromfile(cls, h_grid_def, v_grid_def=None, mask_file=None, calc_areas=True,
                 dx_un = dx_ext[3::2, 1::2] + dx_ext[3::2, 2::2]
                 dy_ue = dy_ext[1::2, 3::2] + dy_ext[2::2, 3::2]
 
-                angle_dx = f.variables['angle_dx'][:]
-                
+                angle_dx = f.variables["angle_dx"][:]
+
                 angle_t = angle_dx[1::2, 1::2]
                 angle_u = angle_dx[2::2, 2::2]