diff --git a/.github/workflows/build_workflow.yml b/.github/workflows/build_workflow.yml index aac1e6f..d724de7 100644 --- a/.github/workflows/build_workflow.yml +++ b/.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 diff --git a/dev-environment.txt b/dev-environment.txt index d6beb7d..84df991 100644 --- a/dev-environment.txt +++ b/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 diff --git a/docs/conf.py b/docs/conf.py index 6160950..d071ff8 100644 --- a/docs/conf.py +++ b/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), } diff --git a/docs/developers_guide/api.md b/docs/developers_guide/api.md index 57104e7..ce46f13 100644 --- a/docs/developers_guide/api.md +++ b/docs/developers_guide/api.md @@ -39,4 +39,10 @@ Properties Descriptor.vertex_patches Descriptor.transform +Helper Functions +---------------- +.. autosummary:: + :toctree: generated/ + + utils.get_invalid_patches ``` diff --git a/docs/index.md b/docs/index.md index 480669d..819e8e9 100644 --- a/docs/index.md +++ b/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} diff --git a/docs/user_guide/quick_start.md b/docs/user_guide/quick_start.md index becc90c..39fb922 100644 --- a/docs/user_guide/quick_start.md +++ b/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: . + +### 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 diff --git a/docs/user_guide/wrapping.md b/docs/user_guide/wrapping.md new file mode 100644 index 0000000..5ca2db9 --- /dev/null +++ b/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. + + + +## Supported Map Projections for Spherical Meshes + +Currently, the only support map projection are: +- +- +- +- +- +- +- +- +- +- +- +- +- +- +- +- +- + +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). diff --git a/mosaic/datasets.py b/mosaic/datasets.py index de8936e..cf20409 100644 --- a/mosaic/datasets.py +++ b/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 ---------- diff --git a/mosaic/descriptor.py b/mosaic/descriptor.py index 5cfed64..4258057 100644 --- a/mosaic/descriptor.py +++ b/mosaic/descriptor.py @@ -1,5 +1,7 @@ from functools import cached_property +from typing import Literal +import cartopy.crs as ccrs import numpy as np import xarray as xr from cartopy.crs import CRS @@ -19,6 +21,13 @@ "verticesOnCell", "edgesOnVertex"] +SUPPORTED_SPHERICAL_PROJECTIONS = (ccrs._RectangularProjection, + ccrs._WarpedRectangularProjection, + ccrs.Stereographic, + ccrs.Mercator, + ccrs._CylindricalProjection, + ccrs.InterruptedGoodeHomolosine) + def attr_to_bool(attr: str): """ Format attribute strings and return a boolean value """ @@ -105,6 +114,8 @@ def __init__( #: ``projection`` kwargs are provided. self.transform = transform + #: Boolean whether parent mesh is (planar) periodic in at least one dim + self.is_periodic = attr_to_bool(mesh_ds.is_periodic) #: Boolean whether parent mesh is spherical self.is_spherical = attr_to_bool(mesh_ds.on_a_sphere) @@ -120,6 +131,10 @@ def __init__( # if both a projection and transform were provided to the constructor self.projection = projection + # method ensures is periodic, avoiding AttributeErrors if non-periodic + self.x_period = self._parse_period(mesh_ds, "x") + self.y_period = self._parse_period(mesh_ds, "y") + def _create_minimal_dataset(self, ds: Dataset) -> Dataset: """ Create a xarray.Dataset that contains the minimal subset of @@ -190,10 +205,20 @@ def projection(self) -> CRS: @projection.setter def projection(self, projection: CRS) -> None: + # We don't support all map projections for spherical meshes, yet... + if (projection is not None and self.is_spherical and + not isinstance(projection, SUPPORTED_SPHERICAL_PROJECTIONS)): + + raise ValueError(f"Invalid projection: {type(projection).__name__}" + f" is not supported - consider using " + f"a rectangular projection.") + + reprojecting = False # Issue warning if changing the projection after initialization # TODO: Add heuristic size (i.e. ``self.ds.nbytes``) above which the # warning is raised if hasattr(self, "_projection"): + reprojecting = True print(("Reprojecting the descriptor can be inefficient " "for large meshes")) @@ -212,6 +237,14 @@ def projection(self, projection: CRS) -> None: self._projection = projection + # if periods are None (i.e. projection was not set at instantiation) or + # the descriptor is being reprojected; update the periods + if (hasattr(self, "_x_period") and hasattr(self, "_x_period")): + if (not self.x_period and not self.y_period) or reprojecting: + # dummy value b/c `_projection` attr will be used by setters + self.x_period = None + self.y_period = None + @property def latlon(self) -> bool: """ @@ -228,6 +261,57 @@ def latlon(self, value) -> None: self._latlon = value + def _parse_period(self, ds, dim: Literal["x", "y"]): + """ Parse period attribute, return None for non-periodic meshes """ + + attr = f"{dim}_period" + try: + period = float(ds.attrs[attr]) + except KeyError: + period = None + + # in the off chance mesh is periodic but does not have period attribute + if self.is_periodic and attr not in ds.attrs: + raise AttributeError((f"Mesh file: \"{ds.encoding['source']}\"" + f"does not have attribute `{attr}` despite" + f"being a planar periodic mesh.")) + if period == 0.0: + return None + else: + return period + + @property + def x_period(self) -> float | None: + """ Period along x-dimension, is ``None`` for non-periodic meshes """ + return self._x_period + + @x_period.setter + def x_period(self, value) -> None: + # needed to avoid AttributeError for non-periodic meshes + if not (self.is_periodic and self.is_spherical): + self._x_period = None + if not self.is_periodic and self.is_spherical and self.projection: + x_limits = self.projection.x_limits + self._x_period = np.abs(x_limits[1] - x_limits[0]) + else: + self._x_period = value + + @property + def y_period(self) -> float | None: + """ Period along y-dimension, is ``None`` for non-periodic meshes """ + return self._y_period + + @y_period.setter + def y_period(self, value) -> None: + # needed to avoid AttributeError for non-periodic meshes + if not (self.is_periodic and self.is_spherical): + self._y_period = None + if not self.is_periodic and self.is_spherical and self.projection: + y_limits = self.projection.y_limits + self._y_period = np.abs(y_limits[1] - y_limits[0]) + else: + self._y_period = value + @cached_property def cell_patches(self) -> ndarray: """:py:class:`~numpy.ndarray` of patch coordinates for cell centered @@ -244,7 +328,13 @@ def cell_patches(self) -> ndarray: patch. Nodes are ordered counter clockwise around the cell center. """ patches = _compute_cell_patches(self.ds) - patches = self._fix_antimeridian(patches, "Cell") + patches = self._wrap_patches(patches, "Cell") + + # cartopy doesn't handle nans in patches, so store a mask of the + # invalid patches to set the dataarray at those locations to nan. + if self.projection: + self._cell_pole_mask = self._compute_pole_mask("Cell") + return patches @cached_property @@ -263,7 +353,13 @@ def edge_patches(self) -> ndarray: corresponding node will be collapsed to the edge coordinate. """ patches = _compute_edge_patches(self.ds) - patches = self._fix_antimeridian(patches, "Edge") + patches = self._wrap_patches(patches, "Edge") + + # cartopy doesn't handle nans in patches, so store a mask of the + # invalid patches to set the dataarray at those locations to nan. + if self.projection: + self._edge_pole_mask = self._compute_pole_mask("Edge") + return patches @cached_property @@ -289,7 +385,13 @@ def vertex_patches(self) -> ndarray: position. """ patches = _compute_vertex_patches(self.ds) - patches = self._fix_antimeridian(patches, "Vertex") + patches = self._wrap_patches(patches, "Vertex") + + # cartopy doesn't handle nans in patches, so store a mask of the + # invalid patches to set the dataarray at those locations to nan. + if self.projection: + self._vertex_pole_mask = self._compute_pole_mask("Vertex") + return patches def _transform_coordinates(self, projection, transform): @@ -304,70 +406,74 @@ def _transform_coordinates(self, projection, transform): self.ds[f"x{loc}"].values = transformed_coords[:, 0] self.ds[f"y{loc}"].values = transformed_coords[:, 1] - def _fix_antimeridian(self, patches, loc, projection=None): - """Correct vertices of patches that cross the antimeridian. - - NOTE: Can this be a decorator? + def _wrap_patches(self, patches, loc): + """Wrap patches for spherical and planar-periodic meshes """ - # coordinate arrays are transformed at initalization, so using the - # transform size limit, not the projection - if not projection: - projection = self.projection - - # should be able to come up with a default size limit here, or maybe - # it's already an attribute(?) Should also factor in a precomputed - # axis period, as set in the attributes of the input dataset - if projection: - # convert to numpy array to that broadcasting below will work - x_center = np.array(self.ds[f"x{loc}"]) - - # get distance b/w the center and vertices of the patches - # NOTE: using data from masked patches array so that we compute - # mask only corresponds to patches that cross the boundary, - # (i.e. NOT a mask of all invalid cells). May need to be - # carefull about the fillvalue depending on the transform - half_distance = x_center[:, np.newaxis] - patches[..., 0].data - - # get the size limit of the projection; - size_limit = np.abs(projection.x_limits[1] - - projection.x_limits[0]) / (2 * np.sqrt(2)) - - # left and right mask, with same number of dims as the patches - l_mask = (half_distance > size_limit)[..., np.newaxis] - r_mask = (half_distance < -size_limit)[..., np.newaxis] + def _find_boundary_patches(patches, loc, coord, period): """ - # Old approach masks out all patches that cross the antimeridian. - # This is unnessarily restrictive. New approach corrects - # the x-coordinates of vertices that lie outside the projections - # bounds, which isn't perfect either - - patches.mask |= l_mask - patches.mask |= r_mask + Find the patches that cross the periodic boundary and what + direction they cross the boundary (i.e. their ``sign``). This + method assumes the patch centroids are not periodic """ + # get axis index we are inquiring over + axis = 0 if coord == "x" else 1 + # get requested coordinate of patch centroids + center = self.ds[f"{coord}{loc.title()}"].values.reshape(-1, 1) + # get difference b/w centroid and nodes of patches + diff = patches[..., axis] - center + + mask = np.abs(diff) > np.abs(period) / (2. * np.sqrt(2.)) + sign = np.sign(diff) + + return mask, sign + + def _wrap_1D(patches, mask, sign, axis, period): + """Correct patch periodicity along a single dimension""" + patches[..., axis][mask] -= np.sign(sign[mask]) * period - l_boundary_mask = ~np.any(l_mask, axis=1) | l_mask[..., 0] - r_boundary_mask = ~np.any(r_mask, axis=1) | r_mask[..., 0] - # get valid half distances for the patches that cross boundary - l_offset = np.ma.MaskedArray(half_distance, l_boundary_mask) - r_offset = np.ma.MaskedArray(half_distance, r_boundary_mask) - - # For vertices that cross the antimeridian reset the x-coordinate - # of invalid vertex to be the center of the patch plus the - # mean valid half distance. - # - # NOTE: this only fixes patches on the side of plot where they - # cross the antimeridian, leaving an empty zipper like pattern - # mirrored over the y-axis. - patches[..., 0] = np.ma.where( - ~l_mask[..., 0], patches[..., 0], - x_center[:, np.newaxis] + l_offset.mean(1)[..., np.newaxis]) - patches[..., 0] = np.ma.where( - ~r_mask[..., 0], patches[..., 0], - x_center[:, np.newaxis] + r_offset.mean(1)[..., np.newaxis]) + # TODO: clip spherical wrapped patches to projection limits + return patches + + # Stereographic projections do not need wrapping, so exit early + if isinstance(self.projection, ccrs.Stereographic): + return patches + + if self.x_period: + # find the patches that are periodic in x-direction + x_mask, x_sign = _find_boundary_patches( + patches, loc, "x", self.x_period + ) + + if np.any(x_mask): + # using the sign of the difference correct patches x coordinate + patches = _wrap_1D(patches, x_mask, x_sign, 0, self.x_period) + + if self.y_period: + # find the patches that are periodic in y-direction + y_mask, y_sign = _find_boundary_patches( + patches, loc, "y", self.y_period + ) + + if np.any(y_mask): + # using the sign of the difference correct patches y coordinate + patches = _wrap_1D(patches, y_mask, y_sign, 1, self.y_period) return patches + def _compute_pole_mask(self, loc) -> ndarray: + """ """ + limits = self.projection.y_limits + centers = self.ds[f"y{loc.title()}"].values + + # TODO: determine threshold for ``isclose`` computation + at_pole = np.any( + np.isclose(centers.reshape(-1, 1), limits, rtol=1e-2), axis=1 + ) + past_pole = np.abs(centers) > np.abs(limits[1]) + + return (at_pole | past_pole) + def _compute_cell_patches(ds: Dataset) -> ndarray: """Create cell patches (i.e. Primary cells) for an MPAS mesh.""" @@ -437,7 +543,6 @@ def _compute_vertex_patches(ds: Dataset) -> ndarray: # get a mask of active nodes cellMask = cellsOnVertex == -1 edgeMask = edgesOnVertex == -1 - unionMask = cellMask & edgeMask # get the coordinates needed to patch construction xCell = ds.xCell.values @@ -457,14 +562,15 @@ def _compute_vertex_patches(ds: Dataset) -> ndarray: nodes[:, 1::2, 1] = np.where(cellMask, yVertex, yCell[cellsOnVertex]) # ------------------------------------------------------------------------- - # NOTE: While the condition below probably only applies to the final edge - # node we apply it to all, since the conditions above ensure the - # patches will still be created correctly + # NOTE: The condition below will only be true for meshes run through the + # MPAS mesh converter after culling. A bug in the converter alters + # the ordering of edges, causing problems for vertex patches + # + # If final cell and edge nodes are missing collapse both back to first edge + # Ensures patches encompasses the full kite area and are properly closed. # ------------------------------------------------------------------------- - # if cell and edge missing collapse edge node to the first edge. - # Because edges lead the vertices this ensures the patch encompasses - # the full kite area and is properly closed. - nodes[:, ::2, 0] = np.where(unionMask, nodes[:, 0:1, 0], nodes[:, ::2, 0]) - nodes[:, ::2, 1] = np.where(unionMask, nodes[:, 0:1, 1], nodes[:, ::2, 1]) + condition = (cellMask & edgeMask)[:, -1:] + nodes[:, 4:, 0] = np.where(condition, nodes[:, 0:1, 0], nodes[:, 4:, 0]) + nodes[:, 4:, 1] = np.where(condition, nodes[:, 0:1, 1], nodes[:, 4:, 1]) return nodes diff --git a/mosaic/polypcolor.py b/mosaic/polypcolor.py index d17d9f1..17111e7 100644 --- a/mosaic/polypcolor.py +++ b/mosaic/polypcolor.py @@ -1,3 +1,4 @@ +import numpy as np from cartopy.mpl.geoaxes import GeoAxes from matplotlib.axes import Axes from matplotlib.collections import PolyCollection @@ -48,12 +49,18 @@ def polypcolor( if "nCells" in c.dims: verts = descriptor.cell_patches + if descriptor.projection and np.any(descriptor._cell_pole_mask): + c = c.where(~descriptor._cell_pole_mask, np.nan) elif "nEdges" in c.dims: verts = descriptor.edge_patches + if descriptor.projection and np.any(descriptor._edge_pole_mask): + c = c.where(~descriptor._edge_pole_mask, np.nan) elif "nVertices" in c.dims: verts = descriptor.vertex_patches + if descriptor.projection and np.any(descriptor._vertex_pole_mask): + c = c.where(~descriptor._vertex_pole_mask, np.nan) transform = descriptor.transform diff --git a/mosaic/utils.py b/mosaic/utils.py new file mode 100644 index 0000000..5db86e8 --- /dev/null +++ b/mosaic/utils.py @@ -0,0 +1,44 @@ +from typing import Literal + +import numpy as np +from shapely import LineString, is_valid + +from mosaic.descriptor import Descriptor + + +def get_invalid_patches( + descriptor: Descriptor, location: Literal["cell", "edge", "vertex"] +) -> None | np.ndarray: + """Helper function to identify problematic patches. + + Returns the indices of the problematic patches as determined by + :py:func:`shapely.is_valid`. + + Parameters + ---------- + descriptor : :py:class:`Descriptor` + The ``Descriptor`` object you want to check the patches of + location : string + The patch location to check: "cell" or "edge" or "vertex" + + Returns + ------- + :py:class:`~numpy.ndarray` or None + Indices of invalid patches. If no patches are invalid then returns None + """ + # extract the patches + patches = descriptor.__getattribute__(f"{location}_patches") + # get the pole mask b/c we know patches will be invalid there + pole_mask = descriptor.__getattribute__(f"_{location}_pole_mask") + # convert the patches to a list of shapely geometries + geoms = [LineString(patch) for patch in patches[~pole_mask]] + # check if the shapely geometries are valid + valid = is_valid(geoms) + + if np.all(valid): + # no invalid patches, so return None + return None + else: + # return the indices of the invalid patches + index_array = np.arange(patches.shape[0]) + return index_array[~pole_mask][~valid] diff --git a/tests/test_spherical.py b/tests/test_spherical.py new file mode 100644 index 0000000..ec2dedb --- /dev/null +++ b/tests/test_spherical.py @@ -0,0 +1,85 @@ +import cartopy.crs as ccrs +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import pytest +from shapely import LinearRing, is_valid + +import mosaic + +matplotlib.use('Agg', force=True) + +SUPPORTED_PROJECTIONS = [ + ccrs.PlateCarree(), + ccrs.LambertCylindrical(), + ccrs.Mercator(), + ccrs.Miller(), + ccrs.Robinson(), + ccrs.Stereographic(), + ccrs.RotatedPole(), + ccrs.InterruptedGoodeHomolosine(), + ccrs.EckertI(), + ccrs.EckertII(), + ccrs.EckertIII(), + ccrs.EckertIV(), + ccrs.EckertV(), + ccrs.EckertVI(), + ccrs.EqualEarth(), + ccrs.NorthPolarStereo(), + ccrs.SouthPolarStereo(), +] + + +def id_func(projection): + return type(projection).__name__ + + +class TestSphericalWrapping: + + ds = mosaic.datasets.open_dataset("QU.240km") + + @pytest.fixture(scope="module", params=SUPPORTED_PROJECTIONS, ids=id_func) + def setup_descriptor(self, request): + descriptor = mosaic.Descriptor(self.ds, request.param, ccrs.Geodetic()) + return descriptor + + @pytest.mark.timeout(30) + @pytest.mark.parametrize("patch", ["Cell", "Edge", "Vertex"]) + @pytest.mark.filterwarnings("ignore:numpy.ndarray size changed") + def test_timeout(self, setup_descriptor, tmp_path, patch): + # do the test setup with the parameterized projection + descriptor = setup_descriptor + + projection = descriptor.projection + + # get the projection name + proj_name = type(projection).__name__ + + # setup with figure with the parameterized projection + fig, ax = plt.subplots(subplot_kw=dict(projection=projection)) + + # get the appropriate dataarray for the parameterized patch location + da = self.ds[f"indexTo{patch}ID"] + + # just testing that this doesn't hang, not for correctness + mosaic.polypcolor(ax, descriptor, da, antialiaseds=True) + + # save the figure so that patches are rendered + fig.savefig(f"{tmp_path}/{proj_name}-{patch}.png") + plt.close() + + @pytest.mark.parametrize("patch", ["Cell", "Edge", "Vertex"]) + def test_valid_patches(self, setup_descriptor, patch): + # do the test setup with the parameterized projection + descriptor = setup_descriptor + + # extract the patches + patches = descriptor.__getattribute__(f"{patch.lower()}_patches") + # get the pole mask b/c we know patches will be invalid there + pole_mask = descriptor.__getattribute__(f"_{patch.lower()}_pole_mask") + + # convert the patches to a list of shapely geometries + geoms = [LinearRing(patch) for patch in patches[~pole_mask]] + + # assert that all the patches are valid + assert np.all(is_valid(geoms))