Python/pygrass: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(+Interface to listing maps (g.list/g.mlist))
m (updated URL to GRASS 83 manual)
 
(12 intermediate revisions by 3 users not shown)
Line 1: Line 1:
''pygrass'' is a library that extends the GRASS capabilities to allow users to access the low-level GRASS API.
''PyGRASS'' is a library that extends the GRASS GIS capabilities to allow users to access the low-level GRASS API.


== Workshops ==
== Documentation ==


=== XIV Italian meeting of GRASS and GFOSS users ===
Detailed PyGRASS documentation: http://grass.osgeo.org/grass83/manuals/libpython/
 
Two workshops were presented during the [http://geomorfolab.arch.unige.it/genova2013/ XIV Italian meeting of GRASS] at Genova.
The workshops use '''ipython notebook''' to show the python and pygrass API, all the material are available on github ([https://github.com/zarch/workshop-python python], [https://github.com/zarch/workshop-pygrass pygrass]). All the execute examples are reported here:
 
* python
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/00_intro.ipynb      Introduction]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/01_data_types.ipynb Data Types]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/02_syntax.ipynb    Syntax]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/03_function.ipynb  Function]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/04_class.ipynb      Class]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/05_other.ipynb      Other]
 
* pygrass
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/00_Modules.ipynb    Modules]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/01_GIS_objects.ipynb GIS]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/02_Vector.ipynb      Vector]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/03_Raster.ipynb      Raster]


== How to test library ==
== How to test library ==


''pygrass'' has doctest inside its code. You can run doctest inside the [http://grass.osgeo.org/sampledata/north_carolina/nc_basic_spm_grass7.tar.gz North Carolina basic location], using the mapset '''user1'''.  
''pygrass'' has 'doctest' included in its code. You can run doctest in the [http://grass.osgeo.org/sampledata/north_carolina/nc_basic_spm_grass7.tar.gz North Carolina basic location], using the mapset '''user1'''.  


To test the single module (file) you have to move to the source code of pygrass and run the following code (this is an example of functions.py):
To test a single module (file) you have to change into the source code directory of pygrass and run the following code (this is an example of functions.py):


<source lang="bash">
<source lang="bash">
Line 32: Line 15:
</source>
</source>


=== Already tested modules ===
== References ==


==== Working ====
* Zambelli, P., Gebbert, S., Ciolli, M., 2013. ''Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS)''. ISPRS International Journal of Geo-Information 2, 201–219. ([http://dx.doi.org/10.3390/ijgi2010201 DOI] | [http://www.mdpi.com/2220-9964/2/1/201/pdf PDF])


The files are found in GRASS GIS 7, lib/python/pygrass/
* [http://grass.osgeo.org/grass71/manuals/libpython/ pyGrass documentation] (updated weekly from SVN)
* functions.py
* gis/region.py
* gis/__init__.py
* modules/__init.py
* vector/abstract.py
* vector/__init__.py
* vector/basic.py
* vector/table.py
* vector/geometry.py
* raster/abstract.py
* ...


See https://svn.osgeo.org/grass/grass/trunk/lib/python/
== Sample PyGRASS scripts ==


==== Not working ====
=== Training material ===


* ...
Two workshops were presented during the [http://geomorfolab.arch.unige.it/genova2013/ XIV Italian meeting of GRASS] at Genova.
The workshops use '''ipython notebook''' to show the python and pygrass API, all the material are available on github ([https://github.com/zarch/workshop-python python], [https://github.com/zarch/workshop-pygrass pygrass]). All the execute examples are reported here:


== References ==
* python training (Genova)
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/00_intro.ipynb      Introduction]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/01_data_types.ipynb Data Types]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/02_syntax.ipynb    Syntax]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/03_function.ipynb  Function]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/04_class.ipynb      Class]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-python/master/05_other.ipynb      Other]


* Zambelli, P., Gebbert, S., Ciolli, M., 2013. ''Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS)''. ISPRS International Journal of Geo-Information 2, 201–219. ([http://dx.doi.org/10.3390/ijgi2010201 DOI] | [http://www.mdpi.com/2220-9964/2/1/201/pdf PDF])
* python training (EURAC)
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/00_intro.ipynb Introduction]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/01_data_types.ipynb Data Types]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/02_syntax.ipynb Syntax]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/03_version_control_system.ipynb Version Control System (hg)]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/04_objects.ipynb Objects]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/05_python_debugger.ipynb Debugger]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/06_functions.ipynb Functions]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/07_numpy.ipynb Numpy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/08_matplotlib.ipynb Matplotlib]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/09_scipy.ipynb Scipy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/10_pandas.ipynb Pandas]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/11_sympy.ipynb Sympy]
** [http://nbviewer.ipython.org/urls/bitbucket.org/zarch/python-course/raw/teacher/12_user_interfaces.ipynb GUI]


* <strike>[http://www.ing.unitn.it/~zambelli/projects/pygrass/ (old) pyGrass documentation]</strike>
* pygrass training
* [http://grass.osgeo.org/programming7/pygrass/ pyGrass documentation] (updated weekly from SVN)
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/00_Modules.ipynb Modules]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/01_GIS_objects.ipynb GIS]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/02_Vector.ipynb Vector]
** [http://nbviewer.ipython.org/urls/raw.github.com/zarch/workshop-pygrass/master/03_Raster.ipynb Raster]


== Sample PyGRASS scripts ==
=== Interface to listing maps (g.list/g.mlist) ===


=== Interface to listing maps (g.list/g.mlist) ===
Here we call a GRASS module that writes to stdout and do not call a Python function that returns a Python object, therefore we can save stdout and then parse it with:
In this case we are calling a GRASS module that write to the stdout and not a Python function that return a Python object, therefore you can save the stdout and then parse it with:


<source lang="python">
<source lang="python">
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import general as g
import subprocess as sub
import subprocess as sub
gl = g.mlist
 
gl(type='rast', pattern='elev-*', stdout=sub.PIPE)
gl = g.list  # GRASS 7
gl(type='raster', pattern='elev-*', stdout=sub.PIPE)
gl.outputs.stdout.split()
gl.outputs.stdout.split()
['elev-000-000', 'elev-000-001', 'elev-000-002', 'elev-001-000',
['elev-000-000', 'elev-000-001', 'elev-000-002', 'elev-001-000',
Line 89: Line 84:
</source>
</source>


that uses the "G_list" function through ctypes and therefore is faster.
... that uses the "G_list" function through ctypes and therefore is faster.


If you choose the first solution is better if you use directly the
Note: If you choose the first solution it is better if you use directly the
function in [http://grass.osgeo.org/programming7/namespacepython_1_1script_1_1core.html python.script.core.mlist_grouped() etc.].
function in [http://grass.osgeo.org/programming7/namespacepython_1_1script_1_1core.html python.script.core.mlist_grouped() etc.].
=== Interface to copying maps (g.copy) ===
<source lang="python">
from grass.pygrass.modules.shortcuts import general as g
vectinmap = 'zipcodes_wake'
vectoutmap = 'zip_elev_zonal'
g.copy(vector=(vectinmap,vectoutmap))
</source>


=== Sample script for opening, query and closing of a vector map ===
=== Sample script for opening, query and closing of a vector map ===
Line 152: Line 156:
#USER HAS TO SET THOSE
#USER HAS TO SET THOSE
#QUIET REPORTS
#QUIET REPORTS
QIET=True
QUIET=True
#OVERWRITE EXISTING FILES
#OVERWRITE EXISTING FILES
OVR=False
OVR=False
Line 283: Line 287:
print "Prefix IN:\t",pref
print "Prefix IN:\t",pref
print "Prefix OUT:\t", outpref
print "Prefix OUT:\t", outpref
i.landsat_toar(input_prefix=pref,output_prefix=outpref,metfile=metadata[0],sensor=LSENSOR,quiet=QIET,overwrite=OVR)
i.landsat_toar(input_prefix=pref,output_prefix=outpref,metfile=metadata[0],sensor=LSENSOR,quiet=QUIET,overwrite=OVR)
#Atmospheric Correction
#Atmospheric Correction
print "Atmospheric Correction"
print "Atmospheric Correction"
Line 334: Line 338:
print "\t> ",b
print "\t> ",b
print "\t> ",b_out
print "\t> ",b_out
i.atcorr(input=b, elevation="dem", visibility="vis", parameters=prm, output=b_out, flags="ra", range=[0,1],quiet=QIET,overwrite=OVR)
i.atcorr(input=b, elevation="dem", visibility="vis", parameters=prm, output=b_out, flags="ra", range=[0,1],quiet=QUIET,overwrite=OVR)
r.timestamp(map=b_out,date=timestamp,finish_=False)
r.timestamp(map=b_out,date=timestamp,finish_=False)
satbandno = satbandno + 1
satbandno = satbandno + 1
Line 362: Line 366:
b_ndvi=pref[:-2]+".surf.ndvi"
b_ndvi=pref[:-2]+".surf.ndvi"
print "NDVI:\t",b_ndvi
print "NDVI:\t",b_ndvi
i.vi(red=b3,nir=b4,output=b_ndvi,viname="ndvi",quiet=QIET,overwrite=OVR,finish_=False)
i.vi(red=b3,nir=b4,output=b_ndvi,viname="ndvi",quiet=QUIET,overwrite=OVR,finish_=False)


b_in=[b1,b2,b3,b4,b5,b7]
b_in=[b1,b2,b3,b4,b5,b7]
b_albedo=pref[:-2]+".surf.albedo"
b_albedo=pref[:-2]+".surf.albedo"
print "Albedo:\t",b_albedo
print "Albedo:\t",b_albedo
i.albedo(input=b_in,output=b_albedo,flags="lc",quiet=QIET,overwrite=OVR,finish_=False)
i.albedo(input=b_in,output=b_albedo,flags="lc",quiet=QUIET,overwrite=OVR,finish_=False)
b_emissivity=pref[:-2]+".surf.emissivity"
b_emissivity=pref[:-2]+".surf.emissivity"
print "Emissivity:\t",b_emissivity
print "Emissivity:\t",b_emissivity
i.emissivity(input=b_ndvi, output=b_emissivity,quiet=QIET,overwrite=OVR,finish_=False)
i.emissivity(input=b_ndvi, output=b_emissivity,quiet=QUIET,overwrite=OVR,finish_=False)
</source>
 
=== Multiprocessing example: parallelized SHAPE file import ===


<source lang="python">
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Multiprocessing example: parallelized SHAPE file import
Created on Fri, Oct 18, 2013, posted to grass-user@lists.osgeo.org
@author: Pietro Zambelli, freely inspired by: http://stackoverflow.com/a/16071616
        http://lists.osgeo.org/pipermail//grass-user/2013-October/069130.html
"""
# Directory containing SHAPE files to import
DIR = '/data/shp'
###
from multiprocessing import Queue, Process, cpu_count
from os.path import split
from subprocess import Popen
from grass.pygrass.functions import findfiles
def spawn(func):
    def fun(q_in, q_out):
        while True:
            path, cmdstr = q_in.get()
            if path is None:
                break
            q_out.put(func(path, cmdstr))
    return fun
def mltp_importer(dirpath, match, cmdstr, func, nprocs=cpu_count()):
    q_in = Queue(1)
    q_out = Queue()
    procs = [Process(target=spawn(func), args=(q_in, q_out))
            for _ in range(nprocs)]
    for proc in procs:
        proc.daemon = True
        proc.start()
    # set the parameters
    sent = [q_in.put((path, cmdstr)) for path in findfiles(dirpath, match)]
    # set the end of the cycle
    [q_in.put((None, None)) for proc in procs]
    [proc.join() for proc in procs]
    return [q_out.get() for _ in range(len(sent))]
def importer(path, cmdstr):
    name = split(path)[-1][:-4]
    print name
    popen = Popen(cmdstr.format(path=path, name=name), shell=True)
    popen.wait()
    return path, name, False if popen.returncode else True
CMD = 'v.in.ogr dsn={path} layer={name} output={name}'
processed = mltp_importer(DIR, '*.shp', CMD, importer)
# check for errors
errors = [p for p in processed if not p[2]]
if errors:
    # do something (print list of failed SHP files at end)
    pass
</source>
</source>


{{Python}}
{{Python}}

Latest revision as of 21:30, 1 December 2023

PyGRASS is a library that extends the GRASS GIS capabilities to allow users to access the low-level GRASS API.

Documentation

Detailed PyGRASS documentation: http://grass.osgeo.org/grass83/manuals/libpython/

How to test library

pygrass has 'doctest' included in its code. You can run doctest in the North Carolina basic location, using the mapset user1.

To test a single module (file) you have to change into the source code directory of pygrass and run the following code (this is an example of functions.py):

python -m doctest functions.py

References

  • Zambelli, P., Gebbert, S., Ciolli, M., 2013. Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS). ISPRS International Journal of Geo-Information 2, 201–219. (DOI | PDF)

Sample PyGRASS scripts

Training material

Two workshops were presented during the XIV Italian meeting of GRASS at Genova. The workshops use ipython notebook to show the python and pygrass API, all the material are available on github (python, pygrass). All the execute examples are reported here:

Interface to listing maps (g.list/g.mlist)

Here we call a GRASS module that writes to stdout and do not call a Python function that returns a Python object, therefore we can save stdout and then parse it with:

from grass.pygrass.modules.shortcuts import general as g
import subprocess as sub

gl = g.list   # GRASS 7
gl(type='raster', pattern='elev-*', stdout=sub.PIPE)
gl.outputs.stdout.split()
['elev-000-000', 'elev-000-001', 'elev-000-002', 'elev-001-000',
'elev-001-001', 'elev-001-002', 'elev-002-000', 'elev-002-001',
'elev-002-002']

or you can use the "glist" method of the Mapset class:

from grass.pygrass.gis import Mapset
m = Mapset()
m.glist('rast', pattern='elev-*')
['elev-002-000',  'elev-000-000',  'elev-000-002',  'elev-000-001',
'elev-001-001',  'elev-002-002',  'elev-002-001',  'elev-001-000',
'elev-001-002']

... that uses the "G_list" function through ctypes and therefore is faster.

Note: If you choose the first solution it is better if you use directly the function in python.script.core.mlist_grouped() etc..

Interface to copying maps (g.copy)

from grass.pygrass.modules.shortcuts import general as g

vectinmap = 'zipcodes_wake'
vectoutmap = 'zip_elev_zonal'
g.copy(vector=(vectinmap,vectoutmap))

Sample script for opening, query and closing of a vector map

#!/usr/bin/env python
 
# Example for pyGRASS usage - vector API

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.vector import VectorTopo
 
g.message("Assessing vector topology...")

vectmap = 'zipcodes_wake'
zipcodes = VectorTopo(vectmap)

# Open the map with topology:
zipcodes.open()

# query number of topological features
areas   = zipcodes.number_of("areas")
islands = zipcodes.number_of("islands")
print 'Map: <' + vectmap + '> with %d areas and %d islands' % (areas, islands)

# http://www.ing.unitn.it/~zambelli/projects/pygrass/attributes.html
# (note that above documentation is slightly outdated)
dblink = zipcodes.dblinks[0]
print 'DB name:'
print dblink.database
table = dblink.table()
print 'Column names:'
print table.columns.names()
print 'Column types:'
print table.columns.types()

zipcodes.close()

Sample script for Landsat 7 processing

#!/usr/bin/env python
#Date: 7th February, 2013
#Public domain, GRASS GIS

#Run this script in the terminal as:
#python python-grass.py
#For debugging, run as:
#python -i python-grass.py
#PURPOSE
#This script processes LANDSAT 7 ETM+ images
#1 - unzip *.gz files
#2 - import files in GRASS GIS Location of your choice (r.in.gdal)
#3 - DN to Top of Atmosphere reflectance (i.landsat.toar)
#4 - TOA reflectance to Surface reflectance (i.atcorr)
#5 - NDVI (i.vi), Albedo (i.albedo), Emissivity (i.emissivity)

#USER HAS TO SET THOSE
#QUIET REPORTS
QUIET=True
#OVERWRITE EXISTING FILES
OVR=False
#Define Landsat 7 sensor for i.landsat.toar
LSENSOR="tm7"
#Setup the path to the Landsat 7 Directories
rsdatapath="/home/yann/RSDATA/Myanmar/L7"
#Setup your GRASS GIS working directory
gisdb="/home/yann/GRASSDATA"
location="L7_Myanmar"
print location
mapset="PERMANENT"
#DEM input to atmospheric correction
inDEM=rsdatapath+"/dem.tif"
#set L7 Metadata wildcards
wldc_mtl="*_MTL.txt"
#wldc_met="*.met"
#Visibility distance [Km]
vis=18

#Set python path to enable finding of grass.script lib

#END OF USER CHANGES
###DO NOT CHANGE ANYTHING AFTER THIS LINE !
###
#Load necessary libraries
import os, glob, time, re
from grass import script as g
from grass.script import setup as gsetup
gisbase=os.environ['GISBASE']
gsetup.init(gisbase,gisdb,location,mapset)
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules.shortcuts import display as d
#Needed for floor()
from math import *
#env = os.environ.copy()
#env['GRASS_MESSAGE_FORMAT'] = 'gui'
#Function to get a list of L7 Directories in the rsdatapath
def fn(path):
	for top, dirs, files in os.walk(path):
		return [os.path.join(top, dir) for dir in dirs]

#START PROCESS
### PART 0: PRE-PROCESSING STUFF ###
#import DEM for atmospheric correction
#r.in.gdal(input=inDEM,output="dem",overwrite=OVR)
#r.mapcalc(expression="dem=25",overwrite=OVR)
#create a visibility map
r.mapcalc(expression="vis=18",overwrite=OVR)
#Find the central location of the Landsat file from metadata
metadata=[]
fileList=[]
L7Dirs=fn(rsdatapath)
for L7Dir in L7Dirs:
	#Ungzip all of your Landsat7 images in all your directories
#	print "Ungzip Landsat files in\t",L7Dir
#	p=os.system("gzip -d -q "+L7Dir+"/*.gz")
	#Using pthreads on multi-core machines
	#p=os.system("pigz -d "+L7Dir+"/*.gz")
	#Wait ten seconds for gzip to create the tif images
#	time.sleep(10)
	print "Import in GRASS GIS"
	for L7f in glob.glob(os.path.join(L7Dir,"*.[tT][iI][fF]")):
		f1=L7f.replace(L7Dir+"/","")
		f2=f1.replace(".TIF","")
		f3=f2.replace("_B10",".1")
		f4=f3.replace("_B20",".2")
		f5=f4.replace("_B30",".3")
		f6=f5.replace("_B40",".4")
		f7=f6.replace("_B50",".5")
		f8=f7.replace("_B61",".61")
		f9=f8.replace("_B62",".62")
		f10=f9.replace("_B70",".7")
		f11=f10.replace("_B80",".8")
		f12=f11.replace("L72","L71")
		L7r=f12.replace("_VCID_","")
		print "\t> ",L7r
		r.in_gdal(input=L7f,output=L7r,flags="e",overwrite=OVR)
		fileList.append(L7r)

	#reproject the DEM World map for the new PERMANENT location extents
	r.proj(input="dem",location="Myanmar",memory=10000,resolution=90.0,overwrite=OVR)
	#Get list of metadata files
	for metaf in glob.glob(os.path.join(L7Dir,wldc_mtl)):
		metadata.append(metaf)
	print "Metadata in:\n",metadata[0]
	with open(metadata[0],"r") as f:
		data=f.read()

	f.close()
	dt=data.split("\n")
	for idx in range(len(dt)):
		found=dt[idx].find("CORNER_UL_LON_PRODUCT = ")
		if found > 0:
			ulx=float(dt[idx].replace("CORNER_UL_LON_PRODUCT = ",""))
		found=dt[idx].find("CORNER_UL_LAT_PRODUCT = ")
		if found > 0:
			uly=float(dt[idx].replace("CORNER_UL_LAT_PRODUCT = ",""))
		found=dt[idx].find("CORNER_LR_LON_PRODUCT = ")
		if found > 0:
			lrx=float(dt[idx].replace("CORNER_LR_LON_PRODUCT = ",""))
		found=dt[idx].find("CORNER_LR_LAT_PRODUCT = ")
		if found > 0:
			lry=float(dt[idx].replace("CORNER_LR_LAT_PRODUCT = ",""))
		#Two cases, because a string starting by 0 is not converted well
		found=dt[idx].find("SCENE_CENTER_TIME = ")
		if found > 0:
			tim=dt[idx].replace("SCENE_CENTER_TIME = ","")
			print "timestamp=",str(tim)
			gmth=float(tim.split(":")[0])
			gmtm=float(tim.split(":")[1])
			gmtdec=float(tim.split(":")[2].replace("Z",""))
			print gmth,gmtm,gmtdec
		found=dt[idx].find("DATE_ACQUIRED = ")
		if found > 0:
			dat=dt[idx].replace("DATE_ACQUIRED = ","")
			date=dat.split("-")
		found=dt[idx].find("SUN_ELEVATION = ")
                if found > 0:
                        sunza=90-float(dt[idx].replace("SUN_ELEVATION = ",""))


	#L7 DN->Rad->Ref
	print "Convert DN to Rad to TOARef"
	f1=sorted(fileList)[0]
	f2=f1.replace(L7Dir+"/","")
	pref=f2[:-1]
	outpref=pref[:-2]+".toar."
	print "Prefix IN:\t",pref
	print "Prefix OUT:\t", outpref
	i.landsat_toar(input_prefix=pref,output_prefix=outpref,metfile=metadata[0],sensor=LSENSOR,quiet=QUIET,overwrite=OVR)
	#Atmospheric Correction
	print "Atmospheric Correction"
	# Basic script for i.atcorr for L 7 ETM+
	#Geometrical conditions (L7ETM+)
	geom=7
	#Sensor height (satellite is -1000)
	sens_height=-1000
	#Here we suppose you have altitude (DEM) and Visibility (VIS) maps ready
	#---------------------------------------------
	#Visibility dummy value (overwritten by VIS raster input)
	vis=15
	#Altitude dummy value (in Km should be negative in this param file)
	#(overwritten by DEM raster input)
	alt=-1.200
	#datetime of satellite overpass (month, day, GMT decimal hour)
	mdh=str(int(date[1]))+" "+str(int(date[2]))+" "+str("%.2f" % (gmth+gmtm/60.0+gmtdec/3600.0))
	print "MM DD hh.ddd:\t",mdh
	# Central Lat/Long
	Long=ulx+(lrx-ulx)/2.0
	Lat=lry+(uly-lry)/2.0
	print "Center:\t(",Long, ",", Lat,")"
	#Atmospheric mode
	atm_mode=1 #Tropical
	#Aerosol model
	aerosol_mode=2 #Sri Lanka is a small island (maritime)
	#satellite band number (L5TM [25,26,27,28,29,30], L7ETM+ [61,62,63,64,65,66,67])
	satbandno=61 #Band 1 of L7ETM+ is first to undergo atmospheric correction
	#make time stamp for use in time series analysis
	dat1=str(date[2])+"-"+str(date[1])+"-"+str(date[0])
	dat=dat1.replace(" ","")
	timestamp=time.strftime("%d %b %Y",time.strptime(dat,"%d-%m-%Y"))
	print "Timestamp:\t",timestamp
	for idx in [1,2,3,4,5,7]:
		b=pref[:-2]+".toar."+str(idx)
		b_out=b.replace(".toar.",".surf.")
		param=[]
		param.append(str(geom)+" - geometrical conditions=Landsat 7 ETM+\n")
		param.append(str(mdh)+" "+str("%.2f" % Long)+" "+str("%.2f" % Lat)+" - month day hh.ddd longitude latitude (hh.ddd is in decimal hours GMT)\n")
		param.append(str(atm_mode)+" - atmospheric mode=tropical\n")
		param.append(str(aerosol_mode)+" - aerosols model=maritime\n")
		param.append(str(vis)+" - visibility [km] (aerosol model concentration), not used as there is raster input\n")
		param.append(str(alt)+" - mean target elevation above sea level [km] (here 600m asl), not used as there is raster input\n")
		param.append(str(sens_height)+" - sensor height (here, sensor on board a satellite)\n")
		param.append(str(satbandno)+" - i th band of ETM+ Landsat 7 (atcorr internal no)\n")
		f=open(os.path.join(L7Dir,"param_L7.txt"),"w")
		f.writelines(param)
		f.close()
		prm=os.path.join(L7Dir,"param_L7.txt")
		print "\t> ",b
		print "\t> ",b_out
		i.atcorr(input=b, elevation="dem", visibility="vis", parameters=prm, output=b_out, flags="ra", range=[0,1],quiet=QUIET,overwrite=OVR)
		r.timestamp(map=b_out,date=timestamp,finish_=False)
		satbandno = satbandno + 1

	#Allocate surface reflectance names
	b1=pref[:-2]+".surf.1"
	b2=pref[:-2]+".surf.2"
	b3=pref[:-2]+".surf.3"
	b4=pref[:-2]+".surf.4"
	b5=pref[:-2]+".surf.5"
	b61=pref[:-2]+".toar.61"
	b62=pref[:-2]+".toar.62"
	b7=pref[:-2]+".surf.7"
	b8=pref[:-2]+".surf.8"

	### PART 1: BASIC STUFF ###
	b_in=pref[:-2]+".toar."
	b_clouds=pref[:-2]+".toar.acca"
	print "Clouds:\t\t>",b_clouds
	i.landsat_acca(input_prefix=b_in,output=b_clouds,overwrite=OVR)
	#png_clouds=L7Dir+"/"+pref[:-2]+".clouds.png"
	#d.mon(start='png',output=png_clouds,width=800,height=800)
	print "MASK:\t\tON"
	#Should always be rewritten!
	r.mask(raster=b_clouds,flags="i",overwrite=True)

	b_ndvi=pref[:-2]+".surf.ndvi"
	print "NDVI:\t",b_ndvi
	i.vi(red=b3,nir=b4,output=b_ndvi,viname="ndvi",quiet=QUIET,overwrite=OVR,finish_=False)

	b_in=[b1,b2,b3,b4,b5,b7]
	b_albedo=pref[:-2]+".surf.albedo"
	print "Albedo:\t",b_albedo
	i.albedo(input=b_in,output=b_albedo,flags="lc",quiet=QUIET,overwrite=OVR,finish_=False)
	
	b_emissivity=pref[:-2]+".surf.emissivity"
	print "Emissivity:\t",b_emissivity
	i.emissivity(input=b_ndvi, output=b_emissivity,quiet=QUIET,overwrite=OVR,finish_=False)

Multiprocessing example: parallelized SHAPE file import

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Multiprocessing example: parallelized SHAPE file import

Created on Fri, Oct 18, 2013, posted to grass-user@lists.osgeo.org

@author: Pietro Zambelli, freely inspired by: http://stackoverflow.com/a/16071616
         http://lists.osgeo.org/pipermail//grass-user/2013-October/069130.html
"""

# Directory containing SHAPE files to import
DIR = '/data/shp'

###
from multiprocessing import Queue, Process, cpu_count
from os.path import split
from subprocess import Popen

from grass.pygrass.functions import findfiles

def spawn(func):
    def fun(q_in, q_out):
        while True:
            path, cmdstr = q_in.get()
            if path is None:
                break
            q_out.put(func(path, cmdstr))
    return fun


def mltp_importer(dirpath, match, cmdstr, func, nprocs=cpu_count()):
    q_in = Queue(1)
    q_out = Queue()
    procs = [Process(target=spawn(func), args=(q_in, q_out))
             for _ in range(nprocs)]
    for proc in procs:
        proc.daemon = True
        proc.start()

    # set the parameters
    sent = [q_in.put((path, cmdstr)) for path in findfiles(dirpath, match)]
    # set the end of the cycle
    [q_in.put((None, None)) for proc in procs]
    [proc.join() for proc in procs]
    return [q_out.get() for _ in range(len(sent))]


def importer(path, cmdstr):
    name = split(path)[-1][:-4]
    print name
    popen = Popen(cmdstr.format(path=path, name=name), shell=True)
    popen.wait()
    return path, name, False if popen.returncode else True

CMD = 'v.in.ogr dsn={path} layer={name} output={name}'

processed = mltp_importer(DIR, '*.shp', CMD, importer)
# check for errors
errors = [p for p in processed if not p[2]]
if errors:
    # do something (print list of failed SHP files at end)
    pass