Cice grid generation

.github/workflows/ci.yml

@@ -1,7 +1,7 @@
 # This workflow will install Python dependencies, run tests and lint with a single version of Python
 # For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python
 
-name: OM3Utils
+name: esmgrids
 
 on:
   push:

README.md

@@ -1,6 +1,4 @@
-
-[![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 +17,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/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))
@@ -58,7 +58,7 @@ def __init__(self, **kwargs):
             self.x_t = x_t
             self.y_t = y_t
 
-        if dx_t is not None and len(dx_t.shape) == 1:
+        if len(dx_t.shape) == 1:
             # Tile dx_t
             self.dx_t = np.tile(x_t, (dy_t.shape[0], 1))
             self.dy_t = np.tile(y_t, (dx_t.shape[0], 1))

esmgrids/cice_grid.py

@@ -1,7 +1,9 @@
 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 *
 
 
 class CiceGrid(BaseGrid):
@@ -65,66 +67,105 @@ def fromfile(cls, h_grid_def, mask_file=None, description="CICE tripolar"):
             description=description,
         )
 
-    def write(self, grid_filename, mask_filename):
+    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
+
+        return self.grid_f.createVariable(
+            name,
+            "f8",
+            dimensions=("ny", "nx"),
+            compression="zlib",
+            complevel=1,
+            chunksizes=chunksizes,
+        )
+
+    def write(self):
         """
         Write out CICE grid to netcdf.
         """
 
-        f = nc.Dataset(grid_filename, "w")
+        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)
-        f.createDimension("nc", 4)
+        # ny is the grid_latitude but doesn't have a value other than its index
 
         # Make all CICE grid variables.
-        ulat = f.createVariable("ulat", "f8", dimensions=("ny", "nx"))
+        # 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.units = "radians"
-        ulat.title = "Latitude of U points"
-        ulon = f.createVariable("ulon", "f8", dimensions=("ny", "nx"))
+        ulat.long_name = "Latitude of U points"
+        ulat.standard_name = "latitude"
+        ulon = self.create_2d_grid_var("ulon")
         ulon.units = "radians"
-        ulon.title = "Longitude of U points"
-        tlat = f.createVariable("tlat", "f8", dimensions=("ny", "nx"))
+        ulon.long_name = "Longitude of U points"
+        ulon.standard_name = "longitude"
+        tlat = self.create_2d_grid_var("tlat")
         tlat.units = "radians"
-        tlat.title = "Latitude of T points"
-        tlon = f.createVariable("tlon", "f8", dimensions=("ny", "nx"))
+        tlat.long_name = "Latitude of T points"
+        tlat.standard_name = "latitude"
+        tlon = self.create_2d_grid_var("tlon")
         tlon.units = "radians"
-        tlon.title = "Longitude of T points"
-
-        if self.clon_t is not None:
-            clon_t = f.createVariable("clon_t", "f8", dimensions=("nc", "ny", "nx"))
-            clon_t.units = "radians"
-            clon_t.title = "Longitude of T cell corners"
-            clat_t = f.createVariable("clat_t", "f8", dimensions=("nc", "ny", "nx"))
-            clat_t.units = "radians"
-            clat_t.title = "Latitude of T cell corners"
-            clon_u = f.createVariable("clon_u", "f8", dimensions=("nc", "ny", "nx"))
-            clon_u.units = "radians"
-            clon_u.title = "Longitude of U cell corners"
-            clat_u = f.createVariable("clat_u", "f8", dimensions=("nc", "ny", "nx"))
-            clat_u.units = "radians"
-            clat_u.title = "Latitude of U cell corners"
-
-        htn = f.createVariable("htn", "f8", dimensions=("ny", "nx"))
+        tlon.long_name = "Longitude of T points"
+        tlon.standard_name = "longitude"
+
+        htn = self.create_2d_grid_var("htn")
         htn.units = "cm"
-        htn.title = "Width of T cells on North side."
-        hte = f.createVariable("hte", "f8", dimensions=("ny", "nx"))
+        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.units = "cm"
-        hte.title = "Width of T cells on East side."
+        hte.long_name = "Width of T cells on East side."
+        hte.coordinates = "tlat ulon"
+        hte.grid_mapping = "crs"
 
-        angle = f.createVariable("angle", "f8", dimensions=("ny", "nx"))
+        angle = self.create_2d_grid_var("angle")
         angle.units = "radians"
-        angle.title = "Rotation angle of U cells."
-        angleT = f.createVariable("angleT", "f8", dimensions=("ny", "nx"))
+        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.units = "radians"
-        angleT.title = "Rotation angle of T cells."
+        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 = f.createVariable("tarea", "f8", dimensions=("ny", "nx"))
+        area_t = self.create_2d_grid_var("tarea")
         area_t.units = "m^2"
-        area_t.title = "Area of T cells."
-        area_u = f.createVariable("uarea", "f8", dimensions=("ny", "nx"))
+        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.units = "m^2"
-        area_u.title = "Area of U cells."
+        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[:]
@@ -135,12 +176,6 @@ def write(self, grid_filename, mask_filename):
         ulat[:] = np.deg2rad(self.y_u)
         ulon[:] = np.deg2rad(self.x_u)
 
-        if self.clon_t is not None:
-            clon_t[:] = np.deg2rad(self.clon_t)
-            clat_t[:] = np.deg2rad(self.clat_t)
-            clon_u[:] = np.deg2rad(self.clon_u)
-            clat_u[:] = np.deg2rad(self.clat_u)
-
         # Convert from m to cm.
         htn[:] = self.dx_tn[:] * 100.0
         hte[:] = self.dy_te[:] * 100.0
@@ -150,9 +185,68 @@ def write(self, grid_filename, mask_filename):
 
         f.close()
 
-        with nc.Dataset(mask_filename, "w") as f:
-            f.createDimension("nx", self.num_lon_points)
-            f.createDimension("ny", self.num_lat_points)
-            mask = f.createVariable("kmt", "f8", dimensions=("ny", "nx"))
-            # CICE uses 0 as masked, whereas internally we use 1 as masked.
-            mask[:] = 1 - self.mask_t
+    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)
+
+        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
+        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
+
+    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))

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):
@@ -86,9 +86,8 @@ def fromfile(cls, h_grid_def, v_grid_def=None, mask_file=None, calc_areas=True,
                 dy_ue = dy_ext[1::2, 3::2] + dy_ext[2::2, 3::2]
 
                 angle_dx = f.variables["angle_dx"][:]
-                # The angle of rotation is a corner cell centres and applies to
-                # both t and u cells.
-                angle_t = angle_dx[2::2, 2::2]
+
+                angle_t = angle_dx[1::2, 1::2]
                 angle_u = angle_dx[2::2, 2::2]
 
                 area = f.variables["area"][:]

esmgrids/util.py

@@ -1,6 +1,11 @@
 import numpy as np
 import pyproj
 from shapely.geometry import shape
+import subprocess
+import os
+from datetime import datetime
+from netCDF4 import Dataset
+
 
 proj_str = "+proj=laea +lat_0={} +lon_0={} +ellps=sphere"
 
@@ -39,3 +44,45 @@ def calc_area_of_polygons(clons, clats):
     assert np.min(areas) > 0
 
     return areas
+
+
+def is_git_repo():
+    """
+    Return True/False depending on whether or not the current directory is a git repo.
+    """
+
+    return subprocess.call(["git", "-C", ".", "status"], stderr=subprocess.STDOUT, stdout=open(os.devnull, "w")) == 0
+
+
+def git_info():
+    """
+    Return the git repo origin url, relative path to this file, and latest commit hash.
+    """
+
+    url = subprocess.check_output(["git", "remote", "get-url", "origin"]).decode("ascii").strip()
+    top_level_dir = subprocess.check_output(["git", "rev-parse", "--show-toplevel"]).decode("ascii").strip()
+    rel_path = os.path.relpath(__file__, top_level_dir)
+    hash = subprocess.check_output(["git", "rev-parse", "HEAD"]).decode("ascii").strip()
+
+    return url, rel_path, hash
+
+
+def create_nc(filename):
+
+    f = Dataset(filename, "w")
+
+    f.timeGenerated = f"{datetime.now()}"
+    f.created_by = f"{os.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
+
+
+def md5sum(filename):
+    from hashlib import md5
+    from mmap import mmap, ACCESS_READ
+
+    with open(filename) as file, mmap(file.fileno(), 0, access=ACCESS_READ) as file:
+        return md5(file).hexdigest()

pyproject.toml

@@ -23,6 +23,8 @@ dependencies = [
     "netcdf4",
     "shapely",
     "pyproj",
+    "xarray",
+    "ocean_model_grid_generator@git+https://github.com/nikizadehgfdl/ocean_model_grid_generator@790069b31f9791864ccd514a2b8f53f385a0452e"
 ]
 
 [build-system]