Repository: regionmask/regionmask

Issue:
warn on overlapping regions?
There are several ways to do this. My suggestions is to go via the area:
```python
total_area = sum(p.area() for p in polygons)
flattened_area = shapely.ops.unary_union(polygons).area

overlap = (total_area - flattened_area) / flattened_area * 100

if overlap > overlap_threshold:
    warnings.warn("overlap")

```
This allows to calculate a relative overlap. One difficulty is that regions might very slightly overlap (see ar6 regions) & we don't want to warn in this case. Questions/ problems:

- The warning might be thrown for a subset of regions but not for all regions (because the total area changes)
- Should this be done in `Regions.__init__` or in `mask`?
- `Regions.__init__`:
  - only done once at instantiation time (+)
  - `overlap_threshold` does not round-trip in `to_geodataframe` (-)
- `mask`: 
  - no roundtrip issue (+)
  - done every time a mask is created (-) [could be partially mitigated by caching it]
  - adds another keyword argument to `mask` (-)
  - we cannot set a region-specific threshold


This is relevant for #244 & #228

TODO: check overlap fraction for all existing regions.

Discussion:
codecov-commenter: ## [Codecov](https://app.codecov.io/gh/regionmask/regionmask/pull/439?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) Report
> Merging [#439](https://app.codecov.io/gh/regionmask/regionmask/pull/439?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) (9d2238e) into [main](https://app.codecov.io/gh/regionmask/regionmask/commit/ae1817082296c2b056730b57f8d48f90fdbc5499?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) (ae18170) will **decrease** coverage by `0.06%`.
> The diff coverage is `93.75%`.

:exclamation: Your organization needs to install the [Codecov GitHub app](https://github.com/apps/codecov/installations/select_target) to enable full functionality.

```diff
@@            Coverage Diff             @@
##             main     #439      +/-   ##
==========================================
- Coverage   98.10%   98.05%   -0.06%     
==========================================
  Files          18       18              
  Lines        1478     1490      +12     
==========================================
+ Hits         1450     1461      +11     
- Misses         28       29       +1     
```

| [Flag](https://app.codecov.io/gh/regionmask/regionmask/pull/439/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | Coverage Δ | |
|---|---|---|
| [bare-minimum](https://app.codecov.io/gh/regionmask/regionmask/pull/439/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | `90.53% <93.75%> (+<0.01%)` | :arrow_up: |
| [min-all-deps](https://app.codecov.io/gh/regionmask/regionmask/pull/439/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | `95.83% <93.75%> (-0.04%)` | :arrow_down: |
| [shapely_lt_2](https://app.codecov.io/gh/regionmask/regionmask/pull/439/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | `?` | |
| [unittests](https://app.codecov.io/gh/regionmask/regionmask/pull/439/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask) | `98.05% <93.75%> (-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.

| [Files Changed](https://app.codecov.io/gh/regionmask/regionmask/pull/439?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://app.codecov.io/gh/regionmask/regionmask/pull/439?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#diff-cmVnaW9ubWFzay9jb3JlL3JlZ2lvbnMucHk=) | `97.53% <ø> (ø)` | |
| [regionmask/core/mask.py](https://app.codecov.io/gh/regionmask/regionmask/pull/439?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=regionmask#diff-cmVnaW9ubWFzay9jb3JlL21hc2sucHk=) | `99.10% <93.75%> (-0.28%)` | :arrow_down: |


mathause: I think this is ready (the failures are not related).

The only open question is: should I warn when `mask_3D` detects overlapping regions? Because this is technically a breaking change.

The warning could be:
```python
warning.warn(
    "Detected overlapping regions. From v0.11.0 these are correctly considered in regionmask, "
    "even when `overlap` is not explicitly set. "
    "To suppress this warning, set `overlap=True` (or `False`)."
)
```

(I think I would remove the warning again some versions later)

Patch:
diff --git a/docs/notebooks/overlap.ipynb b/docs/notebooks/overlap.ipynb
index 5eca01b..650b9d5 100644
--- a/docs/notebooks/overlap.ipynb
+++ b/docs/notebooks/overlap.ipynb
@@ -24,7 +24,7 @@
    "source": [
     "# Overlapping regions\n",
     "\n",
-    "Two or more on regions can share the same area - they overlap, as for example region 3 and 4 of the [PRUDENCE regions](../defined_scientific.html#prudence-regions) This notebook illustrates how overlapping regions can be treated in regionmask."
+    "Two or more on regions can share the same area - they overlap. One example are region 3 and 4 of the [PRUDENCE regions](../defined_scientific.html#prudence-regions). This notebook illustrates how overlapping regions can be treated in regionmask."
    ]
   },
   {
@@ -38,29 +38,41 @@
    "cell_type": "raw",
    "metadata": {},
    "source": [
-    ".. warning::\n",
-    "   Per default regionmask assumes regions are *not* overlapping!"
+    ".. important::\n",
+    "   Starting with v0.11.0 regionmask checks if gridpoints belong to more than one region (unless ``overlap=False``). Thus, it creates \n",
+    "   the correct mask (or raises an error) per default."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Thus, when creating your own `Regions` you need to tell regionmask if they are overlapping.\n",
+    "## In detail"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When creating your own `Regions` you can tell regionmask if they are overlapping:\n",
     "\n",
     "```python\n",
-    "region = regionmask.Regions(..., overlap=True)\n",
-    "region = regionmask.from_geopandas(..., overlap=True)\n",
-    "```"
+    "region = regionmask.Regions(..., overlap=overlap)\n",
+    "region = regionmask.from_geopandas(..., overlap=overlap)\n",
+    "```\n",
+    "\n",
+    "where ``overlap`` must be one of ``None`` (default), ``True`` or ``False``."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If you have two overlapping regions and `overlap=False` regionmask will _silently_ assign the gridpoints of the overlapping regions to the one with the higher number, e.g., region 4 for PRUDENCE (this may change in a future version).\n",
-    "\n",
-    "Note that `overlap` is correctly defined in `regionmask.defined_regions`."
+    "- For `overlap=None` regionmask checks if any grid point belongs to more than one region. If so, 2D masks raise and error (because overlapping regions cannot be represented by a 2D mask). `mask_3D(...)` returns the correct mask.\n",
+    "- For `overlap=True` 2D masks raise an error and ``mask_3D(...)`` returns the correct mask. \n",
+    "- If you have two overlapping regions and `overlap=False` regionmask will _silently_ assign the gridpoints of the overlapping regions to the one with the higher number.\n",
+    "- Setting ``overlap`` to ``True`` or ``False`` is (slightly) faster than the default of ``None``.\n",
+    "- `overlap` is correctly set for regions in `regionmask.defined_regions`."
    ]
   },
   {
@@ -100,16 +112,29 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import xarray as xr\n",
-    "import numpy as np\n",
-    "\n",
     "import cartopy.crs as ccrs\n",
-    "import matplotlib.pyplot as plt\n",
-    "\n",
+    "import matplotlib.patheffects as pe\n",
+    "import numpy as np\n",
     "from matplotlib import colors as mplc\n",
-    "from shapely.geometry import Polygon\n",
-    "\n",
-    "import matplotlib.patheffects as pe"
+    "from shapely.geometry import box"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Define grid:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_US = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(\n",
+    "    *(-160, -29, 2), *(76, 13, -2)\n",
+    ")"
    ]
   },
   {
@@ -141,7 +166,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def plot_region_vh(mask):\n",
+    "def plot_region_vh(mask, region):\n",
     "\n",
     "    fg = mask.plot(\n",
     "        subplot_kws=dict(projection=ccrs.PlateCarree()),\n",
@@ -153,9 +178,9 @@
     "        lw=0.25,\n",
     "    )\n",
     "\n",
-    "    for ax in fg.axes.flatten():\n",
-    "        region_vh[[0]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#6a3d9a\"))\n",
-    "        region_vh[[1]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#ff7f00\"))\n",
+    "    for ax in fg.axs.flatten():\n",
+    "        region[[0]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#6a3d9a\"))\n",
+    "        region[[1]].plot(ax=ax, add_label=False, line_kws=dict(color=\"#ff7f00\"))\n",
     "\n",
     "        ax.set_extent([-105, -75, 25, 55], ccrs.PlateCarree())\n",
     "        ax.plot(\n",
@@ -167,7 +192,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Define the polygons:"
+    "Define overlapping polygons:"
    ]
   },
   {
@@ -176,17 +201,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "coords_v = np.array([[-90.0, 50.0], [-90.0, 28.0], [-100.0, 28.0], [-100.0, 50.0]])\n",
-    "coords_h = np.array([[-80.0, 50.0], [-80.0, 40.0], [-100.0, 40.0], [-100.0, 50.0]])"
+    "coords_v = box(-100, 28, -90, 50)\n",
+    "coords_h = box(-100, 40, -80, 50)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Default behavior (`overlap=False`)\n",
+    "### Default behavior (`overlap=None`)\n",
     "\n",
-    "Fe first test what happens if we keep the default value of `overlap=False`:"
+    "First test what happens if we keep the default value of `overlap=None`:"
    ]
   },
   {
@@ -211,10 +236,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_US = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(\n",
-    "    *(-160, -29, 2), *(76, 13, -2)\n",
-    ")\n",
-    "\n",
     "mask_vh = region_vh.mask_3D(ds_US)"
    ]
   },
@@ -231,15 +252,71 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_region_vh(mask_vh)"
+    "plot_region_vh(mask_vh, region_vh)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The small gray points show the gridpoint center and the vertical and horizontal lines are the gridpoint boundaries. The colored rectangles are the two regions. Here regionmask detects that certain gridpoints belong to two regions and corretly assigns them."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Setting `overlap=False`\n",
+    "\n",
+    "We test what happens for `overlap=False`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "region_vh_overlap_false = regionmask.Regions([coords_v, coords_h], overlap=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Create a mask"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_vh_overlap_false = region_vh_overlap_false.mask_3D(ds_US)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot the masked regions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_region_vh(mask_vh_overlap_false, region_vh_overlap_false)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The small gray points show the gridpoint center and the vertical and horizontal lines are the gridpoint boundaries. The colored rectangles are the two regions. The vertical region has the number 1 and the horizontal region the number 2.\n",
-    "We can see that only the gridpoints in the lower part of the vertical (magenta) region were assigned to it. All gridpoints of the overlapping part are now assigned to the horizontal (orange) region. As mentioned the gridpoints are assigned to the region with the higher number By switching the order of the regions you could have the common gridpoints assigned to the vertical region."
+    "We can see that only the gridpoints in the lower part of the vertical (magenta) region were assigned to it. All gridpoints of the overlapping part are now assigned to the horizontal (orange) region - the gridpoints are assigned to the region with the higher number. By switching the order of the regions you could have the common gridpoints assigned to the vertical region."
    ]
   },
   {
@@ -248,7 +325,7 @@
    "source": [
     "### Setting `overlap=True`\n",
     "\n",
-    "As mentioned regionmask assumes regions are not overlapping, so you need to pass `overlap=True` to the constructor:"
+    "Passing `overlap=True` has the same effect as the default:"
    ]
   },
   {
@@ -291,14 +368,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_region_vh(mask_overlap)"
+    "plot_region_vh(mask_overlap, region_vh)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now the gridpoints in the overlapping part are assigned to both regions."
+    "And the gridpoints in the overlapping part are assigned to both regions."
    ]
   },
   {
@@ -362,8 +439,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "lon = np.arange(-12, 33, 0.5)\n",
-    "lat = np.arange(72, 33, -0.5)\n",
+    "lon = np.arange(-12, 33, 0.25)\n",
+    "lat = np.arange(72, 33, -0.25)\n",
     "\n",
     "mask_prudence = prudence.mask_3D(lon, lat)"
    ]
@@ -385,7 +462,7 @@
     ")\n",
     "\n",
     "\n",
-    "for ax in fg.axes.flatten():\n",
+    "for ax in fg.axs.flatten():\n",
     "    regionmask.defined_regions.prudence.plot(\n",
     "        ax=ax, add_label=False, resolution=\"50m\", line_kws=dict(lw=0.75)\n",
     "    )"
@@ -415,7 +492,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.10.11"
   }
  },
  "nbformat": 4,
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 31584a3..5f7168d 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -42,6 +42,8 @@ Breaking Changes
 Enhancements
 ~~~~~~~~~~~~
 
+- regionmask now checks if regions are overlapping (unless ``overlap=False`` is explicitly set) - check
+  the documentation on :doc:`overlapping regions<notebooks/overlap>` for details (:pull:`439`).
 - Can now pass the ``use_cf`` parameter to :py:func:`mask_geopandas` and :py:func:`mask_3D_geopandas`.
   This could also be counted as a bug fix as these functions could not control the behavior
   of finding the coords otherwise (:pull:`427`).
diff --git a/regionmask/core/_geopandas.py b/regionmask/core/_geopandas.py
index 9417d73..dadbfa6 100644
--- a/regionmask/core/_geopandas.py
+++ b/regionmask/core/_geopandas.py
@@ -189,7 +189,7 @@ def _from_geopandas(
     )
 
 
-def _prepare_gdf_for_mask(geodataframe, method, numbers):
+def _prepare_gdf_for_mask(geodataframe, numbers):
 
     from geopandas import GeoDataFrame, GeoSeries
 
@@ -208,6 +208,9 @@ def _prepare_gdf_for_mask(geodataframe, method, numbers):
     return polygons, numbers
 
 
+# TODO: switch order of use_cf and overlap once the deprecation is finished
+
+
 @_deprecate_positional_args("0.10.0")
 def mask_geopandas(
     geodataframe,
@@ -220,11 +223,17 @@ def mask_geopandas(
     method=None,
     wrap_lon=None,
     use_cf=None,
+    overlap=None,
 ):
 
-    polygons, numbers = _prepare_gdf_for_mask(
-        geodataframe, method=method, numbers=numbers
-    )
+    if overlap:
+        raise ValueError(
+            "Creating a 2D mask with overlapping regions yields wrong results. "
+            "Please use ``mask_3D_geopandas(...)`` instead. "
+            "To create a 2D mask anyway, set ``overlap=False``."
+        )
+
+    polygons, numbers = _prepare_gdf_for_mask(geodataframe, numbers=numbers)
 
     return _mask_2D(
         polygons=polygons,
@@ -235,11 +244,12 @@ def mask_geopandas(
         lat_name=lat_name,
         method=method,
         wrap_lon=wrap_lon,
+        overlap=overlap,
         use_cf=use_cf,
     )
 
 
-mask_geopandas.__doc__ = _inject_mask_docstring(is_3D=False, gp_method=True)
+mask_geopandas.__doc__ = _inject_mask_docstring(is_3D=False, is_gpd=True)
 
 
 @_deprecate_positional_args("0.10.0")
@@ -254,13 +264,11 @@ def mask_3D_geopandas(
     numbers=None,
     method=None,
     wrap_lon=None,
-    overlap=False,
     use_cf=None,
+    overlap=None,
 ):
 
-    polygons, numbers = _prepare_gdf_for_mask(
-        geodataframe, method=method, numbers=numbers
-    )
+    polygons, numbers = _prepare_gdf_for_mask(geodataframe, numbers=numbers)
 
     mask_3D = _mask_3D(
         polygons=polygons,
@@ -272,11 +280,11 @@ def mask_3D_geopandas(
         lat_name=lat_name,
         method=method,
         wrap_lon=wrap_lon,
-        as_3D=overlap,
+        overlap=overlap,
         use_cf=use_cf,
     )
 
     return mask_3D
 
 
-mask_3D_geopandas.__doc__ = _inject_mask_docstring(is_3D=True, gp_method=True)
+mask_3D_geopandas.__doc__ = _inject_mask_docstring(is_3D=True, is_gpd=True)
diff --git a/regionmask/core/mask.py b/regionmask/core/mask.py
index 3cc87e0..357ff80 100644
--- a/regionmask/core/mask.py
+++ b/regionmask/core/mask.py
@@ -108,14 +108,17 @@ drop : boolean, default: True
 """
 
 _OVERLAP_DOCSTRING = """\
-overlap : bool, default: False
-    Indicates if (some of) the regions overlap. If True ``mask_3D_geopandas`` will
-    ensure overlapping regions are correctly assigned to grid points.
-
-    If False (default) assumes non-overlapping regions. Grid points will silently be
-    assigned to the region with the higher number (this may change in a future version).
-
-    There is (currently) no automatic detection of overlapping regions.
+overlap : bool | None, default: None
+    Indicates if (some of) the regions overlap.
+
+    - If True ``mask_3D_geopandas`` ensures overlapping regions are correctly assigned
+      to grid points, while ``mask_geopandas`` raises an Error (because overlapping
+      regions cannot be represented by a 2 dimensional mask).
+    - If False assumes non-overlapping regions. Grid points are silently assigned to the
+      region with the higher number.
+    - If None (default) checks if any gridpoint belongs to more than one region.
+      If this is the case ``mask_3D_geopandas`` correctly assigns them and ``mask_geopandas``
+      raises an Error.
 
 """
 
@@ -130,15 +133,15 @@ flag : str, default: "abbrevs"
 """
 
 
-def _inject_mask_docstring(is_3D, gp_method):
+def _inject_mask_docstring(*, is_3D, is_gpd):
 
     dtype = "boolean" if is_3D else "float"
     nd = "3D" if is_3D else "2D"
     drop_doc = _DROP_DOCSTRING if is_3D else ""
-    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 ""
+    numbers_doc = _NUMBERS_DOCSTRING if is_gpd else ""
+    gp_doc = _GP_DOCSTRING if is_gpd else ""
+    overlap = _OVERLAP_DOCSTRING if is_gpd else ""
+    flags = _FLAG_DOCSTRING if not (is_gpd or is_3D) else ""
     see_also = "mask" if is_3D else "mask_3D"
 
     mask_docstring = _MASK_DOCSTRING_TEMPLATE.format(
@@ -296,8 +299,21 @@ def _mask_2D(
     method=None,
     wrap_lon=None,
     use_cf=None,
+    overlap=None,
 ):
 
+    # NOTE: this is already checked in Regions.mask, and mask_geopandas
+    # double check here if this method is ever made public
+    # if overlap:
+    #     raise ValueError(
+    #         "Creating a 2D mask with overlapping regions yields wrong results. "
+    #         "Please use ``region.mask_3D(...)`` instead. "
+    #         "To create a 2D mask anyway, set ``overlap=False``."
+    #     )
+
+    # if as_3D is not explicitely given - set it to True
+    as_3D = overlap is None
+
     mask = _mask(
         polygons=polygons,
         numbers=numbers,
@@ -308,8 +324,13 @@ def _mask_2D(
         method=method,
         wrap_lon=wrap_lon,
         use_cf=use_cf,
+        as_3D=as_3D,
     )
 
+    # only happens for (overlap == None)
+    if as_3D:
+        mask = _3D_to_2D_mask(mask, numbers)
+
     if np.all(np.isnan(mask)):
         msg = "No gridpoint belongs to any region. Returning an all-NaN mask."
         warnings.warn(msg, UserWarning, stacklevel=3)
@@ -329,10 +350,12 @@ def _mask_3D(
     lat_name=None,
     method=None,
     wrap_lon=None,
-    as_3D=False,
+    overlap=None,
     use_cf=None,
 ):
 
+    as_3D = overlap or overlap is None
+
     mask = _mask(
         polygons=polygons,
         numbers=numbers,
@@ -347,16 +370,25 @@ def _mask_3D(
     )
 
     if as_3D:
-        mask_3D = _unpack_3D_mask(mask, numbers, drop)
+        mask_3D = _3D_to_3D_mask(mask, numbers, drop)
     else:
-        mask_3D = _unpack_2D_mask(mask, numbers, drop)
+        mask_3D = _2D_to_3D_mask(mask, numbers, drop)
+
+    if overlap is None and (mask_3D.sum("region") > 1).any():
+        warnings.warn(
+            "Detected overlapping regions. As of v0.11.0 these are correctly taken into "
+            "account. Note, however, that a different mask is returned than with older "
+            "versions of regionmask. To suppress this warning, set `overlap=True` (to "
+            "restore the old, incorrect, behaviour, set `overlap=False`)."
+        )
 
     mask_3D.attrs = {"standard_name": "region"}
 
     return mask_3D
 
 
-def _unpack_2D_mask(mask, numbers, drop):
+def _2D_to_3D_mask(mask, numbers, drop):
+    # TODO: unify with _3D_to_3D_mask
 
     isnan = np.isnan(mask.values)
 
@@ -392,7 +424,8 @@ def _unpack_2D_mask(mask, numbers, drop):
     return mask_3D
 
 
-def _unpack_3D_mask(mask_3D, numbers, drop):
+def _3D_to_3D_mask(mask_3D, numbers, drop):
+    # TODO: unify with _2D_to_3D_mask
 
     any_masked = mask_3D.any(mask_3D.dims[1:])
 
@@ -425,6 +458,35 @@ def _unpack_3D_mask(mask_3D, numbers, drop):
     return mask_3D
 
 
+def _3D_to_2D_mask(mask_3D, numbers):
+
+    # NOTE: very similar to regionmask.core.utils.flatten_3D_mask
+
+    is_masked = mask_3D.sum("region")
+
+    if (is_masked > 1).any():
+        raise ValueError(
+            "Found overlapping regions for ``overlap=None``. Please create a 3D mask. "
+            "You may want to explicitely set ``overlap`` to ``True`` or ``False``."
+        )
+
+    # reshape because region is the first dim
+    numbers = np.asarray(numbers)
+
+    # I transform so it broadcasts (mask_3D can acually be 2D for unstructured grids,
+    # so reshape would need some logic)
+    mask_2D = (mask_3D.T * numbers).T.sum("region")
+
+    # mask all gridpoints not belonging to any region
+
+    # older xarray versions do not have `keep_attrs` argument (needed to keep the name)
+    # mask_2D = xr.where(is_masked, mask_2D, np.nan, keep_attrs=True)
+
+    mask_2D.values = np.where(is_masked.values, mask_2D, np.nan)
+
+    return mask_2D
+
+
 def _determine_method(lon, lat):
     """find method to be used -> prefers faster methods"""
 
diff --git a/regionmask/core/regions.py b/regionmask/core/regions.py
index db6afe5..920ddf2 100644
--- a/regionmask/core/regions.py
+++ b/regionmask/core/regions.py
@@ -36,16 +36,17 @@ class Regions:
         Name of the collection of regions. Default: "unnamed"
     source : string, optional
         Source of the region definitions. Default: "".
-    overlap : bool, default: False
-        Indicates if (some of) the regions overlap. If True ``mask_3D`` will ensure
-        overlapping regions are correctly assigned to grid points while ``mask`` will
-        error (because overlapping regions cannot be represented by a 2D mask).
-
-        If False (default) assumes non-overlapping regions. Grid points will
-        silently be assigned to the region with the higher number (this may change
-        in a future version of regionmask).
-
-        There is (currently) no automatic detection of overlapping regions.
+    overlap : bool | None, default: None
+        Indicates if (some of) the regions overlap.
+
+        - If True ``Regions.mask_3D`` ensures overlapping regions are correctly assigned
+          to grid points, while ``Regionsmask`` raises an Error  (because overlapping
+          regions cannot be represented by a 2D mask).
+        - If False assumes non-overlapping regions. Grid points are silently assigned to
+          the region with the higher number.
+        - If None (default) checks if any gridpoint belongs to more than one region.
+          If this is the case ``Regions.mask_3D`` correctly assigns them and
+          ``Regions.mask`` raises an Error.
 
     Examples
     --------
@@ -65,7 +66,7 @@ class Regions:
     >>> r = Regions(outlines, numbers, names, abbrevs, name)
     >>> r
     <regionmask.Regions 'Example'>
-    overlap:  False
+    overlap:  None
     <BLANKLINE>
     Regions:
     0 uSq1 Unit Square1
@@ -85,7 +86,7 @@ class Regions:
     >>> r = Regions(outlines, numbers, names, abbrevs, name)
     >>> r
     <regionmask.Regions 'Example'>
-    overlap:  False
+    overlap:  None
     <BLANKLINE>
     Regions:
     1 uSq1 Unit Square1
@@ -97,7 +98,7 @@ class Regions:
     >>> r = Regions(outlines)
     >>> r
     <regionmask.Regions 'unnamed'>
-    overlap:  False
+    overlap:  None
     <BLANKLINE>
     Regions:
     0 r0 Region0
@@ -114,7 +115,7 @@ class Regions:
         abbrevs=None,
         name="unnamed",
         source=None,
-        overlap=False,
+        overlap=None,
     ):
 
         if isinstance(outlines, (np.ndarray, Polygon, MultiPolygon)):
@@ -345,6 +346,7 @@ class Regions:
             method=method,
             wrap_lon=wrap_lon,
             use_cf=use_cf,
+            overlap=self.overlap,
         )
 
         if flag not in [None, "abbrevs", "names"]:
@@ -369,7 +371,7 @@ class Regions:
 
         return mask_2D
 
-    mask.__doc__ = _inject_mask_docstring(is_3D=False, gp_method=False)
+    mask.__doc__ = _inject_mask_docstring(is_3D=False, is_gpd=False)
 
     @_deprecate_positional_args("0.10.0")
     def mask_3D(
@@ -395,7 +397,7 @@ class Regions:
             lat_name=lat_name,
             method=method,
             wrap_lon=wrap_lon,
-            as_3D=self.overlap,
+            overlap=self.overlap,
             use_cf=use_cf,
         )
 
@@ -409,7 +411,7 @@ class Regions:
 
         return mask_3D
 
-    mask_3D.__doc__ = _inject_mask_docstring(is_3D=True, gp_method=False)
+    mask_3D.__doc__ = _inject_mask_docstring(is_3D=True, is_gpd=False)
 
     def to_dataframe(self):
         """Convert this region into a pandas.DataFrame, excluding polygons.


Test patch:
diff --git a/regionmask/tests/test_geopandas.py b/regionmask/tests/test_geopandas.py
index 1df8750..a06f5a9 100644
--- a/regionmask/tests/test_geopandas.py
+++ b/regionmask/tests/test_geopandas.py
@@ -278,9 +278,11 @@ def test_mask_3D_geopandas_numbers(geodataframe_clean, drop, method):
     xr.testing.assert_equal(result, expected)
 
 
+@pytest.mark.filterwarnings("ignore:Detected overlapping regions")
 @pytest.mark.parametrize("drop", [True, False])
 @pytest.mark.parametrize("method", ["rasterize", "shapely"])
-def test_mask_3D_overlap_geopandas_numbers(geodataframe_clean_overlap, drop, method):
+@pytest.mark.parametrize("overlap", [True, None])
+def test_mask_3D_geopandas_overlap(geodataframe_clean_overlap, drop, method, overlap):
 
     expected = expected_mask_3D(drop, overlap=True).drop_vars(["names", "abbrevs"])
     result = mask_3D_geopandas(
@@ -290,12 +292,52 @@ def test_mask_3D_overlap_geopandas_numbers(geodataframe_clean_overlap, drop, met
         drop=drop,
         method=method,
         numbers="numbers",
-        overlap=True,
+        overlap=overlap,
     )
 
     xr.testing.assert_equal(result, expected)
 
 
+@pytest.mark.parametrize("drop", [True, False])
+@pytest.mark.parametrize("method", ["rasterize", "shapely"])
+def test_mask_3D_geopandas_overlap_false(geodataframe_clean_overlap, drop, method):
+
+    # NOTE: we need to adapt the expected da
+    expected = expected_mask_3D(drop, overlap=True).drop_vars(["names", "abbrevs"])
+    result = mask_3D_geopandas(
+        geodataframe_clean_overlap,
+        dummy_ds.lon,
+        dummy_ds.lat,
+        drop=drop,
+        method=method,
+        numbers="numbers",
+        overlap=False,
+    )
+
+    if drop:
+        expected = expected.drop_isel(region=0)
+    else:
+        expected[{"region": 0}] = False
+
+    xr.testing.assert_equal(result, expected)
+
+
+@pytest.mark.parametrize("method", ["rasterize", "shapely"])
+@pytest.mark.parametrize("overlap", [None, True])
+def test_mask_2D_geopandas_overlap(geodataframe_clean_overlap, method, overlap):
+
+    # overlap None / True raises a different error
+    with pytest.raises(ValueError):
+        mask_geopandas(
+            geodataframe_clean_overlap,
+            dummy_ds.lon,
+            dummy_ds.lat,
+            method=method,
+            numbers="numbers",
+            overlap=overlap,
+        )
+
+
 @pytest.mark.parametrize("drop", [True, False])
 def test_mask_3D_geopandas_warns_empty(geodataframe_clean, drop):
 
diff --git a/regionmask/tests/test_mask.py b/regionmask/tests/test_mask.py
index ec7b92c..d987e8a 100644
--- a/regionmask/tests/test_mask.py
+++ b/regionmask/tests/test_mask.py
@@ -319,20 +319,49 @@ def test_mask_wrong_method():
 @pytest.mark.parametrize("method", MASK_METHODS)
 def test_mask_2D_overlap_error(method):
 
+    match = "Creating a 2D mask with overlapping regions yields wrong results"
+    with pytest.raises(ValueError, match=match):
+        dummy_region_overlap.mask(dummy_ds, method=method)
+
+
+@pytest.mark.parametrize("method", MASK_METHODS)
+def test_mask_2D_overlap_false(method):
+
     # make a copy to ensure dummy_region_overlap.overlap is not overwritten
     region = copy.copy(dummy_region_overlap)
     expected = expected_mask_2D(a=1, b=2)
 
-    match = "Creating a 2D mask with overlapping regions yields wrong results"
-    with pytest.raises(ValueError, match=match):
-        region.mask(dummy_ds, method=method)
-
     region.overlap = False
     result = region.mask(dummy_ds, method=method)
 
     xr.testing.assert_equal(result, expected)
 
 
+@pytest.mark.parametrize("method", MASK_METHODS)
+def test_mask_2D_overlap_none(method):
+
+    # make a copy to ensure dummy_region_overlap.overlap is not overwritten
+    region = copy.copy(dummy_region_overlap)
+
+    region.overlap = None
+
+    with pytest.raises(
+        ValueError, match="Found overlapping regions for ``overlap=None``"
+    ):
+        region.mask(dummy_ds, method=method)
+
+
+@pytest.mark.parametrize("method", MASK_METHODS)
+def test_mask_2D_overlap_default(method):
+
+    region = Regions(dummy_region_overlap.polygons)
+
+    with pytest.raises(
+        ValueError, match="Found overlapping regions for ``overlap=None``"
+    ):
+        region.mask(dummy_ds, method=method)
+
+
 @pytest.mark.parametrize("drop", [True, False])
 @pytest.mark.parametrize("method", MASK_METHODS)
 def test_mask_3D_overlap(drop, method):
@@ -343,6 +372,35 @@ def test_mask_3D_overlap(drop, method):
     xr.testing.assert_equal(result, expected)
 
 
+@pytest.mark.parametrize("drop", [True, False])
+@pytest.mark.parametrize("method", MASK_METHODS)
+def test_mask_3D_overlap_one(drop, method):
+
+    # make a copy to ensure dummy_region_overlap.overlap is not overwritten
+    region = copy.copy(dummy_region_overlap)
+
+    region.overlap = None
+
+    expected = expected_mask_3D(drop=drop, overlap=True)
+    result = dummy_region_overlap.mask_3D(dummy_ds, drop=drop, method=method)
+
+    xr.testing.assert_equal(result, expected)
+
+
+@pytest.mark.parametrize("drop", [True, False])
+@pytest.mark.parametrize("method", MASK_METHODS)
+def test_mask_3D_overlap_default(drop, method):
+
+    region = Regions(dummy_region_overlap.polygons)
+
+    expected = expected_mask_3D(drop=drop, overlap=True)
+
+    with pytest.warns(UserWarning, match="Detected overlapping regions"):
+        result = region.mask_3D(dummy_ds, drop=drop, method=method)
+
+    xr.testing.assert_equal(result, expected)
+
+
 @pytest.mark.parametrize("method", MASK_METHODS)
 def test_mask_3D_overlap_empty(method):
 
@@ -938,7 +996,7 @@ def test_mask_whole_grid_overlap(method, outline, lon):
 
 def test_inject_mask_docstring():
 
-    result = _inject_mask_docstring(True, True)
+    result = _inject_mask_docstring(is_3D=True, is_gpd=True)
 
     assert "3D" in result
     assert "2D" not in result
@@ -948,7 +1006,7 @@ def test_inject_mask_docstring():
     assert "overlap" in result
     assert "flag" not in result
 
-    result = _inject_mask_docstring(False, False)
+    result = _inject_mask_docstring(is_3D=False, is_gpd=False)
 
     assert "2D" in result
     assert "float" in result