04_script_to_grass_module.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:e27262c8d34a63fdbbb3cf4316318128f05b5adbb631f29baf90c0255f1fe5d3"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Changing a script into a GRASS GIS module"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, a Python script must contain the so-called *[shebang](https://en.wikipedia.org/wiki/Shebang_%28Unix%29)* (ie., the definition of the interpreter in the first line of the script) for which we here use `#!/usr/bin/env` to call the right Python interpreter on Unix systems. Second, the code should be written as functions rather than directly on the top level. The standard template is the follows. However, note that the first line here, in IPython Notebook, is `%%python` magic which tells IPython to execute the cell content as Python script. Later, we will rewrite the script call, in the way which is closer to how the script is called in GRASS GIS."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "#!/usr/bin/env python\n",
      "\n",
      "def main():\n",
      "    pass\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The construction `if __name__ == '__main__':` is a standard way in Python of marking code which should be executed when the file is used as a script. It is a best practice to use the above mentioned template. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "GRASS GIS parser"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Every GRASS GIS module must use the GRASS parser mechanism. This very advanced parser helps to check the user interput, format the help text and optionally create a graphical user interface for the new module.\n",
      "In Python, this means calling the `parser()` function from `grass.script` package. This function parses the special comments written at the beginning of the Python file (below the 'shebang'), processes the parameters provided in command line when using the module and provides these data accordingly within the module. These special comments start with `#%` and can be referred to as *interface definition comments*."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Minimal template"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The *interface definition comment* should contain at least the description of the module and two keywords as shown below.  [Existing GRASS GIS Python scripts](https://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/scripts/) may help to understand the best practice. These values are defined in section `module` which contains the `description` and `keyword` keys."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.py\n",
      "#!/usr/bin/env python\n",
      "\n",
      "#%module\n",
      "#% description: Compute and analyze viewsheds\n",
      "#% keyword: raster\n",
      "#% keyword: viewshed\n",
      "#%end\n",
      "\n",
      "import grass.script as gscript\n",
      "\n",
      "def main():\n",
      "    gscript.parser()\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the GRASS GIS parser requires a new line character at the end of file. In this case, we use two empty lines because IPython Notebook will remove one."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To run the script, we need to set its permissions to 'executable'. We do so using `chmod` command."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!chmod u+x r.viewshed.points.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can execute the script as a GRASS GIS module with `--help` to get the command line help output."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!./r.viewshed.points.py --help"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Adding parameters as standard options"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we will add parameters to our module. Let's say, we want to pass the name of an elevation raster map and of a vector points map. In the *interface definition comment* we hence add a section `option` for each of these maps and define respectively standard (or predefined) options. For the elevation raster map we use `G_OPT_R_INPUT` and for the points vector map `G_OPT_V_INPUT`. Both these standard options are named `input` by default, so we have to give them unique names. Here we will use `elevation` for the raster and `points` for the vector map. The option's name is defined using a key named `key`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.py\n",
      "#!/usr/bin/env python\n",
      "\n",
      "#%module\n",
      "#% description: Compute and analyze viewsheds\n",
      "#% keyword: raster\n",
      "#% keyword: viewshed\n",
      "#%end\n",
      "#%option G_OPT_R_INPUT\n",
      "#% key: elevation\n",
      "#%end\n",
      "#%option G_OPT_V_INPUT\n",
      "#% key: points\n",
      "#%end\n",
      "\n",
      "import grass.script as gscript\n",
      "\n",
      "def main():\n",
      "    options, flags = gscript.parser()\n",
      "    elevation = options['elevation']\n",
      "    points = options['points']\n",
      "    print(elevation, points)\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!./r.viewshed.points.py --help"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!./r.viewshed.points.py elevation=elevation points=bridges"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Adding parameters as custom options"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "When using standard options, all items are predefined by the parser but we can modify e.g. `key` or other items to our needs. However, for other cases, we have to define the options completely from scratch. In our example, we introduce a non-standard key with name `max_distance` which has its `type` set to `double`. We use additionally `description` to document the meaning of this option. In case we would need a longer description we may use the `label` key for a short and concise description and the `description` key for more lengthy and wordy description. In addition to the `description` we add also `key_desc` which is a short hint for a proper syntax. Again, the [existing GRASS GIS Python scripts](https://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/scripts/) may help to quickly understand the best practice.\n",
      "\n",
      "Our `max_distance` parameter option has as default value `-1` which stands for *infinity* in the case of viewsheds. The user is not required to provide this option, so we set `required` to `no`. Also, there can be only one value provided for this option, so we set `multiple` to `no` as well."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.py\n",
      "#!/usr/bin/env python\n",
      "\n",
      "#%module\n",
      "#% description: Compute and analyze viewsheds\n",
      "#% keyword: raster\n",
      "#% keyword: viewshed\n",
      "#%end\n",
      "#%option G_OPT_R_INPUT\n",
      "#% key: elevation\n",
      "#%end\n",
      "#%option G_OPT_V_INPUT\n",
      "#% key: points\n",
      "#%end\n",
      "#%option\n",
      "#% key: max_distance\n",
      "#% key_desc: value\n",
      "#% type: double\n",
      "#% description: Maximum visibility radius. By default infinity (-1)\n",
      "#% answer: -1\n",
      "#% multiple: no\n",
      "#% required: no\n",
      "#%end\n",
      "\n",
      "import grass.script as gscript\n",
      "\n",
      "def main():\n",
      "    options, flags = gscript.parser()\n",
      "    elevation = options['elevation']\n",
      "    points = options['points']\n",
      "    max_distance = float(options['max_distance'])\n",
      "    print(elevation, points, max_distance)\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!./r.viewshed.points.py elevation=elevation points=bridges max_distance=50"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Adding parameters as flags"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.py\n",
      "#!/usr/bin/env python\n",
      "\n",
      "#%module\n",
      "#% description: Compute and analyze viewsheds\n",
      "#% keyword: raster\n",
      "#% keyword: viewshed\n",
      "#%end\n",
      "#%option G_OPT_R_INPUT\n",
      "#% key: elevation\n",
      "#%end\n",
      "#%option G_OPT_V_INPUT\n",
      "#% key: points\n",
      "#%end\n",
      "#%option\n",
      "#% key: max_distance\n",
      "#% key_desc: value\n",
      "#% type: double\n",
      "#% description: Maximum visibility radius. By default infinity (-1)\n",
      "#% answer: -1\n",
      "#% multiple: no\n",
      "#% required: no\n",
      "#%end\n",
      "#%flag\n",
      "#% key: c\n",
      "#% description: Consider the curvature of the earth (current ellipsoid)\n",
      "#%end\n",
      "\n",
      "\n",
      "import grass.script as gscript\n",
      "\n",
      "def main():\n",
      "    options, flags = gscript.parser()\n",
      "    elevation = options['elevation']\n",
      "    points = options['points']\n",
      "    max_distance = float(options['max_distance'])\n",
      "    use_curvature = flags['c']\n",
      "    print(elevation, points, max_distance)\n",
      "    if use_curvature:\n",
      "        print(use_curvature, \"is true\")\n",
      "    else:\n",
      "        print(use_curvature, \"is false\")\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!./r.viewshed.points.py -c elevation=elevation points=bridges max_distance=50"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Compilation and formal requirements"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As mentioned before, a GRASS GIS module must use the GRASS parser to handle its command line parameters.\n",
      "A module reads and/or writes geospatial data as GRASS raster or vector maps.\n",
      "A module should not overwrite existing data unless specified by the user using the `--overwrite` flag.\n",
      "For raster and vector maps and other files, the GRASS parser automatically checks their existence and ends the module execution with a proper error message in case the output already exists.\n",
      "\n",
      "Furthermore, a file with documentation which uses [simple HTML syntax](https://trac.osgeo.org/grass/wiki/Submitting/Docs#HTMLPages) must be provided.\n",
      "This documentation is then distributed with the addon and it is also automatically available online ([GRASS addons docs](http://grass.osgeo.org/grass70/manuals/addons/)).\n",
      "A template with the main sections follows."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.html\n",
      "<h2>DESCRIPTION</h2>\n",
      "\n",
      "A long description with details about the method, implementation, usage or whatever is appropriate.\n",
      "\n",
      "<h2>EXAMPLES</h2>\n",
      "\n",
      "Examples of how the module can be used alone or in combination with other modules.\n",
      "Possibly using the GRASS North Carolina State Plane Metric sample Location.\n",
      "At least one screenshot (PNG format) of the result should provided to show the user what to expect.\n",
      "\n",
      "<h2>REFERENCES</h2>\n",
      "\n",
      "Reference or references to paper related to the module or references which algorithm the module is based on.\n",
      "\n",
      "<h2>SEE ALSO</h2>\n",
      "\n",
      "List of related or similar GRASS GIS modules or modules used together with this module as well as any related websites, or\n",
      "related pages at the GRASS GIS User wiki.\n",
      "\n",
      "<h2>AUTHORS</h2>\n",
      "\n",
      "List of author(s), their organizations and funding sources."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All GRASS GIS modules must be distributed under GNU GPL license (version 2 or later).\n",
      "There is a specified way how the first comment in module's main file should look like.\n",
      "Here is how it can look like for our module:\n",
      "\n",
      "    \"\"\"\n",
      "    MODULE:    r.viewshed.points\n",
      "\n",
      "    AUTHOR(S): John Smith <email at some domain>\n",
      "\n",
      "    PURPOSE:   Computes viewshed for multiple points\n",
      "\n",
      "    COPYRIGHT: (C) 2015 John Smith\n",
      "\n",
      "               This program is free software under the GNU General Public\n",
      "               License (>=v2). Read the file COPYING that comes with GRASS\n",
      "               for details.\n",
      "    \"\"\""
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Although Python is not a compiled language like C, we need to 'compile' also the Python based GRASS GIS module in order to include it in our GRASS installation, and to create HTML documentation and GUI. For this a *Makefile* needs to be written which follows a standard template as well. The included 'Script.make' takes care of processing everything, given that the Python script, the HTML documentation and optional screenshot(s) in PNG format are present in the same directory."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file Makefile\n",
      "MODULE_TOPDIR = ../..\n",
      "\n",
      "PGM = r.viewshed.points\n",
      "\n",
      "include $(MODULE_TOPDIR)/include/Make/Script.make\n",
      "\n",
      "default: script"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "No we will compile the module which will also add it to our GRASS GIS installation. To do this, we need to have administrator rights (on Linux, use sudo in command line). First we need to create one directory required for the compilation:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!mkdir `grass70 --config path`/scriptstrings"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then we can compile:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!make MODULE_TOPDIR=`grass70 --config path`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In these two command, we are using 'backticks' syntax to include output of command into another command."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Handling temporary maps and region"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In scripts we have to often create temporary maps, which needs to be removed after the script is finished, or if there is an error. \n",
      "\n",
      "Similarly, when we need to change computational region in a script, we use temporary region, so that we don't affect the current region settings. This also allows running multiple scripts in parallel, each with its own region settings.\n",
      "We will show you how to handle such cases."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Temporary maps"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will show you how to deal with temporary maps on a simple example script. We will generate temporary random maps, each with unique name. Unique name can be achieved by using process ID in combination with some other prefix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "import os\n",
      "import grass.script as gscript\n",
      "\n",
      "def main():\n",
      "    for i in range(5):\n",
      "        name = 'tmp_raster_' + str(os.getpid()) + str(i)\n",
      "        gscript.mapcalc(\"{name} = rand(0, 10)\".format(name=name), seed=1)\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To avoid removing temporary maps manually with [g.remove](http://grass.osgeo.org/grass70/manuals/g.remove.html),"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gscript.run_command('g.remove', type='raster', pattern=\"tmp_raster_*\", flags='f')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "we write a function which removes those temporary maps after at the end of the script. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "import os\n",
      "import grass.script as gscript\n",
      "\n",
      "TMP_MAPS = []\n",
      "\n",
      "def main():\n",
      "    for i in range(5):\n",
      "        name = 'tmp_raster_' + str(os.getpid()) + str(i)\n",
      "        TMP_MAPS.append(name)\n",
      "        gscript.mapcalc(\"{name} = rand(0, 10)\".format(name=name), seed=1)\n",
      "    cleanup()\n",
      "        \n",
      "def cleanup():\n",
      "    gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "However in case of error, the cleanup function won't be called:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "import os\n",
      "import grass.script as gscript\n",
      "\n",
      "TMP_MAPS = []\n",
      "\n",
      "def main():\n",
      "    for i in range(5):\n",
      "        name = 'tmp_raster_' + str(os.getpid()) + str(i)\n",
      "        TMP_MAPS.append(name)\n",
      "        gscript.mapcalc(\"{name} = rand(0, 10)\".format(name=name), seed=1)\n",
      "    # we raise intentionally error here:\n",
      "    raise TypeError\n",
      "    cleanup()\n",
      "\n",
      "def cleanup():\n",
      "    gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gscript.list_pairs(type='raster', pattern='tmp_raster_*')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Therefore we will use Python _atexit_ module to ensure the cleanup functions is called everytime, even when the script ends with error."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "import os\n",
      "import atexit\n",
      "import grass.script as gscript\n",
      "\n",
      "TMP_MAPS = []\n",
      "\n",
      "def main():\n",
      "    for i in range(5):\n",
      "        name = 'tmp_raster_' + str(os.getpid()) + str(i)\n",
      "        TMP_MAPS.append(name)\n",
      "        gscript.mapcalc(\"{name} = rand(0, 10)\".format(name=name), seed=1)\n",
      "    # we raise intentionally error here:\n",
      "    raise TypeError\n",
      "\n",
      "def cleanup():\n",
      "    gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    atexit.register(cleanup)\n",
      "    main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Temporary region"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using temporary region is simple - library function `use_temp_region` is called in the beginning of the script and we can then call [g.region](http://grass.osgeo.org/grass70/manuals/g.region.html) anywhere in the script. First, we set the region to match a raster map."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import grass.script as gscript\n",
      "gscript.run_command('g.region', raster='elevation')\n",
      "print gscript.region()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As shown in this example,\n",
      "the region set in the script doesn't affect the current region set in the previous step."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%python\n",
      "import grass.script as gscript\n",
      "\n",
      "\n",
      "def main():\n",
      "    gscript.run_command('g.region', n=100, s=0, w=0, e=100, res=1)\n",
      "    print gscript.region()\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    gscript.use_temp_region()\n",
      "    main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print gscript.region()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Complete module example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we will convert the script from previous section into a fully-fledged module.\n",
      "\n",
      "The input is a vector map with points and elevation map. The output is a new vector map of vector points with viewshed characteristics saved in attribute table and the raster viewsheds (basename specified in G_OPT_R_BASENAME_OUTPUT). Additionally, module allows to set maximum distance and observer elevation and consider curvature when computing viewshed.\n",
      "\n",
      "__Tip:__ use `r.viewshed --script` to copy the definition of options `max_distance` and `observer_elevation`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file r.viewshed.points.py\n",
      "#!/usr/bin/env python\n",
      "\n",
      "#%module\n",
      "#% description: Compute and analyze viewsheds\n",
      "#% keyword: raster\n",
      "#% keyword: viewshed\n",
      "#%end\n",
      "#%option G_OPT_R_INPUT\n",
      "#% key: elevation\n",
      "#%end\n",
      "#%option G_OPT_V_INPUT\n",
      "#% key: points\n",
      "#%end\n",
      "#%option G_OPT_V_OUTPUT\n",
      "#% key: output_points\n",
      "#%end\n",
      "#%option G_OPT_R_BASENAME_OUTPUT\n",
      "#% key: viewshed_basename\n",
      "#%end\n",
      "#%option\n",
      "#% key: observer_elevation\n",
      "#% type: double\n",
      "#% required: no\n",
      "#% multiple: no\n",
      "#% key_desc: value\n",
      "#% description: Viewing elevation above the ground\n",
      "#% answer: 1.75\n",
      "#%end\n",
      "#%option\n",
      "#% key: max_distance\n",
      "#% key_desc: value\n",
      "#% type: double\n",
      "#% description: Maximum visibility radius. By default infinity (-1)\n",
      "#% answer: -1\n",
      "#% multiple: no\n",
      "#% required: no\n",
      "#%end\n",
      "#%flag\n",
      "#% key: c\n",
      "#% description: Consider the curvature of the earth (current ellipsoid)\n",
      "#%end\n",
      "\n",
      "\n",
      "import os\n",
      "import atexit\n",
      "import grass.script as gscript\n",
      "from grass.pygrass.vector import VectorTopo\n",
      "\n",
      "\n",
      "TMP_MAPS = []\n",
      "\n",
      "\n",
      "def main():\n",
      "    options, flags = gscript.parser()\n",
      "    elevation = options['elevation']\n",
      "    input_points = options['points']\n",
      "    basename = options['viewshed_basename']\n",
      "    output_points = options['output_points']\n",
      "    observer_elev = options['observer_elevation']\n",
      "    max_distance = float(options['max_distance'])\n",
      "    flag_curvature = 'c' if flags['c'] else ''\n",
      "\n",
      "    tmp_prefix = 'tmp_r_viewshed_points_' + str(os.getpid()) + '_'\n",
      "    tmp_viewshed_name = tmp_prefix + 'viewshed'\n",
      "    tmp_viewshed_vector_name = tmp_prefix + 'viewshed_vector'\n",
      "    tmp_visible_points = tmp_prefix + 'points'\n",
      "    tmp_point = tmp_prefix + 'current_point'\n",
      "    TMP_MAPS.extend([tmp_viewshed_name, tmp_viewshed_vector_name, tmp_visible_points, tmp_point])\n",
      "\n",
      "    columns = [('cat', 'INTEGER'),\n",
      "               ('area', 'DOUBLE PRECISION'),\n",
      "               ('n_points_visible', 'INTEGER'),\n",
      "               ('distance_to_closest', 'DOUBLE PRECISION')]\n",
      "\n",
      "    with VectorTopo(input_points, mode='r') as points, \\\n",
      "         VectorTopo(output_points, mode='w', tab_cols=columns) as output:\n",
      "\n",
      "        for point in points:\n",
      "            viewshed_id = str(point.cat)\n",
      "            viewshed_name = basename + '_' + viewshed_id\n",
      "            gscript.run_command('r.viewshed', input=elevation,\n",
      "                                output=tmp_viewshed_name, coordinates=point.coords(),\n",
      "                                max_distance=max_distance, flags=flag_curvature,\n",
      "                                observer_elev=observer_elev, overwrite=True, quiet=True)\n",
      "            gscript.mapcalc(exp=\"{viewshed} = if({tmp}, {vid}, null())\".format(viewshed=viewshed_name,\n",
      "                                                                               tmp=tmp_viewshed_name,\n",
      "                                                                               vid=viewshed_id))\n",
      "\n",
      "            # viewshed size\n",
      "            cells = gscript.parse_command('r.univar', map=viewshed_name,\n",
      "                                          flags='g')['n']\n",
      "            area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']\n",
      "\n",
      "            # visible points\n",
      "            gscript.run_command('r.to.vect', input=viewshed_name, output=tmp_viewshed_vector_name,\n",
      "                                type='area', overwrite=True, quiet=True)\n",
      "            gscript.run_command('v.select', ainput=input_points, atype='point',\n",
      "                                binput=tmp_viewshed_vector_name, btype='area', \n",
      "                                operator='overlap', flags='t', \n",
      "                                output=tmp_visible_points, overwrite=True, quiet=True)\n",
      "            n_points_visible = gscript.vector_info_topo(tmp_visible_points)['points'] - 1\n",
      "\n",
      "            # distance to closest visible point\n",
      "            if float(n_points_visible) >= 1:\n",
      "                gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (point.x, point.y),\n",
      "                                      output=tmp_point, overwrite=True, quiet=True)\n",
      "                distance = gscript.read_command('v.distance', from_=tmp_point, from_type='point', flags='p',\n",
      "                                                to=tmp_visible_points, to_type='point', upload='dist',\n",
      "                                                dmin=1, quiet=True).strip()\n",
      "\n",
      "                distance = float(distance.splitlines()[1].split('|')[1])\n",
      "            else:\n",
      "                distance = 0\n",
      "\n",
      "            #\n",
      "            # write each point with its attributes\n",
      "            #\n",
      "            output.write(point, (area, n_points_visible, distance))\n",
      "            output.table.conn.commit()\n",
      "\n",
      "def cleanup():\n",
      "    gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')\n",
      "    \n",
      "if __name__ == '__main__':\n",
      "    atexit.register(cleanup)\n",
      "    main()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!make MODULE_TOPDIR=`grass70 --config path`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import grass.script as gscript\n",
      "gscript.run_command('r.viewshed.points', elevation='elevation',\n",
      "                    points='input_points', viewshed_basename='viewshed',\n",
      "                    output_points='new_points', flags='c', overwrite=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Submitting the new module to the GRASS GIS addons repository"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is important to follow the __PEP8 style__ guide. You can use the *pep8* tool to control it. To identify problems in your Python code, you can use *pylint*. See also the GRASS GIS [submitting rules and best practices](https://trac.osgeo.org/grass/wiki/Submitting/Python) for Python code. Please note that not all examples here or elsewhere strictly follow all PEP8 rules. This is usually just for saving space or for other documentation-related reasons. Another case when PEP8 cannot be applied are interface definition comments (`#%`) at the beginning of the file.  However in general, you should always respect PEP8 in your scripts.\n",
      "\n",
      "Read about getting a write access to GRASS addons repository at [How To Contribute](https://trac.osgeo.org/grass/wiki/HowToContribute#WriteaccesstotheGRASS-Addons-SVNrepository) and then contact the GRASS Project Steering Committee, the [grass-psc](https://lists.osgeo.org/listinfo/grass-psc) mailing list (you have to subscribed to the list first).\n",
      "\n",
      "GRASS GIS is using a Subversion (SVN) repository provided by OSGeo. To get your code into the repository upon acceptance by the PSC, you need the following commands:\n",
      "* `svn checkout https://svn.osgeo.org/grass/grass-addons` to download locally the repository,\n",
      "* `svn add` to add your files and directories; all files for one module should be in one directory named according to the module, so in our case `r.viewshed.points` directory, put this directory to the proper subdirectory of `grass7` directory according to module's category, e.g. `r.viewshed.points` module shoud go to `grass7/raster`,\n",
      "* set the SVN file properties locally which is easily done with is a script (`tools/module_svn_propset.sh`) found in the `grass-addons` (this script, provided with the files to update, will set the required properties automatically),\n",
      "* `svn commit` to commit a new version and upload it to the remote SVN repository.\n",
      "\n",
      "Now the new module is available to all users and can be installed with [g.extension](http://grass.osgeo.org/grass70/manuals/g.extension.html). Furthermore, it will show up in the subsequent day(s) in the [GRASS GIS 7 addon manual pages](http://grass.osgeo.org/grass70/manuals/addons/)."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}