Merge pull request #187 from xylar/honor-cf-bnds

Use CF bounds for grid-cell corners where available

Author: xylar

ci/recipe/meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "pyremap" %}
-{% set version = "2.3.0" %}
+{% set version = "2.4.0" %}
 {% set python_min = "3.10" %}
 
 package:

docs/developer_guide/api.md

@@ -13,6 +13,8 @@ This page provides an auto-generated summary of the pyremap API.
    :toctree: generated/
 
 
+   get_corners_1d
+   get_corners_2d
    interp_extrap_corner
    interp_extrap_corners_2d
 

docs/mesh_descriptors/lat_lon_2d_grid_descriptor.md

@@ -15,6 +15,24 @@ The `LatLon2DGridDescriptor` class is used for grids where latitude and longitud
 - `read`: Reads a 2D latitude-longitude grid from a file.
 - `to_scrip`: Converts the grid to a SCRIP file.
 
+## Grid-Cell Corners
+As with {py:class}`LatLonGridDescriptor <pyremap.LatLonGridDescriptor>`,
+`read()` uses the CF `bounds` of the latitude and longitude variables to find
+grid-cell corners when they are available.  For 2D coordinates, the bounds
+give the 4 vertices of each cell:
+```
+double lat(y, x) ;
+    lat:units = "degrees_north" ;
+    lat:bounds = "lat_bnds" ;
+double lat_bnds(y, x, nv) ;
+```
+CF does not say which vertex comes first or which direction the 4 vertices are
+traversed in, so pyremap works this out from the bounds themselves.  Both
+latitude and longitude must have bounds, and neighboring cells must share
+vertices, since the grid is described by 2D arrays of corners.  Otherwise,
+corners are interpolated and extrapolated from cell centers and a warning is
+raised.
+
 ## Example
 ```python
 from pyremap import LatLon2DGridDescriptor

docs/mesh_descriptors/lat_lon_grid_descriptor.md

@@ -16,6 +16,24 @@ The `LatLonGridDescriptor` class is used to describe a regular latitude-longitud
 - `create`: Creates a latitude-longitude grid programmatically.
 - `to_scrip`: Converts the grid to a SCRIP file.
 
+## Grid-Cell Corners
+Remapping needs the corners of each grid cell, not just the cell centers.
+When `read()` is used, corners come from the CF `bounds` attribute of the
+latitude and longitude variables if it is present:
+```
+double lat(lat) ;
+    lat:units = "degrees_north" ;
+    lat:bounds = "lat_bnds" ;
+double lat_bnds(lat, nbnd) ;
+```
+The bounds must describe contiguous cells (the upper edge of each cell is the
+lower edge of the next), since a 1D lat/lon grid is described by 1D arrays of
+corners.  If the `bounds` attribute is missing, points to a variable that is
+not in the file, has the wrong shape, or describes cells with gaps or overlaps
+between them, corners are instead interpolated between cell centers and
+extrapolated at the ends of the grid, and a warning is raised in all but the
+first of these cases.
+
 ## Example
 ```python
 from pyremap import LatLonGridDescriptor

docs/mesh_descriptors/projection_grid_descriptor.md

@@ -15,6 +15,13 @@ The `ProjectionGridDescriptor` class describes grids defined by map projections.
 - `create`: Creates a projection grid programmatically.
 - `to_scrip`: Converts the grid to a SCRIP file.
 
+## Grid-Cell Corners
+`read()` uses the CF `bounds` of the `x` and `y` variables to find the corners
+of each grid cell in projection space when they are available and describe
+contiguous cells.  Otherwise, corners are interpolated between cell centers
+and extrapolated at the ends of the grid.  Corners are transformed from
+projection space to latitude and longitude by `to_scrip()`.
+
 ## Example
 ```python
 from pyremap import ProjectionGridDescriptor

pyremap/descriptor/__init__.py

@@ -36,6 +36,12 @@
 from pyremap.descriptor.projection_grid_descriptor import (
     ProjectionGridDescriptor as ProjectionGridDescriptor,
 )
+from pyremap.descriptor.utility import (
+    get_corners_1d as get_corners_1d,
+)
+from pyremap.descriptor.utility import (
+    get_corners_2d as get_corners_2d,
+)
 from pyremap.descriptor.utility import (
     interp_extrap_corner as interp_extrap_corner,
 )

pyremap/descriptor/lat_lon_2d_grid_descriptor.py

@@ -18,7 +18,7 @@
 from pyremap.descriptor.utility import (
     add_history,
     expand_scrip,
-    interp_extrap_corners_2d,
+    get_corners_2d,
     round_res,
     unwrap_corners,
 )
@@ -87,6 +87,11 @@ def read(
         """
         Read the lat-lon grid from a file with the given lat/lon var names.
 
+        Grid-cell corners come from the CF ``bounds`` of the latitude and
+        longitude variables if they are available and neighboring cells share
+        vertices.  Otherwise, corners are interpolated and extrapolated from
+        the cell centers.
+
         Parameters
         ----------
         filename : str, optional
@@ -127,9 +132,10 @@ def read(
         else:
             descriptor.units = 'radians'
 
-        # interp/extrap corners
-        descriptor.lon_corner = interp_extrap_corners_2d(descriptor.lon)
-        descriptor.lat_corner = interp_extrap_corners_2d(descriptor.lat)
+        # use CF bounds if available, otherwise interp/extrap corners
+        descriptor.lat_corner, descriptor.lon_corner = get_corners_2d(
+            ds, lat_var_name, lon_var_name
+        )
 
         descriptor._set_coords(
             lat_var_name,

pyremap/descriptor/lat_lon_grid_descriptor.py

@@ -18,7 +18,7 @@
 from pyremap.descriptor.utility import (
     add_history,
     expand_scrip,
-    interp_extrap_corner,
+    get_corners_1d,
     round_res,
     unwrap_corners,
 )
@@ -121,6 +121,11 @@ def read(
         """
         Read the lat-lon grid from a file with the given lat/lon var names.
 
+        Grid-cell corners come from the CF ``bounds`` of the latitude and
+        longitude variables if they are available and describe contiguous
+        cells.  Otherwise, corners are interpolated and extrapolated from the
+        cell centers.
+
         Parameters
         ----------
         filename : str, optional
@@ -157,9 +162,9 @@ def read(
         else:
             descriptor.units = 'radians'
 
-        # interp/extrap corners
-        descriptor.lon_corner = interp_extrap_corner(descriptor.lon)
-        descriptor.lat_corner = interp_extrap_corner(descriptor.lat)
+        # use CF bounds if available, otherwise interp/extrap corners
+        descriptor.lon_corner = get_corners_1d(ds, lon_var_name)
+        descriptor.lat_corner = get_corners_1d(ds, lat_var_name)
 
         descriptor._set_coords(
             lat_var_name,

pyremap/descriptor/projection_grid_descriptor.py

@@ -19,6 +19,7 @@
 from pyremap.descriptor.utility import (
     add_history,
     expand_scrip,
+    get_corners_1d,
     interp_extrap_corner,
     unwrap_corners,
 )
@@ -98,7 +99,11 @@ def read(
         """
         Given a grid file with x and y coordinates defining the axes of the
         logically rectangular grid, read in the x and y coordinates and
-        interpolate/extrapolate to locate corners.
+        locate the corners.
+
+        Corners come from the CF ``bounds`` of the x and y variables if they
+        are available and describe contiguous cells.  Otherwise, corners are
+        interpolated and extrapolated from the cell centers.
 
         Parameters
         ----------
@@ -135,9 +140,9 @@ def read(
             ds[y_var_name].dims[0],
         )
 
-        # interp/extrap corners
-        descriptor.x_corner = interp_extrap_corner(descriptor.x)
-        descriptor.y_corner = interp_extrap_corner(descriptor.y)
+        # use CF bounds if available, otherwise interp/extrap corners
+        descriptor.x_corner = get_corners_1d(ds, x_var_name)
+        descriptor.y_corner = get_corners_1d(ds, y_var_name)
 
         descriptor.history = add_history(ds=ds)
         return descriptor