GRASS and Python: Difference between revisions
m (→Examples) |
|||
Line 69: | Line 69: | ||
See [[GRASS Python Scripting Library]] for notes and examples. | See [[GRASS Python Scripting Library]] for notes and examples. | ||
==Python extensions for GRASS GIS== | ==Python extensions for GRASS GIS== |
Revision as of 18:24, 14 March 2012
(for discussions on the new GRASS GUI, see here)
Python SIGs
Python Special Interest Groups are focused collaborative efforts to develop, improve, or maintain specific Python resources. Each SIG has a charter, a coordinator, a mailing list, and a directory on the Python website. SIG membership is informal, defined by subscription to the SIG's mailing list. Anyone can join a SIG, and participate in the development discussions via the SIG's mailing list. Below is the list of currently active Python SIGs, with links to their resources.
See more at http://www.python.org/community/sigs/
Writing Python scripts in GRASS
Python is a programming language which is more powerful than shell scripting but easier and more forgiving than C. The Python script can contain simple module description definitions which will be processed with g.parser, as shown in the example below. In this way with no extra coding a GUI can be built, inputs checked, and a skeleton help page can be generated automatically. In addition it adds links to the GRASS message translation system. For code which needs access to the power of C, you can access the GRASS C library functions via the Python "ctypes" interface.
- GRASS Python Scripting Library
- http://grass.osgeo.org/programming6/pythonlib.html (GRASS 6)
- http://grass.osgeo.org/programming7/pythonlib.html (GRASS 7)
Code style: Have a look at SUBMITTING_PYTHON.
Creating Python scripts that call GRASS functionality from outside
In order to use GRASS from outside, some environment variables have to be set.
MS-Windows
GISBASE= C:\GRASS-64 GISRC= C:\Documents and Settings\user\.grassrc6 LD_LIBRARY_PATH= C:\GRASS-64\lib PATH= C:\GRASS-64\etc;C:\GRASS-64\etc\python;C:\GRASS-64\lib;C:\GRASS-64\bin;C:\GRASS-64\extralib;C:\GRASS-64\msys\bin;C:\Python26; PYTHONLIB= C:\Python26 PYTHONPATH= C:\GRASS-64\etc\python GRASS_SH= C:\GRASS-64\msys\bin\sh.exe
Some hints:
- The Python interpreter (python.exe) needs to be in the PATH
- Python needs to be associated with the .py extension
- PATHEXT needs to include .py if you want to be able to omit the extension
- PYTHONPATH needs to be set to %WINGISBASE%\etc\python
1-3 should be taken care of by the Python installer. 4 needs to be done by the startup (currently, this doesn't appear to be the case on MS-Windows).
Note:
Currently (as of 22 Feb 2011) if you want to use Python for scripting GRASS on Windows, the best solution is to delete the bundled version of Python 2.5 from the GRASS installation, install Python and the required add-ons (wxPython, NumPy, PyWin32) from their official installers, then edit the GRASS start-up script to remove any references to the bundled version.
Linux
The variables are set like this:
export GISBASE="/usr/local/grass-6.4.svn/" export PATH="$PATH:$GISBASE/bin:$GISBASE/scripts" export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$GISBASE/lib" # for parallel session management, we use process ID (PID) as lock file number: export GIS_LOCK=$$ # path to GRASS settings file export GISRC="$HOME/.grassrc6"
Running external commands from Python
For information on running external commands from Python, see: http://docs.python.org/lib/module-subprocess.html
Avoid using the older os.* functions. Section 17.1.3 lists equivalents using the Popen() interface, which is more robust (particularly on Windows).
Using the GRASS Python Scripting Library
See GRASS Python Scripting Library for notes and examples.
Python extensions for GRASS GIS
wxPython GUI development for GRASS
- See the wxGUI wiki page
GRASS Python Scripting Library
See GRASS Python Scripting Library (Programmer's manual). See also Converting Bash scripts to Python, and sample Python scripts in GRASS 7
Uses for read, feed and pipe, start and exec commands
All of the *_command functions use make_command to construct a command line for a program which uses the GRASS parser. Most of them then pass that command line to subprocess.Popen() via start_command(), except for exec_command() which uses os.execvpe().
[To be precise, they use grass.Popen(), which just calls subprocess.Popen() with shell=True on Windows and shell=False otherwise. On Windows, you need to use shell=True to be able to execute scripts (including batch files); shell=False only works with binary executables.]
start_command() separates the arguments into those which subprocess.Popen() understands and the rest. The rest are passed to make_command() to construct a command line which is passed as the "args" parameter to subprocess.Popen().
In other words, start_command() is a GRASS-oriented interface to subprocess.Popen(). It should be suitable for any situation where you would use subprocess.Popen() to execute a normal GRASS command (one which uses the GRASS parser, which is almost all of them; the main exception is r.mapcalc in 6.x).
Most of the others are convenience wrappers around start_command(), for common use cases.
- run_command() calls the wait() method on the process, so it doesn't return until the command has finished, and returns the command's exit code. Similar to system().
- pipe_command() calls start_command() with stdout=PIPE and returns the process object. You can use the process' .stdout member to read the command's stdout. Similar to popen(..., "r").
- feed_command() calls start_command() with stdin=PIPE and returns the process object. You can use the process' .stdin member to write to the command's stdout. Similar to popen(..., "w")
- read_command() calls pipe_command(), reads the data from the command's stdout, and returns it as a string. Similar to `backticks` in the shell.
- write_command() calls feed_command(), sends the string specified by the "stdin" argument to the command's stdin, waits for the command to finish and returns its exit code. Similar to "echo ... | command".
- parse_command() calls read_command() and parses its output as key-value pairs. Useful for obtaining information from g.region, g.proj, r.info, etc.
- exec_command() doesn't use start_command() but os.execvpe(). This causes the specified command to replace the current program (i.e. the Python script), so exec_command() never returns. Similar to bash's "exec" command. This can be useful if the script is a "wrapper" around a single command, where you construct the command line and execute the command as the final step.
If you have any other questions, you might want to look at the code ($GISBASE/etc/python/grass/script/core.py). Most of these functions are only a few lines long.
Interfacing with NumPy
The grass.script.array module defines a class array which is a subclass of numpy.memmap with .read()
and .write()
methods to read/write the underlying file via r.out.bin/r.in.bin. Metadata can be read with raster::raster_info():
Example:
import grass.script as grass
import grass.script.array as garray
def main():
map = "elevation.dem"
# read map
a = garray.array()
a.read(map)
# get raster map info
print grass.raster_info(map)['datatype']
i = grass.raster_info(map)
# get computational region info
c = grass.region()
print "rows: %d" % c['rows']
print "cols: %d" % c['cols']
# new array for result
b = garray.array()
# calculate new map from input map and store as GRASS raster map
b[...] = (a / 50).astype(int) * 50
b.write("elev.50m")
The size of the array is taken from the current region (computational region).
The main drawback of using numpy is that you're limited by available
memory. Using a subclass of numpy.memmap
lets you use files which may
be much larger, but processing the entire array in one go is likely to
produce in-memory results of a similar size.
Interfacing with NumPy and SciPy
SciPy offers simple access to complex calculations. Example:
from scipy import stats
import grass.script as grass
import grass.script.array as garray
def main():
map = "elevation.dem"
x = garray.array()
x.read(map)
# Descriptive Statistics:
print "max, min, mean, var:"
print x.max(), x.min(), x.mean(), x.var()
print "Skewness test: z-score and 2-sided p-value:"
print stats.skewtest(stats.skew(x))
Interfacing with NumPy, SciPy and Matlab
One may also use the SciPy - Matlab interface:
### SH: in GRASS ### r.out.mat input=elevation output=elev.mat
### PY ###
import scipy.io as sio
# load data
elev = sio.loadmat('elev.mat')
# retrive the actual array. the data set contains also the spatial reference
elev.get('map_data')
data = elev.get('map_data')
# a first simple plot
import pylab
pylab.plot(data)
pylab.show()
# the contour plot
pylab.contour(data)
# obviously data needs to ne reversed
import numpy as np
data_rev = data[::-1]
pylab.contour(data_rev)
# => this is a quick plot. basemap mapping may provide a nicer map!
#######
Testing and installing Python extensions
Debugging
Make sure the script is executable:
chmod +x /path/to/my.extension.py
During development, a Python script can be debugged using the Python Debugger (pdb):
python -m pdb /path/to/my.extension.py input=my_input_layer output=my_output_layer option=value -f
Installation
Once you're happy with your script, you can put it in the scripts/ folder of your GRASS install. To do so, first create a directory named after your extension, then create a Makefile for it, and a HTML man page:
cd /path/to/grass_src/ cd scripts ls # It is useful to check out the existing scripts and their structure mkdir my.extension cd my.extension cp path/to/my.extension.py . touch my.extension.html touch Makefile
Next step is to edit the Makefile. It is a very simple text file, the only thing to check is to put the right extension name (WITHOUT the .py file extension) after PGM:
MODULE_TOPDIR = ../.. PGM = my.extension include $(MODULE_TOPDIR)/include/Make/Script.make default: script
The HTML file would be generated automatically. If you want to add more precisions in it, you can do it (just make sure you start at DESCRIPTION. See existing scripts.)
You can then run "make" within the my.extension folder. Running "make" in the extension directory places the resulting files in the staging directory (path/to/grass_src/dist.<YOUR_ARCH>/). If you're running GRASS from the staging directory (/path/to/grass_src/bin.<YOUR_ARCH>/grass7), subsequent commands will used the updated files.
# in your extension directory (/path/to/grass_src/scripts/my.extension/) make # Starting GRASS from the staging directory /path/to/grass_src/bin.<YOUR_ARCH>/grass7 my.extension help
You can also run "make install" from the top level directory of your GRASS install (say /usr/local/src/grass_trunk/). Running "make install" from the top level just copies the whole of the dist.<YOUR_ARCH>/ directory to the installation directory (e.g. /usr/local/grass70) and the bin.<YOUR_ARCH>/grass70 bin file to the bin directory (e.g. /usr/local/bin), and fixes any embedded paths in scripts and configuration files.
cd /path/to/grass_src make install # Starting GRASS as usual would work and show your extension available grass7 my.extension help
Python Ctypes Interface
This interface allows calling GRASS library functions from Python scripts. See Python Ctypes Examples for details.
Examples:
- Latest and greatest: [GRASS 7 Python scripts]
- More complicated examples <<-- TODO: update to Ctypes
Sample script for GRASS 6 raster access (use within GRASS, Spearfish session):
#!/usr/bin/env python
## TODO: update example to Ctypes
import os, sys
from grass.lib import grass
if "GISBASE" not in os.environ:
print "You must be in GRASS GIS to run this program."
sys.exit(1)
if len(sys.argv)==2:
input = sys.argv[1]
else:
input = raw_input("Raster Map Name? ")
# initialize
grass.G_gisinit('')
# find map in search path
mapset = grass.G_find_cell2(input, '')
# determine the inputmap type (CELL/FCELL/DCELL) */
data_type = grass.G_raster_map_type(input, mapset)
infd = grass.G_open_cell_old(input, mapset)
inrast = grass.G_allocate_raster_buf(data_type)
rown = 0
while True:
myrow = grass.G_get_raster_row(infd, inrast, rown, data_type)
print rown, myrow[0:10]
rown += 1
if rown == 476:
break
grass.G_close_cell(inrast)
grass.G_free(cell)
Sample script for vector access (use within GRASS, Spearfish session):
#!/usr/bin/python
# run within GRASS Spearfish session
# run this before starting python to append module search path:
# export PYTHONPATH=/usr/src/grass70/swig/python
# check with "import sys; sys.path"
# or:
# sys.path.append("/usr/src/grass70/swig/python")
# FIXME: install the grass bindings in $GISBASE/lib/ ?
import os, sys
from grass.lib import grass
from grass.lib import vector as grassvect
if "GISBASE" not in os.environ:
print "You must be in GRASS GIS to run this program."
sys.exit(1)
if len(sys.argv)==2:
input = sys.argv[1]
else:
input = raw_input("Vector Map Name? ")
# initialize
grass.G_gisinit('')
# find map in search path
mapset = grass.G_find_vector2(input,'')
# define map structure
map = grassvect.Map_info()
# define open level (level 2: topology)
grassvect.Vect_set_open_level (2)
# open existing map
grassvect.Vect_open_old(map, input, mapset)
# query
print 'Vect map: ', input
print 'Vect is 3D: ', grassvect.Vect_is_3d (map)
print 'Vect DB links: ', grassvect.Vect_get_num_dblinks(map)
print 'Map Scale: 1:', grassvect.Vect_get_scale(map)
print 'Number of areas:', grassvect.Vect_get_num_areas(map)
# close map
grassvect.Vect_close(map)
Python-GRASS add-ons
Stand-alone addons:
- Jáchym Čepický's G-ps.map, a GUI to typeset printable maps with ps.map (http://193.84.38.2/~jachym/index.py?cat=gpsmap)
- Jáchym Čepický's v.pydigit, a GUI to v.edit (http://les-ejk.cz/?cat=vpydigit)
- Jáchym Čepický's PyWPS, GRASS-Web Processing Service (http://pywps.wald.intevation.org)
Using GRASS gui.tcl in Python
Here is some example code to use the grass automatically generated guis in python code. This could (should) all be bundled up and abstracted away so that the implementation can be replaced later.
import Tkinter
import os
# Startup (once):
tk = Tkinter.Tk()
tk.eval ("wm withdraw .")
tk.eval ("source $env(GISBASE)/etc/gui.tcl")
# Here you could do various things to change what the gui does
# See gui.tcl and README.GUI
# Make a gui (per dialog)
# This sets up a window for the command.
# This can be different to integrate with tkinter:
tk.eval ('set path ".dialog$dlg"')
tk.eval ('toplevel .dialog$dlg')
# Load the code for this command:
fd = os.popen ("d.vect --tcltk")
gui = fd.read()
# Run it
tk.eval(gui)
dlg = tk.eval('set dlg') # This is used later to get and set
# Get the current command in the gui we just made:
currentcommand = tk.eval ("dialog_get_command " + dlg)
# Set the command in the dialog we just made:
tk.eval ("dialog_set_command " + dlg + " {d.vect map=roads}")
FAQ
- Q: Error message "execl() failed: Permission denied" - what to do?
- A: Be sure that the execute bit of the script is set.
Links
General guides
- Wikibook Python Programming
- Quick Python tutorial for programmers of other languages
- More Python tutorials for programmers
- Python programming style guide
- Python Editors
Programming
- Python and GRASS:
- Library interfaces: [GRASS Python Scripting Library http://grass.osgeo.org/programming7/pythonlib.html]
- Graphical user interface (GIU): [GRASS wxPython-based GUI http://grass.osgeo.org/programming7/wxpythonlib.html]
- PyWPS, GRASS-Web Processing Service: WPS
- Python and OSGeo:
- Python and GDAL/OGR:
- Python bindings to PROJ:
- Python and GIS:
- Python and Statistics:
- RPy - Python interface to the R-statistics programming language
- Bindings:
- SIP (C/C++ bindings generator) http://directory.fsf.org/all/Python-SIP.html
- Cython - C-Extensions for Python (compile where speed is needed)
- Other external projects
Presentations
From FOSS4G2006:
- A Python sweeps in the GRASS - A. Frigeri 2006
- GRASS goes web: PyWPS - J. Cepicky 2006 (see also WPS)