02_pygrass_library.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:cf7043d8fd83cc6eda3d0ace5f17bb17e1c970cbccd82d5e7260ea108df81150"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "PyGRASS"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "PyGRASS is a library originally developed during the Google Summer of Code 2012. The [PyGRASS](http://grass.osgeo.org/grass70/manuals/libpython/pygrass_index.html) library adds two main functionalities:\n",
      "\n",
      "* a Python interface through the ctypes binding of the [C API](http://grass.osgeo.org/programming7/) of GRASS, to read and write natively GRASS GIS 7 data structures.\n",
      "* a GRASS GIS module interface using objects to check the parameters and execute the respective modules.\n",
      "\n",
      "For further discussion about the implementation ideas and performance are presented in the article:\n",
      "\n",
      "[Zambelli, P.; Gebbert, S.; Ciolli, M.\t Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS). *ISPRS Int. J. Geo-Inf.* 2013, 2, 201-219.](http://www.mdpi.com/2220-9964/2/1/201)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Playing with the C API of GRASS GIS through Python"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Standard scripting with GRASS modules in Python may sometime seem discouraging especially when you have to do conceptually simple things like: iterate over vector features or raster rows/columns. With the standard scripting approach you would have to use several workaround such as:\n",
      "- create temporary files,\n",
      "- read and parse these files,\n",
      "- work with the data,\n",
      "- and remove all unnecessary files that were created by the script.\n",
      "\n",
      "Most of GRASS GIS is written in C programming language. Using the [C API](http://grass.osgeo.org/programming7/), all this is much more simple since you can directly work on GRASS GIS data and do just what you need to do. However, you perhaps want to stick to Python. Here, the PyGRASS library introduced several objects that allow to interact directly with the data using the underlying C API of GRASS GIS.\n",
      "Three main areas are covered so far:\n",
      "- vector maps,\n",
      "- raster maps,\n",
      "- mapsets, and region.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Working with vector data"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Creating a new vector map"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Import the necessary classes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from grass.pygrass.vector import VectorTopo\n",
      "from grass.pygrass.vector.geometry import Point"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create an instance of the vector map, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view_points = VectorTopo('view_points')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Open the map in write mode, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view_points.open(mode='w')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create some vector geometry features, like two points:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "point1 = Point(635818.8, 221342.4)\n",
      "point2 = Point(633627.7, 227050.2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Add the above two points to the new vector map with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view_points.write(point1)\n",
      "view_points.write(point2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally close the vector map, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view_points.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now do the same thing using the [context manager syntax](https://www.python.org/dev/peps/pep-0343/) provided by `with`, and creating and setting also the attribute table, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Define the columns of the new vector map\n",
      "cols = [(u'cat',       'INTEGER PRIMARY KEY'),\n",
      "        (u'name',      'TEXT')]\n",
      "\n",
      "with VectorTopo('view_points', mode='w', tab_cols=cols, overwrite=True) as view_points:\n",
      "    # save the point and the attribute\n",
      "    view_points.write(point1, ('pub', ))\n",
      "    view_points.write(point2, ('restaurant', ))\n",
      "    # save the changes to the database\n",
      "    view_points.table.conn.commit()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note: we don't need to close the vector map because it is already closed by the context manager."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Reading an existing vector map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with VectorTopo('view_points', mode='r') as points:\n",
      "    for point in points:\n",
      "        print(point, point.attrs['name'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Working with raster data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A class with similar behaviour has been created for raster maps, import this class with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from grass.pygrass.raster import RasterRow"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Reading an existing raster map"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Instanziate a raster class with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elev = RasterRow('elevation', mapset='PERMANENT')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Similar to the vector approach shown above, we can open a raster map with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elev.open(mode='r')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we have row by row access to the whole raster map. Below we show an example to access the first row:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "row = elev[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "row"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since row instances inherit from the numpy arrays, you can treat them as numpy arrays:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "row * 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "row.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To convert the whole raster map to a numpy array, and then back to a GRASS GIS raster map, we can define two functions like:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "from grass.pygrass.raster.buffer import Buffer\n",
      "from grass.pygrass.gis.region import Region\n",
      "\n",
      "def raster2numpy(rastname, mapset=''):\n",
      "    \"\"\"Return a numpy array from a raster map\"\"\"\n",
      "    with RasterRow(rastname, mapset=mapset, mode='r') as rast:\n",
      "        return np.array(rast)\n",
      "\n",
      "\n",
      "def numpy2raster(array, mtype, rastname, overwrite=False):\n",
      "    \"\"\"Save a numpy array to a raster map\"\"\"\n",
      "    reg = Region()\n",
      "    if (reg.rows, reg.cols) != array.shape:\n",
      "        msg = \"Region and array are different: %r != %r\"\n",
      "        raise TypeError(msg % ((reg.rows, reg.cols), array.shape))\n",
      "    with RasterRow(rastname, mode='w', mtype=mtype, overwrite=overwrite) as new:\n",
      "        newrow = Buffer((array.shape[1],), mtype=mtype)\n",
      "        for row in array:\n",
      "            newrow[:] = row[:]\n",
      "            new.put_row(newrow)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elev = raster2numpy('elevation', mapset='PERMANENT')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elev"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "numpy2raster(elev * 2, mtype='FCELL', rastname='doubled_elev', overwrite=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# import a function to render the results\n",
      "from render import view\n",
      "\n",
      "from grass.pygrass.modules.shortcuts import raster as r\n",
      "\n",
      "# set the colors of the input and output maps to a predefined elevation color table\n",
      "r.colors(map='elevation@PERMANENT,doubled_elev', color='elevation')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view(['elevation@PERMANENT'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "view(['doubled_elev'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Calling GRASS GIS modules"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Import the Module class from the pygrass library, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from grass.pygrass.modules import Module"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Execute a GRASS GIS module with (here: viewshed calculation):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed = Module('r.viewshed', input='elevation', output='viewshed2', coordinates=(635818.8, 221342.4), flags='c')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can play with the python instance asking for inputs and outputs parameters and flags"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.inputs.input"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.inputs.coordinates"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.outputs.output"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.flags.c"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can change the parameters values, parameterwise with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.inputs.coordinates = [(635825, 221350)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Or re-execute the module passing only the parameter(s) that need(s) to be different:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed(output='viewshed_newpoint')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can retrieve the command string that will be execute in both bash and Python style:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.get_bash()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "viewshed.get_python()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Extra parameters of the Module class"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "run\\_"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `run_` parameter defines if the module must be execute or not, the default value is True."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "finish\\_"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `finish_` parameter defines if the module has to wait the end of the process, default is True"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "env\\_"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `env_` parameter defines the environmental variable that will be set and used by the process."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "stdin\\_, stdout\\_, stderr\\_"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The parameters `stdin_`, `stdout_` and `stderr_` are used to pass or capture the input or the output of the GRASS GIS modules."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Example"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import subprocess as sub\n",
      "\n",
      "# Another way to use the Module class objects is through the shortcuts\n",
      "# import the shortcuts for the general modules\n",
      "from grass.pygrass.modules.shortcuts import general as g"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can call and capture the stdout of a GRASS GIS modules, with:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "greg = g.region(flags='p', stdout_=sub.PIPE)\n",
      "print(greg.outputs.stdout)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}