workshop_part_3_parallelization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 5: Intro to Parallelization\n",
    "\n",
    "This part will briefly cover how to run GRASS computations in parallel. It requires computing the first two parts of the workshop. \n",
    "\n",
    "First, create a new mapset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/dix_park/parallelization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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 the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"dix_park\", \"parallelization\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tool-level parallelization\n",
    "There are several [internally parallelized tools](https://grass.osgeo.org/grass-stable/manuals/keywords.html#parallel), either using OpenMP or Python multiprocessing library. We can use `nprocs` option to set the number of cores to be used for processing.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set computational region to match dsm raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute moving window analysis and measure time first with one core, then with 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\"r.neighbors\", input=\"dsm\", output=\"dsm_smoothed\", method=\"average\", size=25, nprocs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\"r.neighbors\", input=\"dsm\", output=\"dsm_smoothed\", method=\"average\", size=25, nprocs=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize original and smoothed raster (turn layers on and off):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors_map = gj.InteractiveMap()\n",
    "neighbors_map.add_raster(\"dsm\")\n",
    "neighbors_map.add_raster(\"dsm_smoothed\")\n",
    "neighbors_map.add_layer_control(position=\"bottomright\")\n",
    "neighbors_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GridModule for tiling\n",
    "Some compute-intensive tasks can benefit from spatially splitting the task into tiles, and then running the task in parallel. [GridModule](https://grass.osgeo.org/grass-stable/manuals/libpython/pygrass.modules.grid.html) can automate this splitting-computing-merging procedure and execute the computation in parallel.\n",
    "\n",
    "In this example, we will interpolate an elevation surface from vector points using IDW interpolation. First, generate points randomly by sampling the digital elevation model. We adjusted the resolution to run the examples faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", res=2)\n",
    "gs.run_command(\"r.random\", flags=\"z\", input=\"ground\", npoints=5000, vector=\"samples\", seed=1, quiet=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measure the time without using GridModule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\"v.surf.idw\", input=\"samples\", output=\"idw\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now with GridModule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile interpolation.py\n",
    "from grass.pygrass.modules.grid import GridModule\n",
    "\n",
    "\n",
    "grid = GridModule(\n",
    "    \"v.surf.idw\",\n",
    "    input=\"samples\",\n",
    "    output=\"idw\",\n",
    "    processes=3,\n",
    "    overlap=20,\n",
    "    quiet=True,\n",
    ")\n",
    "grid.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "!python interpolation.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors_map = gj.Map()\n",
    "neighbors_map.d_rast(map=\"idw\")\n",
    "neighbors_map.d_vect(map=\"samples\", size=1, color=\"black\")\n",
    "neighbors_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running multiple independent computations\n",
    "In this example, our goal is to compute multiple viewsheds and render them into a PNG file.\n",
    "Since these are independent computations, we can run them in parallel.\n",
    "The first part implements this task in Python using _multiprocessing_ library\n",
    "and the second part will run each computation using `grass --exec` interface in separate mapsets that allows us to potentially distribute the computation across multiple nodes on an HPC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First compute a shaded relief raster for visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")\n",
    "gs.run_command(\"r.relief\", input=\"dsm\", output=\"relief\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the viewshed points form part 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewpoints = gs.read_command('v.out.ascii', input='viewpoints@viewshed',\n",
    "                             separator='comma', layer=2, where=\"height < 2\").strip().splitlines()\n",
    "viewpoints = [p.split(\",\") for p in viewpoints]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will run viewshed computation from part 2 using `multiprocessing.Pool`. We define a function that computes the viewshed and returns the name of the output or None in case of error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"dsm\", output=name,\n",
    "                       coordinates=(x, y), max_distance=300, flags=\"b\")\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "proc = 2\n",
    "with Pool(processes=proc) as pool:\n",
    "    maps = pool.map(viewshed, viewpoints)\n",
    "print(maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One trick to speedup viewshed computation is to limit the computation only\n",
    "to the actual area given by the maxdistance option. To do that we will locally modify the computational region\n",
    "and pass the environment to the module directly. The current computational region won't be affected. \n",
    "Additionally, this shows how to include a simple progress bar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    max_distance = 300\n",
    "    # set GRASS_REGION variable using region_env function\n",
    "    os.environ[\"GRASS_REGION\"] = gs.region_env(align=\"dsm\",\n",
    "                                               e=x + max_distance,\n",
    "                                               w=x - max_distance,\n",
    "                                               n=y + max_distance,\n",
    "                                               s=y - max_distance)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"dsm\", output=name, flags=\"b\",\n",
    "                      coordinates=(x, y), max_distance=max_distance)\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "proc = 2\n",
    "with Pool(processes=proc) as pool:\n",
    "    maps = list(tqdm(pool.imap(viewshed, viewpoints), total=len(viewpoints)))\n",
    "print(maps)\n",
    "print(f\"Viewshed num cells: {gs.raster_info(maps[0])['cells']}\")\n",
    "print(f\"DSM num cells: {gs.raster_info('dsm')['cells']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we extend this script to include rendering to file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    max_distance = 300\n",
    "    # set GRASS_REGION variable using region_env function\n",
    "    os.environ[\"GRASS_REGION\"] = gs.region_env(align=\"dsm\",\n",
    "                                               e=x + max_distance,\n",
    "                                               w=x - max_distance,\n",
    "                                               n=y + max_distance,\n",
    "                                               s=y - max_distance)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"dsm\", output=name, flags=\"b\",\n",
    "                      coordinates=(x, y), max_distance=max_distance)\n",
    "        # create visualization\n",
    "        viewshed_map = gj.Map(use_region=True)\n",
    "        viewshed_map.d_rast(map=\"relief\")\n",
    "        viewshed_map.d_rast(map=f\"viewshed_{cat}\", values=1)\n",
    "        viewshed_map.d_vect(map=\"viewpoints@viewshed\", layer=2, cat=cat, size=15, icon=\"basic/pin\")\n",
    "        viewshed_map.save(f\"viewshed_{cat}.png\")\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "nprocs = 2\n",
    "with Pool(processes=nprocs) as pool:\n",
    "    maps = list(tqdm(pool.imap(viewshed, viewpoints), total=len(viewpoints)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at one of the computed and rendered viewsheds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(\"viewshed_15.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this way, we can't distribute the computation across multiple nodes (hundreds of cores).\n",
    "We will do the same thing differently, using `grass --exec` [interface](https://grass.osgeo.org/grass-stable/manuals/grass.html), running each task in a separate mapset. This way, the tasks could be distributed across multiple nodes.\n",
    "\n",
    "`--exec` interface allows GRASS tools and user scripts to be executed in a GRASS non-interactive session. For example, here is a simple call to list all available vectors in PERMANENT mapset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass ~/grassdata/dix_park/PERMANENT --exec g.list type=vector mapset=viewshed -t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will create a Python script `myscript.py` computing and rendering viewsheds similarly as in the previous example. The script requires 3 parameters (x and y coordinate, and category). Note that we can set computational region in a standard way, because each script will run in separate mapset, so the different regions won't interfere with each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile myscript.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "\n",
    "def main(x, y, cat):\n",
    "    max_distance = 300\n",
    "    x, y = float(x), float(y)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    gs.run_command(\"g.region\", align=\"dsm\", e=x + max_distance,\n",
    "                   w=x - max_distance, n=y + max_distance, s=y - max_distance)\n",
    "    gs.run_command(\"r.viewshed\", input=\"dsm\", output=name, coordinates=(x, y),\n",
    "                   observer_elevation=3, max_distance=max_distance, flags=\"b\")\n",
    "    # create visualization\n",
    "    viewshed_map = gj.Map(use_region=True)\n",
    "    viewshed_map.d_rast(map=\"relief@parallelization\")\n",
    "    viewshed_map.d_rast(map=f\"viewshed_{cat}\", values=1)\n",
    "    viewshed_map.d_vect(map=\"viewpoints@viewshed\", layer=2, cat=cat, size=15, icon=\"basic/pin\")\n",
    "    viewshed_map.save(f\"viewshed_{cat}.png\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    args = sys.argv[1:]\n",
    "    main(*args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will generate a file `jobs.sh` with one command per line. We will run each task in a temporary mapset so all computed data will be deleted afterwards. That is fine for our example where we need only the final PNG files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"jobs.sh\", \"w\") as f:\n",
    "    for viewpoint in viewpoints:\n",
    "        f.write(f\"grass --tmp-mapset ~/grassdata/dix_park --exec python myscript.py {viewpoint[0]} {viewpoint[1]} {viewpoint[2]}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the content of the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat jobs.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To execute these commands in parallel, we can use e.g. [GNU Parallel](https://www.gnu.org/software/parallel/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "parallel -j 2 < jobs.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check one of the resulting PNG files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"viewshed_15.png\")"
   ]
  }
 ],
 "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": 4
}