raster: Read raster for mask from env variable

doc/development/style_guide.md

@@ -463,15 +463,21 @@ The `--overwrite` flag can be globally enabled by setting the environment variab
 #### Mask
 
 GRASS GIS has a global mask managed by the _r.mask_ tool and represented by a
-raster called MASK. Raster tools called as a subprocess will automatically
+raster called MASK by default. Raster tools called as a subprocess will automatically
 respect the globally set mask when reading the data. For outputs, respecting of
 the mask is optional.
 
-Tools **should not set or remove the global mask**. If the tool cannot avoid
-setting the mask internally, it should check for presence of the mask and fail
-if the mask is present. The tools should not remove and later restore the
-original mask because that creates confusing behavior for interactive use and
-breaks parallel processing.
+Tools should generally respect the global mask set by a user. If the mask set by the
+user is not respected by a tool, the exact behavior should be described in the
+documentation. On the other hand, ignoring mask is usually the desired behavior
+for import tools which corresponds with the mask being applied only when reading
+existing raster data in a project.
+
+Tools **should not set or remove the global mask** to prevent unintended
+behavior during interactive sessions and to maintain parallel processing
+integrity. If a tool requires a mask for its operation, it should implement
+a temporary mask using _MaskManager_ in Python or by setting the `GRASS_MASK`
+environment variable.
 
 Generally, any mask behavior should be documented unless it is the standard case
 where masked cells do not participate in the computation and are represented as
@@ -577,6 +583,76 @@ gs.run_command("r.slope.aspect", elevation=input_raster, slope=slope, env=env)
 This approach makes the computational region completely safe for parallel
 processes as no region-related files are modified.
 
+#### Changing raster mask
+
+The _MaskManager_ in Python API provides a way for tools to change, or possibly
+to ignore, a raster mask for part of the computation.
+
+In the following example, _MaskManager_ modifies the global system environment
+for the tool (aka _os.environ_) so that custom mask can be applied:
+
+```python
+# Previously user-set mask applies here (if any).
+gs.run_command("r.slope.aspect", elevation=input_raster, aspect=aspect)
+
+with gs.MaskManager():
+    # Only the mask we set here will apply.
+    gs.run_command("r.mask", raster=mask_raster)
+    gs.run_command("r.slope.aspect", elevation=input_raster, slope=slope)
+# Mask is disabled and the mask raster is removed at the end of the with block.
+
+# Previously user-set mask applies here again.
+```
+
+Because tools should generally respect the provided mask, the mask in a tool
+should act as an additional mask. This can be achieved when preparing the new
+mask raster using a tool which reads an existing raster:
+
+```python
+# Here we create an initial mask by creating a raster from vector,
+# but that does not use mask.
+gs.run_command(
+    "v.to.rast", input=input_vector, where="name == 'Town'", output=town_boundary
+)
+# So, we use a raster algebra expression. Mask will be applied if set
+# because in the expression, we are reading an existing raster.
+gs.mapcalc(f"{raster_mask} = {town_boundary}")
+
+with gs.MaskManager():
+    gs.run_command("r.mask", raster=mask_raster)
+    # Both user mask and the town_boundary are used here.
+    gs.run_command("r.slope.aspect", elevation=input_raster, slope=slope)
+```
+
+To disable the mask, which may be needed in processing steps of import tool,
+we can do:
+
+```python
+# Mask applies here if set.
+gs.run_command("r.slope.aspect", elevation=input_raster, aspect=aspect)
+
+with gs.MaskManager():
+    # No mask was set in this context, so the tool runs without a mask.
+    gs.run_command("r.slope.aspect", elevation=input_raster, slope=slope)
+
+# Mask applies again.
+```
+
+If needed, tools can implement optional support for a user-set raster mask by
+passing or not passing the current name of a mask obtained from _r.mask.status_
+and by preparing the internal mask raster beforehand with the user mask active.
+
+If different subprocesses, running in parallel, use different masks,
+it is best to create mask rasters beforehand (to avoid limitations of _r.mask_ and
+the underlying _r.reclass_ tool). The name of the mask raster can then be passed to
+the manager:
+
+```python
+env = os.environ.copy()
+with gs.MaskManager(mask_name=mask_raster, env=env):
+    gs.run_command("r.slope.aspect", elevation=input_raster, slope=slope, env=env)
+```
+
 #### Temporary Maps
 
 Using temporary maps is preferred over using temporary mapsets. This follows the

doc/examples/notebooks/parallelization_tutorial.ipynb

@@ -459,6 +459,73 @@
     "\n",
     "Alternatively, a PostgreSQL or another database backend can be used for attributes to offload the parallel writing to the database system."
    ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "df444278",
+   "metadata": {},
+   "source": [
+    "#### Safely modifying raster mask in a single mapset"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "a8d52c4c",
+   "metadata": {},
+   "source": [
+    "Mask is by default specified per-mapset and shared by all the processes. Additionally, *r.mask* is using *r.reclass* in the background which may cause issues if the mask is derived from the same base map in parallel. The use of *MaskManager* in the following example allows each process to use a different raster. The raster is used directly as a mask to avoid the need to use *r.mask*.\n",
+    "\n",
+    "The following code derives basins (watersheds) based on a threshold from the digital elevation model. Then, for each basin, it computes the topographic index with *r.topidx*."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f72dcca",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%writefile example.py\n",
+    "elevation = \"elev_state_500m\"\n",
+    "gs.run_command(\"g.region\", raster=elevation)\n",
+    "gs.run_command(\"r.watershed\", elevation=elevation, basin=\"basins\", threshold=10000)\n",
+    "\n",
+    "cats = gs.parse_command(\"r.describe\", map=\"basins\", flags=\"1n\", format=\"json\")[\"values\"]\n",
+    "\n",
+    "def topidx(cat):\n",
+    "    # Define output name and mask name.\n",
+    "    output = f\"topidx_{cat}\"\n",
+    "    basin = f\"basin_{cat}\"\n",
+    "    # Extract subwatershed by category into separate raster, creating\n",
+    "    # a 0-or-1 mask which has no NULLs (although NULLs in mask are allowed).\n",
+    "    gs.mapcalc(f\"{basin} = if(isnull(basins), 0, if(basins == {cat}, 1, 0))\")\n",
+    "    # Create a copy of the environment for this process before modifying it.\n",
+    "    env = os.environ.copy()\n",
+    "    # Set the computational region to match non-null area in the new raster.\n",
+    "    env[\"GRASS_REGION\"] = gs.region_env(raster=\"basins\", zoom=basin, env=env)\n",
+    "    # Use mask context manager to specify which raster to use as a mask\n",
+    "    # and pass the environment we are using.\n",
+    "    with gs.MaskManager(mask_name=basin, env=env):\n",
+    "        # Run actual computation with active mask.\n",
+    "        gs.run_command(\"r.topidx\", input=elevation, output=output, env=env)\n",
+    "    return output\n",
+    "\n",
+    "with Pool(processes=4) as pool:\n",
+    "    outputs = pool.map(topidx, cats)\n",
+    "    print(outputs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8c6cfb3b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%run example.py"
+   ]
   }
  ],
  "metadata": {

lib/init/variables.html

@@ -512,6 +512,11 @@ <h3>List of selected internal GRASS environment variables</h3>
     This allows programs such as the GUI to run external commands on an
     alternate region without having to modify the WIND file then change it
     back afterwards.</dd>
+
+  <dt>GRASS_MASK</dt>
+  <dd>[libgis]<br>
+    use the raster map specified by name as mask, instead of a raster called
+    MASK in the current mapset.</dd>
 </dl>
 
 <h2>List of selected GRASS gisenv variables</h2>

lib/raster/auto_mask.c

@@ -3,14 +3,13 @@
  *
  * \brief Raster Library - Auto masking routines.
  *
- * (C) 2001-2008 by the GRASS Development Team
+ * (C) 2001-2024 by Vaclav Petras and the GRASS Development Team
  *
  * This program is free software under the GNU General Public License
  * (>=v2). Read the file COPYING that comes with GRASS for details.
  *
  * \author GRASS GIS Development Team
- *
- * \date 1999-2008
+ * \author Vaclav Petras (environmental variable and refactoring)
  */
 
 #include <stdlib.h>
@@ -31,46 +30,51 @@
  * \return 0 if mask unset or unavailable
  * \return 1 if mask set and available and ready to use
  */
-
 int Rast__check_for_auto_masking(void)
 {
     struct Cell_head cellhd;
 
     Rast__init();
 
     /* if mask is switched off (-2) return -2
-       if R__.auto_mask is not set (-1) or set (>=0) recheck the MASK */
+       if R__.auto_mask is not set (-1) or set (>=0) recheck the mask */
 
+    // TODO: This needs to be documented or modified accordingly.
     if (R__.auto_mask < -1)
         return R__.auto_mask;
 
     /* if(R__.mask_fd > 0) G_free (R__.mask_buf); */
 
-    /* look for the existence of the MASK file */
-    R__.auto_mask = (G_find_raster("MASK", G_mapset()) != 0);
+    /* Decide between default mask name and env var specified one. */
+    char *mask_name = Rast_mask_name();
+    char *mask_mapset = "";
+
+    /* Check for the existence of the mask raster. */
+    R__.auto_mask = (G_find_raster2(mask_name, mask_mapset) != 0);
 
     if (R__.auto_mask <= 0)
         return 0;
 
-    /* check MASK projection/zone against current region */
-    Rast_get_cellhd("MASK", G_mapset(), &cellhd);
+    /* Check mask raster projection/zone against current region */
+    Rast_get_cellhd(mask_name, mask_mapset, &cellhd);
     if (cellhd.zone != G_zone() || cellhd.proj != G_projection()) {
         R__.auto_mask = 0;
         return 0;
     }
 
     if (R__.mask_fd >= 0)
         Rast_unopen(R__.mask_fd);
-    R__.mask_fd = Rast__open_old("MASK", G_mapset());
+    R__.mask_fd = Rast__open_old(mask_name, mask_mapset);
     if (R__.mask_fd < 0) {
         R__.auto_mask = 0;
-        G_warning(_("Unable to open automatic MASK file"));
+        G_warning(_("Unable to open automatic mask <%s>"), mask_name);
         return 0;
     }
 
     /*    R__.mask_buf = Rast_allocate_c_buf(); */
 
     R__.auto_mask = 1;
+    G_free(mask_name);
 
     return 1;
 }

lib/raster/get_row.c

@@ -764,7 +764,7 @@ void Rast_get_d_row_nomask(int fd, DCELL *buf, int row)
  *            two particular types check the functions).
  *    - Step 4:  read or simmulate null value row and zero out cells
  * corresponding to null value cells. The masked out cells are set to null when
- * the mask exists. (the MASK is taken care of by null values (if the null file
+ * the mask exists. (the mask is taken care of by null values (if the null file
  * doesn't exist for this map, then the null row is simulated by assuming that
  * all zero are nulls *** in case of Rast_get_row() and assuming that all data
  * is valid in case of G_get_f/d_raster_row(). In case of deprecated function
@@ -1089,7 +1089,7 @@ static void embed_nulls(int fd, void *buf, int row, RASTER_MAP_TYPE map_type,
 
    Read or simulate null value row and set the cells corresponding
    to null value to 1. The masked out cells are set to null when the
-   mask exists. (the MASK is taken care of by null values
+   mask exists. (the mask is taken care of by null values
    (if the null file doesn't exist for this map, then the null row
    is simulated by assuming that all zeros in raster map are nulls.
    Also all masked out cells become nulls.