notebook_3_viewshed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing and Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Existing Tools\n",
    "\n",
    "To call a GRASS tool from Python, we will use _run_command_ from the _grass.script_ package (typically imported as _gs_):\n",
    "\n",
    "```python\n",
    "gs.run_command(\n",
    "    \"r.viewshed\",\n",
    "    input=elevation,\n",
    "    output=\"viewshed\",\n",
    "    coordinates=(1200, 3000)\n",
    ")\n",
    "```\n",
    "\n",
    "For tools with textual outputs, _read_command_ can be used to capture that textual output in Python (as a string):\n",
    "\n",
    "```python\n",
    "table_data = gs.read_command(\n",
    "    \"r.univar\",\n",
    "    map=\"elevation\",\n",
    "    flags=\"t\",\n",
    "    separator=\"comma\",\n",
    ")\n",
    "```\n",
    "\n",
    "There are other ways of calling GRASS tools from Python, but using functions from the _run_command_ family is the most typical way for GRASS tools.\n",
    "\n",
    "## Processing Parameters\n",
    "\n",
    "The values parsed from the command line are stored in a dictionary returned by the _parse_ function. They can be accessed using `dictionary[\"name\"]` syntax where _name_ is the name of the parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile viewscape.py\n",
    "#!/usr/bin/env python\n",
    "\n",
    "# %module\n",
    "# % description: Compute viewshed and compute statistics about visible parts of sample layers\n",
    "# % keyword: raster\n",
    "# % keyword: statistics\n",
    "# % keyword: viewshed\n",
    "# %end\n",
    "# %option G_OPT_R_ELEV\n",
    "# % description: Name of input elevation raster map\n",
    "# %end\n",
    "# %option G_OPT_M_COORDS\n",
    "# %end\n",
    "# %option G_OPT_R_INPUTS\n",
    "# %end\n",
    "# %option G_OPT_F_OUTPUT\n",
    "# %end\n",
    "\n",
    "import subprocess\n",
    "import sys\n",
    "import csv\n",
    "import io\n",
    "\n",
    "import grass.script as gs\n",
    "\n",
    "\n",
    "def viewshed(\n",
    "    elevation,\n",
    "    coords,\n",
    "    output,\n",
    "    sample_continuous,\n",
    "):\n",
    "    name = \"viewshed\"\n",
    "    gs.run_command(\n",
    "        \"r.viewshed\",\n",
    "        input=elevation,\n",
    "        output=name,\n",
    "        coordinates=coords,\n",
    "        flags=\"cb\",\n",
    "    )\n",
    "    gs.run_command(\"r.null\", map=name, setnull=0)\n",
    "    results = {}\n",
    "    results[\"x\"] = coords[0]\n",
    "    results[\"y\"] = coords[1]\n",
    "    for raster in sample_continuous:\n",
    "        table_data = gs.read_command(\n",
    "            \"r.univar\",\n",
    "            map=raster,\n",
    "            zones=name,\n",
    "            quiet=True,\n",
    "            flags=\"t\",\n",
    "            separator=\"comma\",\n",
    "        )\n",
    "        # While we could use .strip().splitlines()[-1].split(\",\") here,\n",
    "        # using a proper CSV reader is more robust.\n",
    "        reader = csv.DictReader(io.StringIO(table_data))\n",
    "        for row in reader:\n",
    "            for key, value in row.items():\n",
    "                new_key = f\"{raster}_{key}\"\n",
    "                results[new_key] = value\n",
    "    print(results)\n",
    "\n",
    "\n",
    "def main():\n",
    "    options, flags = gs.parser()\n",
    "    coords = options[\"coordinates\"].split(\",\")\n",
    "    sample_continuous = options[\"input\"].split(\",\")\n",
    "    viewshed(\n",
    "        elevation=options[\"input\"],\n",
    "        # We will be storing the coordinate values, so we convert them to floats\n",
    "        # instead of just passing them to r.viewshed as a string (string would\n",
    "        # work just fine for that).\n",
    "        coords=(float(coords[0]), float(coords[1])),\n",
    "        output=options[\"output\"],\n",
    "        sample_continuous=sample_continuous,\n",
    "    )\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we will make the script executable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!chmod u+x viewscape.py\n",
    "!grass ~/grassdata/nc_spm_08_grass7/foss4g --exec g.region raster=elevation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grass ~/grassdata/nc_spm_08_grass7/foss4g --exec ./viewscape.py elevation=elevation coordinates=639889.984967892,226698.08042942305 output=\"data.txt\" input=elevation --o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the Data\n",
    "\n",
    "To view the data in the notebook we will render the raster using _grass.jupyter.Map_. Usually, we would just create a GRASS session in the notebook. However, to keep our development environment as is, we will avoid creating a session in the notebook process, but we will use a subprocess to do the rendering into a PNG image. The following uses the `%%python` magic to execute Python code in a subprocess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%python\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "import grass.script.setup  # Needed only in 8.2 and older.\n",
    "\n",
    "with grass.script.setup.init(\"~/grassdata/nc_spm_08_grass7/foss4g\") as session:\n",
    "    gs.run_command(\"r.colors\", map=\"viewshed\", color=\"grey\")\n",
    "    ortho_map = gj.Map(use_region=True)\n",
    "    ortho_map.d_rast(map=\"elevation\")\n",
    "    ortho_map.d_rast(map=\"viewshed\")\n",
    "    # Save the image (in a standard notebook, we would just display the image now).\n",
    "    ortho_map.save(\"result.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, use _Image_ from _IPython.display_ to display the PNG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(\"result.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using r.univar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grass ~/grassdata/nc_spm_08_grass7/foss4g --exec r.univar map=elevation zones=viewshed -t separator=comma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task: Implement JSON Output\n",
    "\n",
    "```python\n",
    "json.dump(...)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task: Run the Tool with Two Rasters for Sampling\n",
    "\n",
    "The following computes NDVI so we have second continuous raster to sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grass ~/grassdata/nc_spm_08_grass7/foss4g --exec i.vi red=lsat7_2002_30 nir=lsat7_2002_40 viname=ndvi output=ndvi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Task: Implement CSV Output Instead of JSON\n",
    "\n",
    "```python\n",
    "csv.DictWriter(...)\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}