{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "V5FdN92AUR_q"
   },
   "source": [
    "# Install GRASS GIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tW-I7tg_85X6"
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "DEBIAN_FRONTEND=noninteractive sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable && apt update && apt install grass subversion grass-dev parallel htop && apt remove libproj22\n",
    "git clone https://github.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023\n",
    "cd grass-workshop-gis-week-2023\n",
    "sh download_dataset.sh\n",
    "grass --tmp-mapset gis_week_2023 --exec g.extension r.futures\n",
    "grass --tmp-mapset gis_week_2023 --exec g.extension r.mapcalc.tiled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "UCdzD3IbTbFb"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir(\"grass-workshop-gis-week-2023\")"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47fb6e71-636b-48b4-b7ec-3ff20b5557a1",
   "metadata": {},
   "source": [
    "# Urban growth modeling in GRASS GIS: Parallel computing case study\n",
    "The purpose of this notebook is to demonstrate several parallel computing principles and how they are implemented in GRASS GIS.\n",
    "We use FUTURES urban growth model implemented as a GRASS GIS addon [r.futures](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.html).\n",
    "\n",
    "This notebook uses a [prepared dataset](https://doi.org/10.5281/zenodo.6577922) you downloaded at the beginning of the workshop. This dataset is a GRASS GIS project containing:\n",
    " * [NLCD 2001-2019](https://www.mrlc.gov/) (National Land Cover Database; land use and impervious surface descriptor)\n",
    " * [US county boundaries](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html)\n",
    " * [US-PAD protected areas](https://www.usgs.gov/programs/gap-analysis-project/science/pad-us-data-overview)\n",
    " * [USGS DEM](https://www.usgs.gov/3d-elevation-program/about-3dep-products-services)\n",
    " \n",
    "Population files were downloaded from [Zenodo](https://doi.org/10.5281/zenodo.6577903).\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "befd1cb4-24ee-467b-b969-0412d9af1f14",
   "metadata": {},
   "source": [
    "<div  class=\"alert alert-info\">This notebook combines Python 3 and Bash cells. By default a code cell is in Python.\n",
    "We use IPython <a href=\"https://ipython.readthedocs.io/en/stable/interactive/magics.html#cell-magics\">cell magic</a> including `%%bash`, `%%time`, `%%timeit` and `%%writefile`.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d82edc-889d-4a41-aaa8-3e61ec03c489",
   "metadata": {},
   "source": [
    "## Setting up\n",
    "Import Python packages and initialize GRASS GIS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a90175-1e87-4dda-a82f-47da95c5cf9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import sys\n",
    "import pathlib\n",
    "import json\n",
    "import pandas as pd\n",
    "from IPython.display import Image\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_2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63e1921b-0635-4178-b489-692287aeadc8",
   "metadata": {},
   "source": [
    "## Data preprocessing\n",
    "Before we start processing, let's take a look at what data we have and our study area. List dataset layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea67a34b-141d-4b8d-9333-69e29ea07664",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.list type=raster,vector -p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a976ec2-a48d-49f8-9110-412c681a0fd0",
   "metadata": {},
   "source": [
    "Display counties using [grass.jupyter.InteractiveMap](https://grass.osgeo.org/grass83/manuals/libpython/grass.jupyter.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d896545e-649a-4cb3-a463-3253b66befeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.InteractiveMap()\n",
    "m.add_vector(name=\"counties\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f33772fb-4751-4627-bc4a-441846226e48",
   "metadata": {},
   "source": [
    "Set computational region to match counties' extent and aligns with the NLCD raster. Computational region defines the extent and resolution of all raster computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7545a6d0-4594-46ec-ba11-e015fb1f9fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", vector=\"counties\", align=\"nlcd_2019\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbb6bc6c-a25d-4d0f-a8d4-3d82b1f80c83",
   "metadata": {},
   "source": [
    "Now, we will prepocess the data to derive spatial predictors used by the urban growth model in this case study. This diagram shows the workflow and highlights the parts that we will parallelize.\n",
    "\n",
    "![FUTURES data preparation](https://raw.githubusercontent.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023/aaf6ec7c2b249db68cb79b6d1042247a80381c44/FUTURES_case_study_data_prep_opengeohub.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3d84f42-9e42-4667-9587-1eb41e6605a3",
   "metadata": {},
   "source": [
    "### DEM to slope\n",
    "Compute slope with [r.slope.aspect](https://grass.osgeo.org/grass-stable/manuals/r.slope.aspect.html) which uses OpenMP for parallelization. For fun, measure the time difference between using 1 or more cores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "009f8b8d-e2d3-44e8-97af-b2bace75da9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n2 -r3\n",
    "gs.run_command(\n",
    "    \"r.slope.aspect\", elevation=\"DEM\", slope=\"slope\", flags=\"e\", nprocs=1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1240eafa-d06f-4f98-b3df-c1b5848b4699",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n2 -r3\n",
    "gs.run_command(\n",
    "    \"r.slope.aspect\", elevation=\"DEM\", slope=\"slope\", flags=\"e\", nprocs=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c602a93e-af65-4978-8179-02e80f2dd6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"slope\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.d_legend(raster=\"slope\", range=[0, 10])\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb0fc69-b2bc-4962-bcf9-2236835d0507",
   "metadata": {},
   "source": [
    "### Derive predictors and other layers from NLCD data\n",
    "Most of our predictors we will derive from NLCD data (land cover type and impervious descriptor products). See [NLCD legend](https://www.mrlc.gov/data/legends/national-land-cover-database-class-legend-and-description) for classification. With [r.reclass](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) we will create separate layers for water, wetland, forest, roads, and developed land.\n",
    "Note that those rasters are virtual (they behave the same way, but are only pointing to the original NLCD raster),\n",
    "so reclassification is very fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95497c40-7131-4b37-be96-8413d7316426",
   "metadata": {},
   "outputs": [],
   "source": [
    "NLCD_years = [2001, 2004, 2006, 2008, 2011, 2013, 2016, 2019]\n",
    "NLCD_start_end_years = [2001, 2019]\n",
    "# water (1 or no data, class 11)\n",
    "gs.write_command(\n",
    "    \"r.reclass\", input=\"nlcd_2019\", output=\"water\", rules=\"-\", stdin=\"11 = 1\"\n",
    ")\n",
    "# binary wetlands (classes 90 and 95)\n",
    "gs.write_command(\n",
    "    \"r.reclass\",\n",
    "    input=\"nlcd_2019\",\n",
    "    output=\"wetlands\",\n",
    "    rules=\"-\",\n",
    "    stdin=\"90 95 = 1 \\n * = 0\",\n",
    ")\n",
    "# developed land (classes 21 to 24)\n",
    "for year in NLCD_years:\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_{year}\",\n",
    "        output=f\"urban_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"21 22 23 24 = 1\\n* = 0\",\n",
    "    )\n",
    "for year in NLCD_start_end_years:\n",
    "    # forest classes 41 to 43\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_{year}\",\n",
    "        output=f\"forest_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"41 42 43 = 1\",\n",
    "    )\n",
    "    # roads classes 20-23 in the descriptor layer\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_descriptor_{year}\",\n",
    "        output=f\"roads_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"20 21 22 23 = 1\",\n",
    "    )\n",
    "    # urban without roads (classes 24-26 in the descriptor layer)\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_descriptor_{year}\",\n",
    "        output=f\"urban_no_roads_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"24 25 26 = 1\\n* = 0\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e0c984f-0869-4f7c-b274-ae4bcce8f478",
   "metadata": {},
   "source": [
    "Next, we will compute **distance to water, forest, and roads** (from different years). One of the simplest way to compute these independent tasks in parallel is to run them in the background using Bash by appending &.\n",
    "Command `wait` forces to wait for the background processes to finish. Adding time in front of wait will report the elapsed time.\n",
    "Once the distance is computed, we will use raster algebra to transform it logarithmically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e39a059-7334-4f36-8b08-0c9e5e74b5b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.grow.distance input=water distance=dist_to_water -m --q &\n",
    "r.grow.distance input=forest_2001 distance=dist_to_forest_2001 -m --q &\n",
    "r.grow.distance input=forest_2019 distance=dist_to_forest_2019 -m --q &\n",
    "r.grow.distance input=roads_2001 distance=dist_to_roads_2001 -m --q &\n",
    "r.grow.distance input=roads_2019 distance=dist_to_roads_2019 -m --q &\n",
    "time wait\n",
    "r.mapcalc \"log_dist_to_water = log(dist_to_water + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_forest_2001 = log(dist_to_forest_2001 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_forest_2019 = log(dist_to_forest_2019 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_roads_2001 = log(dist_to_roads_2001 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_roads_2019 = log(dist_to_roads_2019 + 1)\" --q &\n",
    "time wait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aa010b9-3d0e-42b1-a0eb-84eb96c68e20",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"log_dist_to_water\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4089ca-0670-40f0-a890-42af4dc2d1af",
   "metadata": {},
   "source": [
    "As another predictor, we compute **wetland density** (percentage of wetland in 1 km squared circular neighborhood). Module [r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html) for moving window analysis is internally parallelized using OpenMP, so we can use `nprocs` option:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46f73dca-2556-47ec-8c1b-6f339584f825",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\n",
    "    \"r.neighbors\",\n",
    "    input=\"wetlands\",\n",
    "    output=\"wetland_density\",\n",
    "    size=37,\n",
    "    method=\"average\",\n",
    "    flags=\"c\",\n",
    "    nprocs=4,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed777a2f-47b9-4f0d-b5cc-b9f5820e5182",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"wetland_density\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77055576-de4f-4d11-9d07-14d6f4ba4837",
   "metadata": {},
   "source": [
    "FUTURES uses a special predictor called development pressure, which can be computed with [r.futures.devpressure](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.devpressure.html), which is internally parallelized.\n",
    "Since we need to compute it for 2 years, if we have enough cores, we can use a hybrid approach which runs both commands as background process and each of them runs in parallel.\n",
    "To do that we split the number of available processes so that each r.futures.devpressure process gets half of the available processes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e684d44-4b27-44f3-b161-4c44ee15fb7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.futures.devpressure input=urban_no_roads_2001 output=devpressure_2001 size=15 gamma=0.5 nprocs=2 scaling_factor=0.1 &\n",
    "r.futures.devpressure input=urban_no_roads_2019 output=devpressure_2019 size=15 gamma=0.5 nprocs=2 scaling_factor=0.1 &\n",
    "time wait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f40180a-0cb4-466b-a781-9db89fe3e334",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"devpressure_2001\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c3c863-107d-4a0b-8d6c-1ec4d83ade26",
   "metadata": {},
   "source": [
    "### Mask\n",
    "Compute mask to avoid urban growth simulation in water, protected areas, and outside the study area.\n",
    "Here we demonstrate data-based parallelization using  [r.mapcalc.tiled](https://grass.osgeo.org/grass-stable/manuals/addons/r.mapcalc.tiled.html). Note that with our relatively small dataset, this approach is going to actually slow down the computation in comparison to simple r.mapcalc call due to overhead from patching the results together. However, for larger dataset, we would be able to see speedup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "944a7d7c-ee58-4223-9b68-ea3a174c0c85",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\n",
    "    \"r.mapcalc.tiled\",\n",
    "    expression=\"masking = if((isnull(protected) &&  isnull(water) && nlcd_2019), 1, null())\",\n",
    "    nprocs=4,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59384732-5c75-4656-ac1a-d7a1c86ada75",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.mapcalc(\"masking = if((isnull(protected) &&  isnull(water) && nlcd_2019), 1, null())\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c3c70a1-8c3d-4ab0-bf31-46d16a5c8a2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"masking\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdeac7a4-b38c-4f99-af56-0c14195a06a7",
   "metadata": {},
   "source": [
    "## FUTURES simulation\n",
    "FUTURES simulation has several components to compute how much land is going to be developed, where, and with what size of patches:\n",
    "\n",
    "![FUTURES diagram](https://raw.githubusercontent.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023/aaf6ec7c2b249db68cb79b6d1042247a80381c44/FUTURES.svg)\n",
    "\n",
    "To keep this workshop on schedule, in this part we will demonstrate the computation of **land demand** using [r.futures.demand](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.demand.html), then we skip a few steps and proceed with the patch growing simulation. If interested, you can find the skipped steps in `FUTURES_potential_and_calibration.ipynb`.\n",
    "\n",
    "\n",
    "\n",
    "### Land demand\n",
    "Here we compute how much land will be developed in each step of the simulation.\n",
    "Logarithmic curves are fit to per-capita land consumption data derived from NLCD time series and observed population for each county. \n",
    "\n",
    "Here we wanted to demonstrate a parallelization approach that splits spatial data into spatial units that are then processed in parallel, and that can be run in a distributed way (on an HPC). Technically, each spatial unit will be processed in a separate GRASS mapset. Distributing of the tasks is done here with GNU Parallel, but that would depend on the specific HPC setup.\n",
    "While US states would be more appropriate spatial units for this approach (given their size), in this small case study we will apply this approach to counties.\n",
    "\n",
    "![FUTURES simulation_with highlighted parts of the workflow](https://raw.githubusercontent.com/ncsu-geoforall-lab/grass-workshop-gis-week-2023/aaf6ec7c2b249db68cb79b6d1042247a80381c44/FUTURES_case_study_simulation_nosketch_opengeohub.svg)\n",
    "\n",
    "First, we split and rasterize counties for further parallelization steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77dcaed-ea27-4921-927a-66741edfb3a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of all county ids\n",
    "county_ids = gs.read_command(\"v.db.select\", map=\"counties\", column=\"FIPS\", format=\"csv\", flags=\"c\").splitlines()\n",
    "print(county_ids)\n",
    "# use temporary region\n",
    "gs.use_temp_region()\n",
    "for fips in county_ids:\n",
    "    gs.run_command(\n",
    "        \"v.extract\",\n",
    "        input=\"counties\",\n",
    "        where=f\"FIPS == '{fips}'\",\n",
    "        output=f\"county_{fips}\",\n",
    "    )\n",
    "    gs.run_command(\"g.region\", vector=f\"county_{fips}\", align=\"nlcd_2019\")\n",
    "    gs.run_command(\n",
    "        \"v.to.rast\",\n",
    "        input=f\"county_{fips}\",\n",
    "        output=f\"county_{fips}\",\n",
    "        use=\"attr\",\n",
    "        attribute_column=\"FIPS\",\n",
    "    )\n",
    "gs.del_temp_region()\n",
    "gs.run_command(\"g.remove\", pattern=\"county_*\", type=\"vector\", flags=\"f\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e498282-b24f-48b7-8549-31dff94f1f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map(use_region=True)\n",
    "m.d_rast(map=\"county_37183\")\n",
    "m.d_vect(map=\"counties\", attribute_column=\"FIPS\", fill_color=\"none\", label_color=\"black\", xref=\"center\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccbed82f-d7ae-4db5-ad83-155d4f4aa05a",
   "metadata": {},
   "source": [
    "Then, we create a Python script that takes a county id as an input parameter,\n",
    "sets the computational region to the county extent, excludes roads from the computation, and runs [r.futures.demand](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.demand.html),\n",
    "creating an output CSV needed as input for r.futures.simulation and a plot for visual check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d4b82e2-c472-477d-9290-d9e13680204b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile demand_for_county.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "\n",
    "# input parameter: state FIPS code\n",
    "fips = sys.argv[1]\n",
    "\n",
    "gs.run_command(\"g.mapsets\", mapset=\"part_2\", operation=\"add\")\n",
    "gs.run_command(\"g.region\", raster=f\"county_{fips}\")\n",
    "gs.mapcalc(\"MASK = if (isnull(roads_2019), 1, null())\")\n",
    "gs.run_command(\n",
    "    \"r.futures.demand\",\n",
    "    subregions=f\"county_{fips}\",\n",
    "    development=[\n",
    "        f\"urban_{year}\" for year in [2001, 2004, 2006, 2008, 2011, 2013, 2016, 2019]\n",
    "    ],\n",
    "    observed_population=\"observed_population_SE_counties_2001-2019.csv\",\n",
    "    projected_population=\"projected_population_SE_counties_2020-2100_SSP2.csv\",\n",
    "    simulation_times=list(range(2019, 2051)),\n",
    "    method=\"logarithmic\",\n",
    "    demand=f\"demand_{fips}.csv\",\n",
    "    plot=f\"demand_{fips}.png\",\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4d24955-09f9-4085-9caa-c0b72c13bedf",
   "metadata": {},
   "source": [
    "For each county we will use Bash scripting to generate grass `--exec` command calling the script within a temporary mapset and append each line to a file `demand_jobs.sh`.\n",
    "Then, we will run these commands in parallel with GNU Parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21b9404e-8fa7-4894-8c0f-3e752f2740d0",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">See <a src=\"https://www.freecodecamp.org/news/bash-scripting-tutorial-linux-shell-script-and-command-line-for-beginners/\">Bash scripting basics</a> \n",
    "or other sources to understand bash commands, for loops, or command substitution.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c436b2d-6ad3-43fe-be52-07139435d5b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# remove the job file (for repeated executions of this cell)\n",
    "rm -f demand_jobs.sh\n",
    "# for each of the county ids\n",
    "for S in $(v.db.select map=\"counties\" column=\"FIPS\" format=\"csv\" -c)\n",
    "do\n",
    "    # write \"grass ...\" command into demand_jobs.sh\n",
    "    # we will use temporary mapset, no need to store it afterwards\n",
    "    echo grass --tmp-mapset gis_week_2023 --exec python demand_for_county.py ${S} >> demand_jobs.sh\n",
    "done\n",
    "# print\n",
    "cat demand_jobs.sh\n",
    "# specify number of parallel jobs, reroute error output\n",
    "parallel -j 2 < demand_jobs.sh 2> log.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a82a1f1-7c59-43c0-bc2b-f6885dde35a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"demand_37183.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae2df84-1d7e-4b51-9802-8ddb6c43a00d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">At this point, we skip several steps of FUTURES modeling workflow to keep this short and relevant. We will use precomputed files included in the repository. If you are interested in the skipped workflow, see FUTURES_potential_and_calibration.ipynb.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8c6b04-c0f7-4b41-9242-06222f6777e1",
   "metadata": {},
   "source": [
    "### Patch Growing Algorithm\n",
    "With [r.futures.simulation](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.simulation.html) we will run the simulation from 2019 until 2050. r.futures.simulation is stochastic, so with different random seeds, we will get somewhat different result. We will use the same parallelization approach as with r.futures.demand.\n",
    "\n",
    "First, create a Python script with county id and random seed as an input. The script sets computational region to match the county, turns on the mask, and runs the simulation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f93a8e-e266-488d-8e68-54e8c5928e7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile simulation_for_county.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "\n",
    "fips, seed = sys.argv[1:3]\n",
    "\n",
    "# g.mapsets for accessing maps from part_2 mapset\n",
    "gs.run_command(\"g.mapsets\", mapset=\"part_2\", operation=\"add\")\n",
    "gs.run_command(\"g.region\", raster=f\"county_{fips}\")\n",
    "gs.run_command(\"r.mask\", raster=\"masking\")\n",
    "gs.run_command(\n",
    "    \"r.futures.simulation\",\n",
    "    developed=\"urban_2019\",\n",
    "    development_pressure=\"devpressure_2019\",\n",
    "    compactness_mean=0.4,\n",
    "    compactness_range=0.1,\n",
    "    discount_factor=0.5,\n",
    "    predictors=[\n",
    "        \"log_dist_to_forest_2019\",\n",
    "        \"log_dist_to_roads_2019\",\n",
    "        \"log_dist_to_water\",\n",
    "        \"slope\",\n",
    "        \"wetland_density\",\n",
    "    ],\n",
    "    n_dev_neighbourhood=15,\n",
    "    devpot_params=\"best_model.csv\",\n",
    "    num_neighbors=4,\n",
    "    seed_search=\"probability\",\n",
    "    development_pressure_approach=\"gravity\",\n",
    "    gamma=0.5,\n",
    "    scaling_factor=0.1,\n",
    "    subregions=f\"county_{fips}\",\n",
    "    demand=f\"demand_{fips}.csv\",\n",
    "    num_steps=31,\n",
    "    output=f\"out_county_{fips}_seed_{seed}\",\n",
    "    patch_sizes=f\"patch_sizes/patch_sizes_{fips}.csv\",\n",
    "    incentive_power=2,\n",
    "    random_seed=seed,\n",
    ")\n",
    "gs.run_command(\"r.mask\", flags=\"r\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "112abbdd-2dc8-48f9-9c14-445e0c43f90f",
   "metadata": {},
   "source": [
    "Similarly, we create a list of commands executing this Python script with a given county and seeds (from 1 to 10).\n",
    "We will create a new mapset for each command. There will be *number of counties* * *number of random seeds* commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0e3b51e-d790-4135-ad6f-5e1a56354fcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "rm -f pga_jobs.sh\n",
    "for SEED in {1..10}\n",
    "do\n",
    "    for COUNTY in $(v.db.select map=\"counties\" column=\"FIPS\" format=\"csv\" -c)\n",
    "    do\n",
    "        # to start fresh delete the mapset in case it exists from repeated runs.\n",
    "        rm -rf gis_week_2023/pga_${COUNTY}_${SEED}\n",
    "        echo grass -c gis_week_2023/pga_${COUNTY}_${SEED} --exec python simulation_for_county.py ${COUNTY} ${SEED} >> pga_jobs.sh\n",
    "    done\n",
    "done\n",
    "cat pga_jobs.sh\n",
    "parallel -j 4 < pga_jobs.sh 2> log.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "332aacb8-a6b2-4aca-bfa0-b86ec77f3cb0",
   "metadata": {},
   "source": [
    "Afterwards, we patch the results from all the counties together. Tool [r.patch](https://grass.osgeo.org/grass-stable/manuals/r.patch.html) is internally parallelized, so we can use that extra speed up if we have available cores. For each seed we will get the list of layers and patch them. This simple loop can be easily parallelized in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88c2370-8a00-471d-8f70-9478debf07fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile patch.py\n",
    "\n",
    "import tqdm\n",
    "from multiprocessing import Pool\n",
    "import grass.script as gs\n",
    "\n",
    "\n",
    "def patch(seed):\n",
    "    maps = gs.read_command(\"g.list\", type=\"raster\", pattern=f\"out_county_*_seed_{seed}\",\n",
    "                           mapset=\"*\", flags=\"m\", separator=\"comma\").strip()\n",
    "    gs.run_command(\"r.patch\", input=maps, output=f\"out_seed_{seed}\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    with Pool(processes=2) as pool:\n",
    "        # a simple way to run it:\n",
    "        # pool.map(patch, list(range(1, 11)))\n",
    "        # more complex, but gives you progress bar:\n",
    "        r = list(tqdm.tqdm(pool.imap(patch, list(range(1, 11))), total=10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72f44dd8-460e-41bb-a998-96e6c4189a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run patch.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d890adfd-96ef-4550-8893-2cccb8b77acb",
   "metadata": {},
   "source": [
    "Set the color ramp of all merged simulation outputs for visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52c6a4fa-8060-4812-ba96-a0db01d07909",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"r.colors\",\n",
    "    map=\"out_seed_1\",\n",
    "    raster=\"out_county_37183_seed_1@pga_37183_1\",\n",
    ")\n",
    "m = gj.InteractiveMap()\n",
    "m.add_raster(name=\"out_seed_1\", opacity=0.8)\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf7dbfe-9e6f-406c-8c95-e7b2fc8a1395",
   "metadata": {},
   "source": [
    "## FUTURES results postprocessing\n",
    "In this last section we will derive some useful outputs from the stochastic simulation runs.\n",
    "\n",
    "### Future Development Probability\n",
    "By aggregating the stochastic runs, we can compute the projected development probability. First we reclassify output\n",
    "to binary developed/undeveloped results. Then we run [r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) in parallel to compute how many times a cell was developed\n",
    "and then divide that by number of runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cd1212d-9142-4982-a993-cafdd7aae7a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile reclass.txt\n",
    "-1 0 = 0\n",
    "1 thru 100 = 1\n",
    "* = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c44a6c-6389-4a35-b8c0-2db5a7b9bfe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "for seed in range(1, 11):\n",
    "    gs.run_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"out_seed_{seed}\",\n",
    "        output=f\"out_seed_{seed}_binary\",\n",
    "        rules=\"reclass.txt\",\n",
    "    )\n",
    "gs.run_command(\n",
    "    \"r.series\",\n",
    "    input=[f\"out_seed_{seed}_binary\" for seed in range(1, 11)],\n",
    "    output=\"probability\",\n",
    "    method=\"sum\",\n",
    "    weights=[0.1] * 10,\n",
    "    nprocs=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8364fc-632b-430a-a429-738b9ce62764",
   "metadata": {},
   "outputs": [],
   "source": [
    "# zoom in to see more details\n",
    "gs.run_command(\n",
    "    \"g.region\",\n",
    "    raster=\"county_37183\",\n",
    "    save=\"zoomin\",\n",
    ")\n",
    "gs.run_command(\"r.colors\", map=\"urban_2019\", color=\"grey\")\n",
    "m = gj.Map(saved_region=\"zoomin\", width=700)\n",
    "m.d_background(color=\"grey\")\n",
    "m.d_rast(map=\"probability\", values=\"0.2-1\")\n",
    "m.d_rast(map=\"urban_2019\", values=1)\n",
    "m.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}