Merge pull request #23 from andrewdnolan/periodic_planar

Add support for planar periodic meshes by adding wrapping functionality that's applicable to both planar periodic and spherical meshes.

Author: andrewdnolan

.github/workflows/build_workflow.yml

@@ -110,7 +110,7 @@ jobs:
           conda activate mosaic_dev
           pip check
           python -c "import mosaic"
-          #pytest tests
+          pytest tests -v
 
       - if: ${{ steps.skip_check.outputs.should_skip != 'true' }}
         name: Build Sphinx Docs

dev-environment.txt

@@ -8,6 +8,7 @@ numpy
 pooch
 pyproj
 scipy
+shapely
 tqdm
 xarray
 
@@ -17,6 +18,7 @@ setuptools >=60
 # Linting and tesing
 pip
 pytest
+pytest-timeout
 isort
 flake8
 mypy

docs/conf.py

@@ -64,6 +64,7 @@ def setup(app):
         'matplotlib': ('https://matplotlib.org/stable', None),
         'numpy': ('https://numpy.org/doc/stable', None),
         'python': ('https://docs.python.org/3/', None),
+        'shapely': ('https://shapely.readthedocs.io/en/stable/', None),
         'xarray': ('https://xarray.pydata.org/en/stable', None),
         }
 

docs/developers_guide/api.md

@@ -39,4 +39,10 @@ Properties
     Descriptor.vertex_patches
     Descriptor.transform
 
+Helper Functions
+----------------
+.. autosummary::
+    :toctree: generated/
+
+    utils.get_invalid_patches
 ```

docs/index.md

@@ -12,6 +12,7 @@ Currently `mosaic` only supports [MPAS](https://mpas-dev.github.io/) meshes, but
 :caption: User Guide:
 
 user_guide/quick_start
+user_guide/wrapping
 ```
 
 ```{toctree}

docs/user_guide/quick_start.md

@@ -35,7 +35,7 @@ ds = mosaic.datasets.open_dataset("QU.240km")
 # define a map projection for our figure
 projection = ccrs.InterruptedGoodeHomolosine()
 # define the transform that describes our dataset
-transform = ccrs.PlateCarree()
+transform = ccrs.Geodetic()
 
 # create the figure and a GeoAxis 
 fig, ax = plt.subplots(1, 1, figsize=(9,7), facecolor="w",
@@ -54,13 +54,16 @@ ax.coastlines()
 fig.colorbar(collection, fraction=0.1, shrink=0.5, label="Cell Index");
 ```
 
-### Planar Non-Periodic 
+For more information about how spherical meshes are handled and a list of supported
+map projections, see: <project:#wrapping>.
+
+### Planar Non-Periodic
 
 In this case the underlying coordinate arrays (i.e. `xCell/yCell`)
 correspond to a South Polar Stereographic projection, which is also the map projection we
 want to us. Therefore, the `projection` and the `transform` will be equivalent
 for this example. When instantiating the `mosaic.Descriptor` object we have to 
-careful to set `use_latlon=False` to ensure the `xCell`/`yCell` coordinate
+be careful to set `use_latlon=False` to ensure the `xCell`/`yCell` coordinate
 arrays are parsed (c.f. `lonCell`/`latCell`). 
 
 ```{code-cell} ipython3
@@ -101,13 +104,13 @@ fig.colorbar(collection, fraction=0.1, label="Thickness [m]");
 In the case where we do not know what projection the coordinate arrays of the
 mesh correspond to we can use the `lonCell`/`latCell` coordinates and `mosaic`
 will handle the transformation to the requested map projection under the hood.
-In this scenario the `transform` should correspond to `ccrs.PlateCarree()`
+In this scenario the `transform` should correspond to `ccrs.Geodetic()`
 and `use_latlon=True` must be set in the `mosaic.Descriptor` object
 instantiation. Nearly all the lines would be the same as the above example,
 with the exception of the transform definition: 
 ```python 
 # define the transform that describes our dataset
-transform = ccrs.PlateCarree()
+transform = ccrs.Geodetic()
 ```
 and the `mosaic.Descriptor` instantiation: 
 ```python

docs/user_guide/wrapping.md

@@ -0,0 +1,78 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+---
+
+# Periodic Mesh Support
+
+We currently make the simplifying assumption that spherical meshes are a
+special instance of periodic meshes, which are periodic across the antimeridian.
+We support both planar periodic and map projected spherical meshes, using the
+same methods, by assuming that the period (in the x and y directions
+respectively) is constant. This assumption is valid for planar periodic meshes
+and for some map projections, but falls apart for certain map projection where
+the antimeridian does not follow a strait line. Therefore we only support a
+subset of `cartopy` projections, which are listed below. Future work will
+develop a more elaborate method of dealing with spherical mesh periodicity,
+which in turn will expand the list supported map projections.
+
+For patches that cross a periodic boundary we simply correct the coordinates to
+remove the periodicity, which enables plotting. Future work will mirror the
+patches across the periodic boundary, so as to more accurately demonstrate the
+periodic nature of the mesh.
+
+<!-- 
+## Planar Periodic Meshes
+
+```{code-cell} ipython3
+---
+mystnb:
+    remove_code_source: true
+---
+import mosaic
+import matplotlib.pyplot as plt
+
+# download and read the mesh from lcrc
+ds = mosaic.datasets.open_dataset("doubly_periodic_4x4")
+
+# create the figure and a GeoAxis 
+fig, ax = plt.subplots(constrained_layout=True,)
+
+descriptor = mosaic.Descriptor(ds)
+
+pc = mosaic.polypcolor(
+    ax, descriptor, ds.indexToCellID, alpha=0.8, antialiaseds=True, ec="k"
+)
+
+ax.scatter(descriptor.ds.xCell, descriptor.ds.yCell, c='k')
+ax.scatter(*descriptor.cell_patches.T, c='tab:blue', marker='^')
+ax.scatter(ds.xVertex, ds.yVertex, ec='tab:orange', fc='none', marker='o', s=5.)
+ax.set_aspect('equal')
+```
+-->
+
+## Supported Map Projections for Spherical Meshes
+
+Currently, the only support map projection are:
+- <inv:#*.PlateCarree>
+- <inv:#*.LambertCylindrical>
+- <inv:#*.Mercator>
+- <inv:#*.Miller>
+- <inv:#*.Robinson>
+- <inv:#*.Stereographic>
+- <inv:#*.RotatedPole>
+- <inv:#*.InterruptedGoodeHomolosine>
+- <inv:#*.EckertI>
+- <inv:#*.EckertII>
+- <inv:#*.EckertIII>
+- <inv:#*.EckertIV>
+- <inv:#*.EckertV>
+- <inv:#*.EckertVI>
+- <inv:#*.EqualEarth>
+- <inv:#*.NorthPolarStereo>
+- <inv:#*.SouthPolarStereo>
+
+It is important to note that planer (non-periodic) meshes are not limited to
+this list of map projections and can choose from the full list of `cartopy`
+[projections](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html).

mosaic/datasets.py

@@ -26,7 +26,12 @@
     "mpasli.AIS8to30": {
         "lcrc_path": "inputdata/glc/mpasli/mpas.ais8to30km/ais_8to30km.20221027.nc",
         "sha256_hash": "sha256:932a1989ff8e51223413ef3ff0056d6737a1fc7f53e440359884a567a93413d2"
-    }
+    },
+
+    "doubly_periodic_4x4": {
+        "lcrc_path": "mpas_standalonedata/mpas-ocean/mesh_database/doubly_periodic_1920km_7680x7680km.151124.nc",
+        "sha256_hash": "sha256:5409d760845fb682ec56e30d9c6aa6dfe16b5d0e0e74f5da989cdaddbf4303c7"
+    },
 }
 
 # create a parsable registry for pooch from human friendly one
@@ -55,8 +60,9 @@ def open_dataset(
 
     * ``"QU.960km"`` : Quasi-uniform spherical mesh, with approximately 960km horizontal resolution
     * ``"QU.240km"`` : Quasi-uniform spherical mesh, with approximately 240km horizontal resolution
-    * ``"mpaso.IcoswISC30E3r5"`` : Icosahedral 30 km MPAS-Ocean mesh with ice shelf cavaties
+    * ``"mpaso.IcoswISC30E3r5"`` : Icosahedral 30 km MPAS-Ocean mesh with ice shelf cavities
     * ``"mpasli.AIS8to30"`` : 8-30 km resolution planar non-periodic MALI mesh of Antarctica
+    * ``"doubly_periodic_4x4"``: Doubly periodic planar mesh that is four cells wide in both the x and y dimensions.
 
     Parameters
     ----------