Repository: regionmask/regionmask

Issue:
add `flag_values` and `flag_meanings` to 2D masks?
This would allow to indicate which region corresponds to which number via `attrs` and would work nicely with `cf_xarray`.

```python
import cf_xarray
import xarray as xr

da = xr.DataArray(
    [[1, 2, 1], [1, 2, np.NaN]],
    dims=("lat", "lon"),
    attrs={
        "flag_values": [1, 2],
        "flag_meanings": "atlantic_ocean pacific_ocean",
        "standard_name": "region",
    },
  name="mask",
)

print(da.cf)

print(da.cf == "atlantic_ocean")

print(da.cf.isin(["indian_ocean", "pacific_ocean"]))

```

See
- https://cf-xarray.readthedocs.io/en/latest/flags.html
- xarray-contrib/cf-xarray#199
- http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#flags

Open questions:
- do we add `abbrevs` or `names` (probably `abbrevs`)
- does it work for `float` arrays (as I encode "not belonging to any region" with `NaN`)?

Discussion:
codecov-commenter: # [Codecov](https://codecov.io/gh/regionmask/regionmask/pull/361?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) Report
> Merging [#361](https://codecov.io/gh/regionmask/regionmask/pull/361?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) (5b77f29) into [main](https://codecov.io/gh/regionmask/regionmask/commit/e5416cc03027c69a0b1a451860c4028d43dadc21?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) (e5416cc) will **decrease** coverage by `0.05%`.
> The diff coverage is `94.73%`.

```diff
@@            Coverage Diff             @@
##             main     #361      +/-   ##
==========================================
- Coverage   98.64%   98.58%   -0.06%     
==========================================
  Files          15       15              
  Lines        1328     1345      +17     
==========================================
+ Hits         1310     1326      +16     
- Misses         18       19       +1     
```

| Flag | Coverage Δ | |
|---|---|---|
| bare-minimum | `91.43% <94.73%> (+0.03%)` | :arrow_up: |
| min-all-deps | `97.78% <94.73%> (-0.05%)` | :arrow_down: |
| unittests | `98.58% <94.73%> (-0.06%)` | :arrow_down: |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#carryforward-flags-in-the-pull-request-comment) to find out more.

| [Impacted Files](https://codecov.io/gh/regionmask/regionmask/pull/361?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | Coverage Δ | |
|---|---|---|
| [regionmask/core/regions.py](https://codecov.io/gh/regionmask/regionmask/pull/361/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#diff-cmVnaW9ubWFzay9jb3JlL3JlZ2lvbnMucHk=) | `97.47% <92.85%> (-0.38%)` | :arrow_down: |
| [regionmask/core/mask.py](https://codecov.io/gh/regionmask/regionmask/pull/361/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#diff-cmVnaW9ubWFzay9jb3JlL21hc2sucHk=) | `100.00% <100.00%> (ø)` | |
| [regionmask/core/utils.py](https://codecov.io/gh/regionmask/regionmask/pull/361/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#diff-cmVnaW9ubWFzay9jb3JlL3V0aWxzLnB5) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask).

jbusecke: Still have this on my list to test, will try to get to it tomorrow. Sorry for the delay.
mathause: No worries - another look is certainly very welcome.

That being said I went over the code again and from my side this is ready to be merged.

Patch:
diff --git a/ci/requirements/docs.yml b/ci/requirements/docs.yml
index 903c5d7..d729bb8 100644
--- a/ci/requirements/docs.yml
+++ b/ci/requirements/docs.yml
@@ -18,6 +18,7 @@ dependencies:
   - xarray
 # dependencies for the examples
   - cftime
+  - cf_xarray
   - netcdf4
 # for regionmask intake example
   - aiohttp
diff --git a/docs/installation.rst b/docs/installation.rst
index 450d831..1abe396 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -23,6 +23,14 @@ For plotting
 - `cartopy <http://scitools.org.uk/cartopy/>`__ (0.17 or later) for plotting on
   geographical maps.
 
+For parsing flag values
+~~~~~~~~~~~~~~~~~~~~~~~
+
+- `cf_xarray <https://cf-xarray.readthedocs.io/en/latest/>`__ (0.6.1 or later) allows
+  rich comparison of ``mask.cf`` for 2D masks with the abbreviations or names of the
+  regions.
+
+
 For faster masking
 ~~~~~~~~~~~~~~~~~~
 
diff --git a/docs/notebooks/mask_2D.ipynb b/docs/notebooks/mask_2D.ipynb
index b64d809..6240dfe 100644
--- a/docs/notebooks/mask_2D.ipynb
+++ b/docs/notebooks/mask_2D.ipynb
@@ -276,7 +276,89 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We want to select the region 'Central North America'. Thus we first need to find out which number this is:"
+    "### Select a specific region - using cf_xarray"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If [cf_xarray](https://cf-xarray.readthedocs.io/en/latest/index.html) is installed and imported we can select any individual region using the [flags](https://cf-xarray.readthedocs.io/en/latest/flags.html) saved in the attributes (`attrs`) of the mask (new in regionmask 0.10.0):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cf_xarray"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Thus to select the region named \"Central North America\" (abbreviated \"CNA\") we can do:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_CNA = mask.cf == \"CNA\"\n",
+    "\n",
+    "# show a subset\n",
+    "mask_CNA.sel(lat=30, lon=slice(240, 260))"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "metadata": {},
+   "source": [
+    ".. warning::\n",
+    "   ``flag_meanings`` cannot contain spaces - they will be replaced by underscores:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_names = regionmask.defined_regions.srex.mask(airtemps, flag=\"names\")\n",
+    "\n",
+    "try:\n",
+    "    mask_names.cf == \"C. North America\"\n",
+    "except ValueError:\n",
+    "    print(\"error raised\")\n",
+    "\n",
+    "mask_names.cf == \"C. North America\".replace(\" \", \"_\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Select a specific region - without cf_xarray"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We first need to find out which number the region 'C. North America' corresponds to:"
    ]
   },
   {
@@ -289,6 +371,15 @@
     "CNA_index"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_CNA = mask == CNA_index"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -300,7 +391,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "`xarray` provides the handy `where` method ([documentation](https://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html)), that filters elements from this object according to a condition:"
+    "To replace all grid points outside of CNA we use the `where` method of `xr.Dataset` ([documentation](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.where.html)), which filters elements from this object according to a condition:"
    ]
   },
   {
@@ -309,7 +400,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "airtemps_CNA = airtemps.where(mask == CNA_index)"
+    "airtemps_CNA = airtemps.where(mask_CNA)"
    ]
   },
   {
@@ -524,7 +615,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rasm_NAS = rasm.where(mask == regionmask.defined_regions.srex.map_keys(\"NAS\"))"
+    "rasm_NAS = rasm.where(mask.cf == \"NAS\")"
    ]
   },
   {
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 877deb6..17b6bdd 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -22,6 +22,10 @@ Breaking Changes
 Enhancements
 ~~~~~~~~~~~~
 
+- 2D masks (e.g. :py:meth:`Regions.mask`) now contain `flag_values` and `flag_values` as
+  attributes (`attrs`). Together with `cf_xarray <https://cf-xarray.readthedocs.io/en/latest/flags.html>`__
+  these can be leveraged to select single (``mask.cf == "CNA"``) or multiple (``mask.cf.isin``)
+  regions (:pull:`361`, :issue:`346`).
 - The masking functions (e.g. :py:meth:`Regions.mask`) now warn if the `units` of the
   coordinates(``lat.attrs["units"]`` ) are given as "radians" (:issue:`279`).
 
@@ -49,6 +53,8 @@ Internal Changes
 
   - using shapely, pygeos (:pull:`343`), and rasterio (:pull:`345`)
   - in the function to determine points at *exactly* -180°E (or 0°E) and -90°N (:pull:`341`)
+- Use importlib.metadata if available (i.e. for python > 3.7) - should lead to a faster
+  import time for regionmask (:pull:`369`).
 
 .. _whats-new.0.9.0:
 
diff --git a/regionmask/__init__.py b/regionmask/__init__.py
index 0c86de5..3482f47 100644
--- a/regionmask/__init__.py
+++ b/regionmask/__init__.py
@@ -1,16 +1,20 @@
 # flake8: noqa
 
-import os
-
-import pkg_resources
-
 from . import core, defined_regions
 from .core._geopandas import from_geopandas, mask_3D_geopandas, mask_geopandas
 from .core.plot import plot_3D_mask
 from .core.regions import Regions, _OneRegion
 
 try:
-    __version__ = pkg_resources.get_distribution("regionmask").version
+    from importlib.metadata import version as _get_version
+except ImportError:
+    # importlib.metadata not available in python 3.7
+    import pkg_resources
+
+    _get_version = lambda pkg: pkg_resources.get_distribution(pkg).version
+
+try:
+    __version__ = _get_version("regionmask")
 except Exception:
     # Local copy or not installed with setuptools.
     # Disable minimum version checks on downstream libraries.
diff --git a/regionmask/core/mask.py b/regionmask/core/mask.py
index 7260c31..99c2092 100644
--- a/regionmask/core/mask.py
+++ b/regionmask/core/mask.py
@@ -55,7 +55,8 @@ wrap_lon : bool | 180 | 360, optional
       and wrapped to `[-180, 180[` if its maximum is larger than 180.
     - ``180``: Wraps longitude coordinates to `[-180, 180[`
     - ``360``: Wraps longitude coordinates to `[0, 360[`
-{overlap}
+
+{overlap}{flags}
 
 Returns
 -------
@@ -97,6 +98,15 @@ overlap : bool, default: False
     There is (currently) no automatic detection of overlapping regions.
 """
 
+_FLAG_DOCSTRING = """\
+flag : str, default "abbrevs"
+    Indicates if the "abbrevs" (abbreviations) or "names" should be added as
+    `flag_values` and `flag_meanings` to the attributes (`attrs`) of the mask. If None
+    nothing is added. Using cf_xarray these can be used to select single
+    (``mask.cf == "CNA"``) or multiple (``mask.cf.isin``) regions. Note that spaces are
+    replaced by underscores.
+"""
+
 
 def _inject_mask_docstring(is_3D, gp_method):
 
@@ -106,6 +116,7 @@ def _inject_mask_docstring(is_3D, gp_method):
     numbers_doc = _NUMBERS_DOCSTRING if gp_method else ""
     gp_doc = _GP_DOCSTRING if gp_method else ""
     overlap = _OVERLAP_DOCSTRING if (gp_method and is_3D) else ""
+    flags = _FLAG_DOCSTRING if not (gp_method or is_3D) else ""
 
     mask_docstring = _MASK_DOCSTRING_TEMPLATE.format(
         dtype=dtype,
@@ -114,6 +125,7 @@ def _inject_mask_docstring(is_3D, gp_method):
         numbers_doc=numbers_doc,
         gp_doc=gp_doc,
         overlap=overlap,
+        flags=flags,
     )
 
     return mask_docstring
@@ -245,6 +257,8 @@ def _mask_2D(
         msg = "No gridpoint belongs to any region. Returning an all-NaN mask."
         warnings.warn(msg, UserWarning, stacklevel=3)
 
+    mask.attrs = {"standard_name": "region"}
+
     return mask
 
 
@@ -280,6 +294,8 @@ def _mask_3D(
     else:
         mask_3D = _unpack_2D_mask(mask, numbers, drop)
 
+    mask_3D.attrs = {"standard_name": "region"}
+
     return mask_3D
 
 
diff --git a/regionmask/core/regions.py b/regionmask/core/regions.py
index 7a7dd7b..1a93c24 100644
--- a/regionmask/core/regions.py
+++ b/regionmask/core/regions.py
@@ -290,6 +290,7 @@ class Regions:
         lat_name="lat",
         method=None,
         wrap_lon=None,
+        flag="abbrevs",
     ):
 
         if self.overlap:
@@ -298,7 +299,7 @@ class Regions:
                 "Set ``region.overlap = False`` to create a 2D mask anyway."
             )
 
-        return _mask_2D(
+        mask_2D = _mask_2D(
             outlines=self.polygons,
             lon_bounds=self.bounds_global[::2],
             numbers=self.numbers,
@@ -310,6 +311,28 @@ class Regions:
             wrap_lon=wrap_lon,
         )
 
+        if flag not in [None, "abbrevs", "names"]:
+            raise ValueError(
+                f"`flag` must be one of `None`, `'abbrevs'` and `'names'`, found {flag}"
+            )
+
+        if flag is not None:
+            # see http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#flags
+
+            # find detected regions (assign ALL regions?)
+            isnan = np.isnan(mask_2D.values)
+            numbers = np.unique(mask_2D.values[~isnan])
+            numbers = numbers.astype(int)
+
+            flag_meanings = getattr(self[numbers], flag)
+            # TODO: check for invalid characters
+            flag_meanings = " ".join(flag.replace(" ", "_") for flag in flag_meanings)
+
+            mask_2D.attrs["flag_values"] = numbers
+            mask_2D.attrs["flag_meanings"] = flag_meanings
+
+        return mask_2D
+
     mask.__doc__ = _inject_mask_docstring(is_3D=False, gp_method=False)
 
     @_deprecate_positional_args("0.10.0")


Test patch:
diff --git a/regionmask/tests/test_mask.py b/regionmask/tests/test_mask.py
index 831ecb4..c478a0b 100644
--- a/regionmask/tests/test_mask.py
+++ b/regionmask/tests/test_mask.py
@@ -84,7 +84,7 @@ def test_mask(method):
 
     expected = expected_mask_2D()
     result = dummy_region.mask(dummy_ds.lon, dummy_ds.lat, method=method)
-    xr.testing.assert_equal(result, expected)
+    xr.testing.assert_identical(result, expected)
 
 
 @pytest.mark.skipif(has_pygeos, reason="Only errors if pygeos is missing")
@@ -105,7 +105,7 @@ def test_mask_xr_keep_name(method):
 
     result = dummy_region.mask(ds.longitude, ds.latitude, method=method)
 
-    xr.testing.assert_equal(result, expected)
+    xr.testing.assert_identical(result, expected)
 
 
 @pytest.mark.parametrize("method", MASK_METHODS)
@@ -333,6 +333,42 @@ def test_mask_3D_overlap_empty(method):
     assert not result.any()
 
 
+def test_mask_flag():
+
+    expected = expected_mask_2D()
+    result = dummy_region.mask(dummy_ds)
+    xr.testing.assert_identical(result, expected)
+
+    result = dummy_region.mask(dummy_ds, flag="names")
+    expected.attrs["flag_meanings"] = "Region0 Region1"
+    xr.testing.assert_identical(result, expected)
+
+    result = dummy_region.mask(dummy_ds, flag=None)
+    del expected.attrs["flag_meanings"]
+    del expected.attrs["flag_values"]
+
+    xr.testing.assert_identical(result, expected)
+
+
+def test_mask_flag_space():
+
+    r = Regions(dummy_region.polygons, names=["name with space", "another", "last"])
+
+    expected = expected_mask_2D()
+    expected.attrs["flag_meanings"] = "name_with_space another"
+
+    result = r.mask(dummy_ds, flag="names")
+    xr.testing.assert_identical(result, expected)
+
+
+def test_mask_flag_only_found():
+
+    result = dummy_region.mask([0.5], [0.5])
+
+    assert result.attrs["flag_meanings"] == "r0"
+    np.testing.assert_equal(result.attrs["flag_values"], np.array([0]))
+
+
 # ======================================================================
 
 # test 2D array
@@ -456,7 +492,7 @@ def test_mask_3D(drop, method):
     expected = expected_mask_3D(drop)
     result = dummy_region.mask_3D(dummy_ds.lon, dummy_ds.lat, drop=drop, method=method)
 
-    xr.testing.assert_equal(result, expected)
+    xr.testing.assert_identical(result, expected)
 
 
 @pytest.mark.parametrize("method", MASK_METHODS)
@@ -826,6 +862,7 @@ def test_inject_mask_docstring():
     assert "drop :" in result
     assert "geodataframe" in result
     assert "overlap" in result
+    assert "flag" not in result
 
     result = _inject_mask_docstring(False, False)
 
@@ -834,3 +871,4 @@ def test_inject_mask_docstring():
     assert "drop :" not in result
     assert "geodataframe" not in result
     assert "overlap" not in result
+    assert "flag" in result
diff --git a/regionmask/tests/utils.py b/regionmask/tests/utils.py
index 759c1fb..0f2d097 100644
--- a/regionmask/tests/utils.py
+++ b/regionmask/tests/utils.py
@@ -39,7 +39,15 @@ def expected_mask_2D(a=0, b=1, fill=np.NaN, flatten=False, coords=None, dims=Non
     if dims is None:
         dims = ("lat", "lon")
 
-    return xr.DataArray(mask, dims=dims, coords=coords, name="region")
+    flag_values = np.array([a, b])
+    flag_meanings = "r0 r1"
+
+    attrs = {"standard_name": "region"}
+
+    attrs["flag_values"] = flag_values
+    attrs["flag_meanings"] = flag_meanings
+
+    return xr.DataArray(mask, dims=dims, coords=coords, name="mask", attrs=attrs)
 
 
 def expected_mask_3D(drop, coords=None, overlap=False):
@@ -66,8 +74,9 @@ def expected_mask_3D(drop, coords=None, overlap=False):
         }
     )
     dims = ("region",) + ("lat", "lon")
+    attrs = {"standard_name": "region"}
 
-    expected = xr.DataArray(mask, coords=coords, dims=dims)
+    expected = xr.DataArray(mask, coords=coords, dims=dims, name="mask", attrs=attrs)
 
     return expected.drop_sel(region=numbers[-1]) if drop else expected