01_intro_to_GRASS_parallelization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7e79c218-1c52-49a8-9d67-e0617cd89648",
   "metadata": {},
   "source": [
    "# Introduction to Parallelization in GRASS GIS\n",
    "The goal of parallelization is to speed up computation by using multiple cores. This notebooks introduces parallelization concepts, existing parallelized tools, and approaches to parallelizing user scripts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69b6064a-679d-49c4-83dd-adc7d3253f4f",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">Throughout the workshop, we will be using shell syntax (Bash shell specifically) and Python. Some things are faster to write in Bash, and some in Python. By default, a cell in this notebook will interpret the code as Python, unless we use <code>%%bash</code> magic to let the notebook know a particular cell contains Bash syntax.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae2cd0a2-2544-46dc-8e3f-cd892a69ef9e",
   "metadata": {},
   "source": [
    "First, let's download and unzip the [prepared dataset](https://doi.org/10.5281/zenodo.8206463) by executing the cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "598fe960-cdd3-45a1-a803-ed3a8a108a3a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "sh download_dataset.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebe80f9a-b1ad-495b-bb15-0c2ae3d4d9e9",
   "metadata": {},
   "source": [
    "Let's start GRASS to run examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee7e7db5-ed77-4611-b539-644e29db2ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import GRASS packages\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "session = gj.init(\"gis_week_2023/part_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38aff518-c001-4553-94dc-c1f695663796",
   "metadata": {},
   "source": [
    "The examples will run in *gis_week_2023* project (previously called location) that contains sample data in the *PERMANENT* mapset (subproject). To keep things organized, this part of the workshop will run in a currently empty mapset *part_1*:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023/aaf6ec7c2b249db68cb79b6d1042247a80381c44/data_structure_embed.svg\" alt=\"GRASS data structure\" width=\"400\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd31203-dbdf-41a1-807c-810947e63b26",
   "metadata": {},
   "source": [
    "### Measuring time\n",
    "It's useful to be able to measure how much time a particular computation takes. We will be using `time` command in shell and `%%timeit` [IPython magic](https://ipython.readthedocs.io/en/stable/interactive/magics.html) for Python cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6490ebfa-7d0b-413b-ac0d-eaebd8a63378",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "time sleep 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d200ac-f5fd-4d3e-b8d0-5981754dc485",
   "metadata": {},
   "source": [
    "Command time gives you 3 different numbers, but we are usually interested in the first one (real time), which is the real-life time it takes for a process to run from start to finish. The other two measure CPU time and you may see the \"user\" time can be larger than the \"real\" time for parallel tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0116139a-b387-4a2a-b5c4-d529b9f516b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit  -n10 -r3\n",
    "import time\n",
    "\n",
    "def function(x):\n",
    "    time.sleep(x)\n",
    "\n",
    "function(0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41678131-68de-48c2-9bdc-7074de13abad",
   "metadata": {},
   "source": [
    "%%timeit magic will return elapsed time executing the Python cell. It typically runs the code multiple times to get a more accurate estimate. You can modify the number of loops (-n) and numer of repeated runs (-r)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbeb3c49-a413-4df8-975f-958db6c91275",
   "metadata": {},
   "source": [
    "### Running GRASS tools in Bash and Python\n",
    "GRASS tools can be executed from Bash and Python (using the [grass.script package](https://grass.osgeo.org/grass-stable/manuals/libpython/script_intro.html) which is part of [GRASS GIS Python API](https://grass.osgeo.org/grass-stable/manuals/libpython/index.html). In the following example, you can extract Wake county from Triangle counties with [v.extract](https://grass.osgeo.org/grass-stable/manuals/v.extract.html) using the following Bash and Python syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8fc5bab-cdc1-4a56-9b17-9acae1e3ffd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "v.extract input=counties output=wake where=\"name == 'Wake'\" --q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a965c09-36bd-4d38-bdbb-0635e665d8d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.extract\", input=\"counties\", output=\"wake\", where=\"name == 'Wake'\", quiet=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bdfb158-7fad-43e9-a8c3-014808492df6",
   "metadata": {},
   "source": [
    "Display Wake county using [grass.jupyter.InteractiveMap](https://grass.osgeo.org/grass83/manuals/libpython/grass.jupyter.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a94b6ac-281b-46df-b7bd-619b92cc44a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "map = gj.InteractiveMap()\n",
    "map.add_vector(name=\"wake\")\n",
    "map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59ff3fed-5763-498d-a1de-9d48e6f85dfc",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Using parallelized tools in GRASS GIS\n",
    "\n",
    "There are many tools in GRASS GIS that are already parallelized ([see the list](https://grass.osgeo.org/grass-stable/manuals/keywords.html#parallel)). Many tools in GRASS Addons are parallelized as well.\n",
    "\n",
    "Generally, there are two types of implementation in GRASS GIS.\n",
    "Multithreading in C tools:\n",
    "   * Threads have low overhead, so they can be spawned more efficiently.\n",
    "   * Tools use OpenMP API. One of the advantages of OpenMP for software distribution is that code works (compiles and runs in serial) also without OpenMP library present on the system.\n",
    "   * Memory is shared, so programmer needs to be cautious about race conditions (e.g., writing into the same variable).\n",
    "   \n",
    "Multiprocessing in Python tools:\n",
    "   * There are multiple ways to implement it, typically tools use `subprocess` and `multiprocessing` package.\n",
    "   * Python tools are often wrappers around GRASS tools implemented in C. For example, tool [r.sun.daily](https://grass.osgeo.org/grass-stable/manuals/addons/r.sun.daily.html) runs [r.sun](https://grass.osgeo.org/grass-stable/manuals/r.sun.html) for multiple days in parallel.\n",
    "   \n",
    "Parallelized tools have `nprocs` parameter to specify number of cores to use. For C tools using OpenMP, GRASS GIS needs to be compiled with OpenMP support to take advantage of it. Both implementations work well on a single machine, but can't be scaled to a distributed system. Scaling to a distributed system is covered at the end of this tutorial.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "198846a3-58a5-4a43-b648-dd81ad5903d9",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">The following examples assume an active GRASS session, meaning GRASS is started (or initialized in the notebook) in a specific project and mapset and tools are executed interactively. Towards the end of this notebook, we will cover running tools non-interactively (sometimes called batch mode).</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bb71615-6e2b-4c72-8974-8d5256613c3e",
   "metadata": {},
   "source": [
    "### Example\n",
    "Let's run our first GRASS tool using multiple cores. Tool [r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html) computes moving window analysis, in this case we will smooth a digital elevation model. r.neighbors use OpenMP for parallelization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d635cf05-1b13-4930-899b-04c04fc57f76",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "time r.neighbors --q input=DEM output=DEM_smooth size=15 method=average nprocs=1\n",
    "time r.neighbors --q input=DEM output=DEM_smooth size=15 method=average nprocs=2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6964b1a6-502c-42d7-be14-a89ef7125ff2",
   "metadata": {},
   "source": [
    "The speedup (ratio of serial time to parallelized time with N cores) typically does not increase linearly with  the number of cores and parallel efficiency (speedup / N cores) decreases when adding cores. See, for example, the [benchmarks for r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html#performance). This behavior is due to the serial parts of the code (see [Amdahl's law](https://en.wikipedia.org/wiki/Amdahl%27s_law)) and computation overhead."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc1c4453-79e2-4c42-b5d8-7a3181732a5b",
   "metadata": {},
   "source": [
    "Let's use GRASS benchmarking library and matplotlib to look at this behavior in more detail. We will compare the time and parallel efficiency for different number of cores and different neighborhood window size of 3 and 9 cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e02f66e-8f96-4baf-8d79-0b10d339386b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import grass.benchmark as bm\n",
    "from grass.pygrass.modules import Module\n",
    "from subprocess import DEVNULL\n",
    "\n",
    "results = []\n",
    "module = Module(\n",
    "        \"r.neighbors\",\n",
    "        input=\"DEM\",\n",
    "        output=\"benchmark\",\n",
    "        size=3,\n",
    "        run_=False,\n",
    "        stdout_=DEVNULL,\n",
    "        overwrite=True,\n",
    "    )\n",
    "results.append(bm.benchmark_nprocs(module, label=\"size 3\", max_nprocs=4, repeat=1))\n",
    "module.inputs.size = 15\n",
    "results.append(bm.benchmark_nprocs(module, label=\"size 15\", max_nprocs=4, repeat=1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c8cd246-deea-410c-b36f-5c586063a8eb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.subplot(1, 3, 1)\n",
    "plt.plot(results[0].nprocs, results[0].times, \"-o\", label=\"size_3\")\n",
    "plt.plot(results[1].nprocs, results[1].times, \"-^\", label=\"size_15\")\n",
    "plt.ylim(bottom=0)\n",
    "plt.xlabel(\"Number of cores\")\n",
    "plt.ylabel(\"Time [s]\")\n",
    "plt.legend()\n",
    "# plt.tight_layout()\n",
    "\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.plot(results[0].nprocs, np.array(results[0].times[0]) / results[0].times, \"-o\", label=\"size_3\")\n",
    "plt.plot(results[1].nprocs, np.array(results[1].times[0]) / results[1].times, \"-o\", label=\"size_15\")\n",
    "# plt.plot(results[1].nprocs, results[1].efficiency, \"-^\", label=\"size_9\")\n",
    "# plt.ylim(bottom=0)\n",
    "plt.xlabel(\"Number of cores\")\n",
    "plt.ylabel(\"Parallel speedup\")\n",
    "# plt.tight_layout()\n",
    "\n",
    "\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.plot(results[0].nprocs, results[0].efficiency, \"-o\", label=\"size_3\")\n",
    "plt.plot(results[1].nprocs, results[1].efficiency, \"-^\", label=\"size_15\")\n",
    "plt.ylim(bottom=0)\n",
    "plt.xlabel(\"Number of cores\")\n",
    "plt.ylabel(\"Parallel efficiency\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2934b6ea-2ac1-4080-b959-9a4054b56213",
   "metadata": {},
   "source": [
    "While computation with larger window size takes longer, it has better parallel efficiency, because a larger proportion of time is spent in the parallel part of the computation. As a rule of thumb, the parallel efficiency should be greater than 0.5 to not waste resources."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8285bf07-0aed-4db7-9055-a313cbae131d",
   "metadata": {},
   "source": [
    "## Parallelization of workflows\n",
    "In a geoprocessing workflow, there are often multiple independent tasks that can be executed in parallel.\n",
    "The strategy how to divide the workflow into parallel tasks generally falls under either data-based or task-based parallelization.\n",
    "\n",
    "Task-based parallelism partitions tasks so that independent tasks can be completed simultaneously.\n",
    "\n",
    "With data-based parallelization, the spatial domain is partitioned for concurrent computations of individual spatial units \n",
    "and once processed, spatial units are mosaicked back into the initial spatial domain (if applicable).\n",
    "\n",
    "<img src=\"https://github.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023/blob/main/parallelization_strategies.png?raw=true\" alt=\"data-based and task-based diagrams\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc5906a5-83b8-490c-8c5e-a3194a8147ad",
   "metadata": {},
   "source": [
    "### Data-based parallelization\n",
    "Data-based parallelism involves spatial domain decomposition, a process of splitting data into smaller datasets that can be processed in parallel.\n",
    "As part of [GRASS GIS Python API](https://grass.osgeo.org/grass-stable/manuals/libpython/index.html), [GridModule](https://grass.osgeo.org/grass-stable/manuals/libpython/pygrass.modules.grid.html) decomposes input data into rectangular tiles, executes a given tool for each tile in parallel, and merges the results. Effectively, tiling is applied virtually (using computational region), determining the spatial extent of analysis for each parallel process. In some cases such as moving window analysis, tiles need to overlap to get correct results. Note that GridModule can be fairly inefficient due to the overhead from merging results back and is therefore best suited for computationally-itensive tasks such as interpolation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c35b1922-f689-43b1-8fa7-75138c4947bf",
   "metadata": {},
   "source": [
    "The following example shows IDW interpolation split into 4 tiles. In this case, specifying an overlap is needed to get correct results without edge artifacts. Here, the number and size of tiles is automatically derived from the number of cores, but can be specified. First we create sample data for this example by taking 10,000 random sample points from the DEM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1c51d47-5948-45ca-91f3-827cae68405e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region res=100\n",
    "r.random -z input=DEM npoints=10000 vector=samples seed=1 --q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "304566ca-5275-4c6c-aa8c-4fd098b51a93",
   "metadata": {},
   "source": [
    "We also set the resolution for our IDW computation and sampling to 100 using the g.region tool."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f55cc524-28b4-4ea1-8396-8645fc790c88",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">Each mapset is using a <strong>computational region</strong> that determines the extent and resolution of raster-based computations.\n",
    "It can be changed with <a href=\"https://grass.osgeo.org/grass-stable/manuals/g.region.html\">g.region</a> tool.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d54534c8-d654-4d39-98c7-39a71654ff14",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%python\n",
    "import time\n",
    "from grass.pygrass.modules.grid import GridModule\n",
    "import time\n",
    "start = time.time()\n",
    "grid = GridModule(\n",
    "    \"v.surf.idw\",\n",
    "    input=\"samples\",\n",
    "    output=\"idw\",\n",
    "    processes=4,\n",
    "    overlap=20,\n",
    "    quiet=True,\n",
    ")\n",
    "grid.run()\n",
    "print(f\"Elapsed time: {time.time() - start} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd2113e1-8572-418e-bc30-8e011e8d416f",
   "metadata": {},
   "source": [
    "The following is the same tool ran in serial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "488daca8-3763-44bb-be34-7aa741bef854",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "time v.surf.idw --q input=samples output=idw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b65dc4-088f-41cc-88d2-828f0c58e5c5",
   "metadata": {},
   "source": [
    "There are tools that already integrate tiling. For example, addon [r.mapcalc.tiled](https://grass.osgeo.org/grass-stable/manuals/addons/r.mapcalc.tiled.html) uses the tiling concept for raster algebra computation. Using this is better for more complex algebra expression and large input data, otherwise the parallel efficiency of this method can be low."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e5a1401-7c1a-4b0c-ba65-b5dbc8234c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region raster=DEM\n",
    "r.mapcalc.tiled expression=\"vertical_sobel = -DEM[-1, -1] - 2 * DEM[0, -1] - DEM[1, -1] + DEM[-1, 1] + 2 * DEM[0, 1] + DEM[1, 1]\" overlap=1 nprocs=4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6774328-e620-4656-b25b-c8e9b8fd2a5e",
   "metadata": {},
   "source": [
    "### Task-based parallelization\n",
    "With task-based parallelism, we identify independent tasks and execute them concurrently.\n",
    "Tasks are typically GRASS processing tools executed as separate processes. Processes, unlike threads, do not share memory. When tasks are limited by disk I/O, parallel processing may have large overhead.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1cf23a3-c647-40b7-b1e6-457862b3e588",
   "metadata": {},
   "source": [
    "#### Examples in Python\n",
    "There are multiple ways to execute tasks in parallel using Python, for example, there are libraries `multiprocessing` and `concurrent.futures`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4e0ef5d-8137-4724-b1a2-0b4f0ed7a452",
   "metadata": {},
   "source": [
    "In the following example viewsheds from different coordinates are computed in parallel using `multiprocessing.Pool` class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb1cfd6f-15b1-4896-9695-a5116167d5ae",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">Note that Python multiprocessing.Pool examples do not work in an interactive interpreter (such as Jupyter Notebook). That's why we will first write a Python script with %%writefile and then use %run magic to execute it.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c01f28-c156-4d29-94ee-caf7a5b7c3bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile example.py\n",
    "import os\n",
    "from multiprocessing import Pool\n",
    "import grass.script as gs\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    gs.run_command(\"r.viewshed\", input=\"DEM\", output=f\"viewshed_{cat}\", coordinates=(x, y), maxdistance=10000)\n",
    "    return f\"viewshed_{cat}\"\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    viewpoints = [(594231, 275545, 1),\n",
    "                  (659805, 259566, 2),\n",
    "                  (646109, 232901, 3),\n",
    "                  (602946, 203226, 4)]\n",
    "    with Pool(processes=2) as pool:\n",
    "        maps = pool.map(viewshed, viewpoints)\n",
    "    print(maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee83549a-0ff1-46e7-ad91-f89c0891eca6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region raster=DEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13c75b70-c157-434c-a14f-7ecf6dc8db26",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run example.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b966bcd6-3f03-4e38-aee8-cf491b12eaaa",
   "metadata": {},
   "source": [
    "#### Examples in Bash\n",
    "In the simplest case, tasks can be executed in parallel from a command line shell by running the geoprocessing tasks in the background (by appending `&`). Command `wait` forces to wait for the background processes to finish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63e997db-8cee-47d2-81bf-36f99e900558",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region raster=DEM res=100\n",
    "v.kernel --q input=samples output=kernel_10000 radius=10000 &\n",
    "v.kernel --q input=samples output=kernel_15000 radius=15000 &\n",
    "wait"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd37389-9b1a-4b40-8ebf-a7b7af891713",
   "metadata": {
    "tags": []
   },
   "source": [
    "Larger number of tasks can be scheduled to run in parallel by tools such as [GNU Parallel](https://www.gnu.org/software/parallel/) and xargs.\n",
    "In this simple example, we use a loop to write commands into a file and execute those commands in parallel, using 2 cores. \n",
    "Whenever a task is finished, a next one is picked from the queued tasks.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca4d707a-82e7-4a0c-886f-dd71b44aa186",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for R in 5000 10000 15000\n",
    "do\n",
    "    echo v.kernel --q input=samples output=kernel_${R} radius=${R} >> commands.sh\n",
    "done\n",
    "time parallel -j 2 < commands.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d6600de-765f-4c6c-927d-f75aa2120074",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">While waiting, let's look at the processors working using <code>htop</code>. Click on the blue button New Launcher, create a terminal, and run htop.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b3138f3-77b5-4e9e-a75d-936fc4b1fd88",
   "metadata": {},
   "source": [
    "See manual pages of GNU Parallel or xargs for more advanced uses. GNU Parallel can be configured to distribute jobs across multiple machines. In that case, use `--exec` interface described below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c5b1a00-77e6-4c26-b057-67b6823e6bfc",
   "metadata": {},
   "source": [
    "### Executing processes on distributed systems\n",
    "To enable parallelization on distributed systems, tasks or scripts need to be executed non-interactively using the `--exec` [interface](https://grass.osgeo.org/grass-stable/manuals/grass.html) interface.\n",
    "For that you need to specify project and mapset.\n",
    "\n",
    "\n",
    "For example, here is a simple call to list all available rasters in PERMANENT mapset:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e03e21ee-2a0b-4433-b992-fbd17d9eb44a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass gis_week_2023/PERMANENT --exec g.list type=raster mapset=PERMANENT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ef52b0d-9ea0-476e-a194-c3c84e5b78c4",
   "metadata": {},
   "source": [
    "Non-interactive tasks need to be run in separate mapsets. One of the previous examples that was running within GRASS session in a single mapset can be rewritten so that each task runs in a newly created mapset. Note that by default newly created mapsets use default computational region for that project (you can use `g.region -s` to modify it). For raster computations, you need to change the computational region for each new mapset if the default one is not desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "197a2701-93a6-4fc7-ba0a-c4c8d6f9e622",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for R in 5000 10000 15000\n",
    "do\n",
    "    # first create a new mapset with -c flag and set computational region\n",
    "    grass -c gis_week_2023/mapset_${R} --exec g.region raster=DEM res=100\n",
    "    # write the command executing v.kernel in the newly created mapset to a file\n",
    "    echo grass gis_week_2023/mapset_${R} --exec v.kernel --q input=samples@part_1 output=kernel_${R} radius=${R} >> exec_commands.sh\n",
    "done\n",
    "parallel -j 2 < exec_commands.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "814eb419-9050-4baf-99fb-03e5cc69f08d",
   "metadata": {},
   "source": [
    "In some cases, only a temporary mapset or project is needed, see [examples](https://grass.osgeo.org/grass-stable/manuals/grass.html#batch-jobs-with-the-exec-interface).\n",
    "Besides individual tools, the `--exec` interface can execute an entire script to enable more complex workflows."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}