First, a Python script must contain the so-called shebang (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.
In [ ]:
#!/usr/bin/env python
def main():
pass
if __name__ == '__main__':
main()
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.
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.
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.
The interface definition comment should contain at least the description of the module and two keywords as shown below. Existing GRASS GIS Python scripts may help to understand the best practice. These values are defined in section module which contains the description and keyword keys.
In [ ]:
%%file r.viewshed.points.py
#!/usr/bin/env python
#%module
#% description: Compute and analyze viewsheds
#% keyword: raster
#% keyword: viewshed
#%end
import grass.script as gscript
def main():
gscript.parser()
if __name__ == '__main__':
main()
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.
To run the script, we need to set its permissions to 'executable'. We do so using chmod command.
In [ ]:
!chmod u+x r.viewshed.points.py
Now we can execute the script as a GRASS GIS module with --help to get the command line help output.
In [ ]:
!./r.viewshed.points.py --help
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.
In [ ]:
%%file r.viewshed.points.py
#!/usr/bin/env python
#%module
#% description: Compute and analyze viewsheds
#% keyword: raster
#% keyword: viewshed
#%end
#%option G_OPT_R_INPUT
#% key: elevation
#%end
#%option G_OPT_V_INPUT
#% key: points
#%end
import grass.script as gscript
def main():
options, flags = gscript.parser()
elevation = options['elevation']
points = options['points']
print(elevation, points)
if __name__ == '__main__':
main()
In [ ]:
!./r.viewshed.points.py --help
In [ ]:
!./r.viewshed.points.py elevation=elevation points=bridges
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 may help to quickly understand the best practice.
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.
In [ ]:
%%file r.viewshed.points.py
#!/usr/bin/env python
#%module
#% description: Compute and analyze viewsheds
#% keyword: raster
#% keyword: viewshed
#%end
#%option G_OPT_R_INPUT
#% key: elevation
#%end
#%option G_OPT_V_INPUT
#% key: points
#%end
#%option
#% key: max_distance
#% key_desc: value
#% type: double
#% description: Maximum visibility radius. By default infinity (-1)
#% answer: -1
#% multiple: no
#% required: no
#%end
import grass.script as gscript
def main():
options, flags = gscript.parser()
elevation = options['elevation']
points = options['points']
max_distance = float(options['max_distance'])
print(elevation, points, max_distance)
if __name__ == '__main__':
main()
In [ ]:
!./r.viewshed.points.py elevation=elevation points=bridges max_distance=50
In [ ]:
%%file r.viewshed.points.py
#!/usr/bin/env python
#%module
#% description: Compute and analyze viewsheds
#% keyword: raster
#% keyword: viewshed
#%end
#%option G_OPT_R_INPUT
#% key: elevation
#%end
#%option G_OPT_V_INPUT
#% key: points
#%end
#%option
#% key: max_distance
#% key_desc: value
#% type: double
#% description: Maximum visibility radius. By default infinity (-1)
#% answer: -1
#% multiple: no
#% required: no
#%end
#%flag
#% key: c
#% description: Consider the curvature of the earth (current ellipsoid)
#%end
import grass.script as gscript
def main():
options, flags = gscript.parser()
elevation = options['elevation']
points = options['points']
max_distance = float(options['max_distance'])
use_curvature = flags['c']
print(elevation, points, max_distance)
if use_curvature:
print(use_curvature, "is true")
else:
print(use_curvature, "is false")
if __name__ == '__main__':
main()
In [ ]:
!./r.viewshed.points.py -c elevation=elevation points=bridges max_distance=50
As mentioned before, a GRASS GIS module must use the GRASS parser to handle its command line parameters.
A module reads and/or writes geospatial data as GRASS raster or vector maps.
A module should not overwrite existing data unless specified by the user using the --overwrite flag.
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.
Furthermore, a file with documentation which uses simple HTML syntax must be provided. This documentation is then distributed with the addon and it is also automatically available online (GRASS addons docs). A template with the main sections follows.
In [ ]:
%%file r.viewshed.points.html
<h2>DESCRIPTION</h2>
A long description with details about the method, implementation, usage or whatever is appropriate.
<h2>EXAMPLES</h2>
Examples of how the module can be used alone or in combination with other modules.
Possibly using the GRASS North Carolina State Plane Metric sample Location.
At least one screenshot (PNG format) of the result should provided to show the user what to expect.
<h2>REFERENCES</h2>
Reference or references to paper related to the module or references which algorithm the module is based on.
<h2>SEE ALSO</h2>
List of related or similar GRASS GIS modules or modules used together with this module as well as any related websites, or
related pages at the GRASS GIS User wiki.
<h2>AUTHORS</h2>
List of author(s), their organizations and funding sources.
All GRASS GIS modules must be distributed under GNU GPL license (version 2 or later). There is a specified way how the first comment in module's main file should look like. Here is how it can look like for our module:
"""
MODULE: r.viewshed.points
AUTHOR(S): John Smith <email at some domain>
PURPOSE: Computes viewshed for multiple points
COPYRIGHT: (C) 2015 John Smith
This program is free software under the GNU General Public
License (>=v2). Read the file COPYING that comes with GRASS
for details.
"""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.
In [ ]:
%%file Makefile
MODULE_TOPDIR = ../..
PGM = r.viewshed.points
include $(MODULE_TOPDIR)/include/Make/Script.make
default: script
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:
In [ ]:
!mkdir `grass70 --config path`/scriptstrings
In [ ]:
!g.gisenv set="MAPSET=PERMANENT"
In [ ]:
!g.gisenv
Then we can compile:
In [ ]:
!make MODULE_TOPDIR=`grass70 --config path`
In these two command, we are using 'backticks' syntax to include output of command into another command.
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.
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. We will show you how to handle such cases.
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:
In [ ]:
%tb
import os
import grass.script as gscript
def main():
for i in range(5):
name = 'tmp_raster_' + str(os.getpid()) + str(i)
gscript.mapcalc("{name} = rand(0, 10)".format(name=name), seed=1)
if __name__ == '__main__':
main()
To avoid removing temporary maps manually with g.remove,
we write a function which removes those temporary maps after at the end of the script.
In [ ]:
%tb
import os
import grass.script as gscript
TMP_MAPS = []
def main():
for i in range(5):
name = 'tmp_raster_' + str(os.getpid()) + str(i)
TMP_MAPS.append(name)
gscript.mapcalc("{name} = rand(0, 10)".format(name=name), seed=1)
cleanup()
def cleanup():
gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')
if __name__ == '__main__':
main()
However in case of error, the cleanup function won't be called:
In [ ]:
import os
import grass.script as gscript
TMP_MAPS = []
def main():
for i in range(5):
name = 'tmp_raster_' + str(os.getpid()) + str(i)
TMP_MAPS.append(name)
gscript.mapcalc("{name} = rand(0, 10)".format(name=name), seed=1)
# we raise intentionally error here:
raise TypeError
cleanup()
def cleanup():
gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')
if __name__ == '__main__':
main()
In [ ]:
gscript.list_pairs(type='raster', pattern='tmp_raster_*')
Therefore we will use Python atexit module to ensure the cleanup functions is called everytime, even when the script ends with error.
In [ ]:
import os
import atexit
import grass.script as gscript
TMP_MAPS = []
def main():
for i in range(5):
name = 'tmp_raster_' + str(os.getpid()) + str(i)
TMP_MAPS.append(name)
gscript.mapcalc("{name} = rand(0, 10)".format(name=name), seed=1)
# we raise intentionally error here:
raise TypeError
def cleanup():
gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')
if __name__ == '__main__':
atexit.register(cleanup)
main()
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 anywhere in the script. First, we set the region to match a raster map.
In [ ]:
import grass.script as gscript
gscript.run_command('g.region', raster='elevation')
print gscript.region()
As shown in this example, the region set in the script doesn't affect the current region set in the previous step.
In [ ]:
import grass.script as gscript
def main():
gscript.run_command('g.region', n=100, s=0, w=0, e=100, res=1)
print gscript.region()
if __name__ == '__main__':
gscript.use_temp_region()
main()
In [ ]:
print gscript.region()
Now we will convert the script from previous section into a fully-fledged module.
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.
Tip: use r.viewshed --script to copy the definition of options max_distance and observer_elevation.
In [ ]:
%%file r.viewshed.points.py
#!/usr/bin/env python
#%module
#% description: Compute and analyze viewsheds
#% keyword: raster
#% keyword: viewshed
#%end
#%option G_OPT_R_INPUT
#% key: elevation
#%end
#%option G_OPT_V_INPUT
#% key: points
#%end
#%option G_OPT_V_OUTPUT
#% key: output_points
#%end
#%option G_OPT_R_BASENAME_OUTPUT
#% key: viewshed_basename
#%end
#%option
#% key: observer_elevation
#% type: double
#% required: no
#% multiple: no
#% key_desc: value
#% description: Viewing elevation above the ground
#% answer: 1.75
#%end
#%option
#% key: max_distance
#% key_desc: value
#% type: double
#% description: Maximum visibility radius. By default infinity (-1)
#% answer: -1
#% multiple: no
#% required: no
#%end
#%flag
#% key: c
#% description: Consider the curvature of the earth (current ellipsoid)
#%end
import os
import atexit
import grass.script as gscript
from grass.pygrass.vector import VectorTopo
TMP_MAPS = []
def main():
options, flags = gscript.parser()
elevation = options['elevation']
input_points = options['points']
basename = options['viewshed_basename']
output_points = options['output_points']
observer_elev = options['observer_elevation']
max_distance = float(options['max_distance'])
flag_curvature = 'c' if flags['c'] else ''
tmp_prefix = 'tmp_r_viewshed_points_' + str(os.getpid()) + '_'
tmp_viewshed_name = tmp_prefix + 'viewshed'
tmp_viewshed_vector_name = tmp_prefix + 'viewshed_vector'
tmp_visible_points = tmp_prefix + 'points'
tmp_point = tmp_prefix + 'current_point'
TMP_MAPS.extend([tmp_viewshed_name, tmp_viewshed_vector_name, tmp_visible_points, tmp_point])
columns = [('cat', 'INTEGER'),
('area', 'DOUBLE PRECISION'),
('n_points_visible', 'INTEGER'),
('distance_to_closest', 'DOUBLE PRECISION')]
with VectorTopo(input_points, mode='r') as points, \
VectorTopo(output_points, mode='w', tab_cols=columns) as output:
for point in points:
viewshed_id = str(point.cat)
viewshed_name = basename + '_' + viewshed_id
gscript.run_command('r.viewshed', input=elevation,
output=tmp_viewshed_name, coordinates=point.coords(),
max_distance=max_distance, flags=flag_curvature,
observer_elev=observer_elev, overwrite=True, quiet=True)
gscript.mapcalc(exp="{viewshed} = if({tmp}, {vid}, null())".format(viewshed=viewshed_name,
tmp=tmp_viewshed_name,
vid=viewshed_id))
# viewshed size
cells = gscript.parse_command('r.univar', map=viewshed_name,
flags='g')['n']
area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']
# visible points
gscript.run_command('r.to.vect', input=viewshed_name, output=tmp_viewshed_vector_name,
type='area', overwrite=True, quiet=True)
gscript.run_command('v.select', ainput=input_points, atype='point',
binput=tmp_viewshed_vector_name, btype='area',
operator='overlap', flags='t',
output=tmp_visible_points, overwrite=True, quiet=True)
n_points_visible = gscript.vector_info_topo(tmp_visible_points)['points'] - 1
# distance to closest visible point
if float(n_points_visible) >= 1:
gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (point.x, point.y),
output=tmp_point, overwrite=True, quiet=True)
distance = gscript.read_command('v.distance', from_=tmp_point, from_type='point', flags='p',
to=tmp_visible_points, to_type='point', upload='dist',
dmin=1, quiet=True).strip()
distance = float(distance.splitlines()[1].split('|')[1])
else:
distance = 0
#
# write each point with its attributes
#
output.write(point, (area, n_points_visible, distance))
output.table.conn.commit()
def cleanup():
gscript.run_command('g.remove', type='raster', name=','.join(TMP_MAPS), flags='f')
if __name__ == '__main__':
atexit.register(cleanup)
main()
In [ ]:
# to sbuild the module we need to remove the anaconda python from the PATH
# which has a newer version of GDAL in conflict with the system one used by GRASS
PATH='/usr/local/bin:/usr/local/grass-7.0.2svn/bin:/usr/local/grass-7.0.2svn/scripts:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin'
In [ ]:
import os
In [ ]:
os.environ['PATH']=PATH
In [ ]:
!make MODULE_TOPDIR=`grass70 --config path`
In [ ]:
!v.info input_points
In [ ]:
!g.region vect=input_points res=100 -ap
In [ ]:
import grass.script as gscript
gscript.run_command('r.viewshed.points', elevation='elevation',
points='input_points', viewshed_basename='viewshed',
output_points='new_points', flags='c', overwrite=True)
In [ ]:
!v.info new_points
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 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.
Read about getting a write access to GRASS addons repository at How To Contribute and then contact the GRASS Project Steering Committee, the grass-psc mailing list (you have to subscribed to the list first).
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:
svn checkout https://svn.osgeo.org/grass/grass-addons to download locally the repository,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,tools/module_svn_propset.sh) found in the grass-addons (this script, provided with the files to update, will set the required properties automatically),svn commit to commit a new version and upload it to the remote SVN repository.Now the new module is available to all users and can be installed with g.extension. Furthermore, it will show up in the subsequent day(s) in the GRASS GIS 7 addon manual pages.