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(\"opengeohub_2023/part_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38aff518-c001-4553-94dc-c1f695663796",
   "metadata": {},
   "source": [
    "The examples will run in *opengeohub_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=\"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": "24bebeaa-aff3-4471-a98a-77acee36766e",
   "metadata": {},
   "source": [
    "##### *Try it yourself*\n",
    "Using bash or Python, construct your own timed cell. You can use a sleep function or your own simple function (i.e. Hello World, 2+2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d8cf4f0-52cb-4628-bf85-962f45d5c8cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it here"
   ]
  },
  {
   "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 print min and max values of a raster with [r.info](https://grass.osgeo.org/grass-stable/manuals/r.info.html) using the following Bash and Python syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4baf2cac-44f5-49d6-9795-e13a8819b54a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.info map=DEM -r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a965c09-36bd-4d38-bdbb-0635e665d8d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.info\", map=\"DEM\", flags=\"r\")"
   ]
  },
  {
   "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": "7dccc8cf-c5d3-4a30-bef5-95051497ff3c",
   "metadata": {},
   "source": [
    "##### *Try it yourself*\n",
    "What's the speedup for the r.neighbors command above with 2 cores? The efficiency?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "759fab91-a0d8-4c6c-b401-b47c5bac100d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it here"
   ]
  },
  {
   "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=\"parallelization_strategies.png\" 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",
    "v.surf.idw --q input=samples output=idw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d68cc5-25a8-4179-8ca0-2dfe02e8e781",
   "metadata": {},
   "source": [
    "##### *Try it yourself*\n",
    "Time the serial computation above. What is the speedup?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb83164-5303-4938-9fb5-95678ebb605d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it here"
   ]
  },
  {
   "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 Location 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 opengeohub_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 GRASS location (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 opengeohub_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 opengeohub_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 location 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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8a5a5bf-53a2-4c47-bcaf-0ea12c585ca1",
   "metadata": {},
   "source": [
    "## Safe execution of parallel tasks (advanced)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dc8fcb0-2234-433c-ba55-978baad3c2a5",
   "metadata": {},
   "source": [
    "There are certain situations you need to avoid when running GRASS tools in parallel:\n",
    "\n",
    " * write output maps/files with identical names\n",
    " * modify computational region when running within the same mapset\n",
    " * modify MASK when running within the same mapset\n",
    " * modify vector attribute database within the same mapset\n",
    " * use [r.reclass](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) to reclassify from the same base map\n",
    "\n",
    "\n",
    "The following examples show how to deal with some of these situations when running parallel tasks in the same mapset is desired."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52e8628e-a307-4928-abde-c0ac981a0b5c",
   "metadata": {},
   "source": [
    "### Safely modifying computational region in a single mapset\n",
    "\n",
    "Sometimes modifying computational region in a script is needed. It is a good practice to not change the global computational region, which effectively modifies a file in a mapset,\n",
    "but only change the environment variable `GRASS_REGION`.\n",
    "Here, we modified the previous viewshed example to compute in parallel viewsheds with different extents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b9d056-0798-4d01-bdf2-5de6d83fa64e",
   "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",
    "    # copy current environment, modify and pass it to r.viewshed\n",
    "    env = os.environ.copy()\n",
    "    maxdistance = 10000\n",
    "    # create region based on the viewpoint and maxdistance\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(e=x + maxdistance, w=x - maxdistance, n=y + maxdistance, s=y - maxdistance, align=\"DEM\")\n",
    "    gs.run_command(\"r.viewshed\", input=\"DEM\", output=f\"viewshed_{cat}\", coordinates=(x, y), max_distance=maxdistance, env=env)\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": "ff2a9fe1-25ae-4db8-8f77-4b586dd51207",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run example.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9433010b-3ad7-44c1-a162-1cdda67cb655",
   "metadata": {},
   "source": [
    "### Safely modifying vectors with attributes in a single mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f9ade9d-847b-46ab-a776-cc2743e28d3a",
   "metadata": {},
   "source": [
    "By default vector maps share a single SQLite database file, however SQLite does not support concurrent write access. That poses a problem when modifying vectors with attributes in parallel. While this can be solved by running the computations in separate mapsets, it is also possible to change the default behavior to write attributes of each vector to the vector's individual SQLite file. This behavior can be activated after a new mapset is created with:\n",
    "\n",
    "```\n",
    " db.connect driver=sqlite database='$GISDBASE/$LOCATION_NAME/$MAPSET/vector/$MAP/sqlite.db'\n",
    "```\n",
    "\n",
    "Alternatively, a PostgreSQL or another database backend can be used for attributes to offload the parallel writing to the database system."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "439de6f9-494d-48ce-8260-ba20c1f87268",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81aaf8df-227e-4049-ab57-152a12c2b8e5",
   "metadata": {
    "tags": []
   },
   "source": [
    "Consider this (slightly simplified) example code taken from the [r.sun documentation](https://grass.osgeo.org/grass-stable/manuals/r.sun.html):\n",
    "```\n",
    "# slope + aspect\n",
    "r.slope.aspect elevation=DEM aspect=aspect slope=slope\n",
    "\n",
    "# calculate global irradiation for day 180\n",
    "r.sun elevation=DEM aspect=aspect slope=slope glob_rad=global_rad day=180 nprocs=2\n",
    "# result: output global (total) irradiation [Wh.m-2.day-1] for given day\n",
    "r.univar global_rad\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43259679-8543-471a-831b-8f94a78128f4",
   "metadata": {},
   "source": [
    "Use [r.sun](https://grass.osgeo.org/grass-stable/manuals/r.sun.html) to compute the solar radiation on the winter and summer soltices (day of year 355 and 172 respectively), the spring equinox (day 80) and your birthday (you can look it up [here](https://asd.gsfc.nasa.gov/Craig.Markwardt/doy2023.html)). Try to parallelize the workflow while considering available resources. What strategies for parallelizing have you used? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bbb23f5-ebf6-43c6-ad17-b8569dfa3249",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your work here (use Bash or Python)"
   ]
  }
 ],
 "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
}